[R-sig-Geo] raser::boxplot : cannot add grouping layer

[Agustin Lobo]

According to the help page of raster::boxplot
x: Raster* object
y: if x is a RasterLayer object, y can be an additional RasterLayer to
group the values of x by 'zone'

but I get an error when I add a raster layer that represents a classification:
r1 <- r2 <- r3 <- raster(ncol=10, nrow=10)
r1[] <- rnorm(ncell(r1), 100, 40)
r2[] <- rnorm(ncell(r1), 80, 10)
r3[] <- rnorm(ncell(r1), 120, 30)
s <- stack(r1, r2, r3)
names(s) <- c('A', 'B', 'C')

rc <- round(r1[[1]]/100)
hist(rc)
summary(rc)

boxplot(s[[1]],rc)

Error in as.data.frame.default(data) :
  cannot coerce class "structure("RasterLayer", package = "raster")"
to a data.frame

Am I not interpreting what this function is supposed to do? I want to
get the boxplot
of the values of s by class in rc
(i can do it having rc as an sp pol data frame and using extract,  but
thought raster::boxplot
would be a more direct way)

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raser::boxplot : cannot add grouping layer

[Agustin Lobo]

Solved just by upgrading raster.
Thanks
Agus

On Thu, Oct 8, 2015 at 9:57 PM, Robert J. Hijmans  wrote:
> Hi Agus,
> Your example works for me without error.
> Can you update.packages() and try this in a clean workspace? (perhaps
> another package is interfering?)
> Robert
>
> On Wed, Oct 7, 2015 at 2:44 AM, Agustin Lobo  wrote:
>> According to the help page of raster::boxplot
>> x: Raster* object
>> y: if x is a RasterLayer object, y can be an additional RasterLayer to
>> group the values of x by 'zone'
>>
>> but I get an error when I add a raster layer that represents a 
>> classification:
>> r1 <- r2 <- r3 <- raster(ncol=10, nrow=10)
>> r1[] <- rnorm(ncell(r1), 100, 40)
>> r2[] <- rnorm(ncell(r1), 80, 10)
>> r3[] <- rnorm(ncell(r1), 120, 30)
>> s <- stack(r1, r2, r3)
>> names(s) <- c('A', 'B', 'C')
>>
>> rc <- round(r1[[1]]/100)
>> hist(rc)
>> summary(rc)
>>
>> boxplot(s[[1]],rc)
>>
>> Error in as.data.frame.default(data) :
>>   cannot coerce class "structure("RasterLayer", package = "raster")"
>> to a data.frame
>>
>> Am I not interpreting what this function is supposed to do? I want to
>> get the boxplot
>> of the values of s by class in rc
>> (i can do it having rc as an sp pol data frame and using extract,  but
>> thought raster::boxplot
>> would be a more direct way)
>>
>> Thanks
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Apply a different threshold for each layer in a brick object

[Agustin Lobo]

Hi!

Given
r <- raster(ncol=10, nrow=10)
r[]=1:ncell(r)
s <- brick(r,r,r,r,r,r)
s <- s * 1:6
summary(s)

applying a common threshold is simple, i.e.
th <- 50
s[s>th] <- th

but what if th is a vector with a different value for each layer?
Is there a faster (raster!) way (i.e. using calc() or stackApply() )
for the following:

s2 <- brick(r,r,r,r,r,r)
s2 <- s2 * 1:6
th <- c(50*(1:6))
for (i in 1:nlayers(s)){
  a <- s2[[i]]
  a[a>th[i]] <- th[i]
  s2[[i]]  <- a
}

Tried with calc() and stackApply(), but get errors.
Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] title in plot of raster::stack or brick objects

[Agustin Lobo]

Given
r <- raster(nrows=10, ncols=10)
r <- setValues(r, 1:ncell(r))
s <- stack(r, sqrt(r), r/sqrt(r))
names(s) <- paste0("S",1:3)
plot(s)

Is there any way to put a main title centered on the multiple plot and
get the names of the layers displayed as well?

plot(s, main="S")
puts the title for the 1st layer only and eliminates the names of the
layers from the plot

Thanks

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] title in plot of raster::stack or brick objects

[Agustin Lobo]

Thanks. Not great though, as the actual position depends on
the shape of the of the display . I.e., if you are on RStudio, open
the zoom window and the main title gets lower than the names of the
layers.

Perhaps the problem would be alleviated if raster::plot() could display
layer names as subtitles (at the bottom of the layers). That would leave more
room for the main title at the main top.

Agus

On Fri, Nov 20, 2015 at 10:42 AM, Jose  wrote:
> Hello,
> try with title():
>
> r <- raster(nrows=10, ncols=10)
> r <- setValues(r, 1:ncell(r))
> s <- stack(r, sqrt(r), r/sqrt(r))
> names(s) <- paste0("S",1:3)
> plot(s)
> title("Title", line=3)
>
> Best
>
> __
> José Manuel Prieto
> https://jmprietob.shinyapps.io/eltiempo/
> es.linkedin.com/in/josemanuelprietoblazquez/
>
> 2015-11-20 10:34 GMT+01:00 Agustin Lobo :
>
>> Given
>> r <- raster(nrows=10, ncols=10)
>> r <- setValues(r, 1:ncell(r))
>> s <- stack(r, sqrt(r), r/sqrt(r))
>> names(s) <- paste0("S",1:3)
>> plot(s)
>>
>> Is there any way to put a main title centered on the multiple plot and
>> get the names of the layers displayed as well?
>>
>> plot(s, main="S")
>> puts the title for the 1st layer only and eliminates the names of the
>> layers from the plot
>>
>> Thanks
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Different values from raster::Moran() and spdep::moran.test()

[Agustin Lobo]

Using the testGrid example provided by
https://stat.ethz.ch/pipermail/r-sig-geo/2014-December/022106.html

lng <- rep(seq(0, 7, by=1), 8)
counter = 1
subCounter = 0
startNum = 0
lat = NULL
while (counter < 65) {
  if (subCounter == 8) {
startNum = startNum + 1
subCounter = 0
  }
  lat = c(lat, startNum)
  subCounter = subCounter + 1
  counter = counter + 1
}
# now add the black/ white chessboard pattern
chess <- rep(c(0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0), 4)
gridNum <- seq(1:64)
testGrid <- data.frame(gridNum, lat, lng, chess)

I define
require(raster)
rtestGrid <- rasterFromXYZ(testGrid[,2:4])
plot(rtestGrid)

and run

Moran(rtestGrid)
which results on:

[1] 0.07438017

while

require(spdep)
nb8q <- cell2nb(8, 8, type="queen", torus=FALSE)
moran.test(testGrid$chess, nb2listw(nb8q, style="B"),alternative="two.sided")

which results on:
Moran's I test under randomisation

data:  testGrid$chess
weights: nb2listw(nb8q, style = "B")

Moran I statistic standard deviate = -0.7714, p-value = 0.4404
alternative hypothesis: two.sided
sample estimates:
Moran I statistic   Expectation  Variance
 -0.06667  -0.015873016   0.004335237

Considering the grid is regular, a negative I value makes sense.

Why is Moran() so different?

Thanks

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Results from spdep::moran.mc() different from raster::Moran() on randomized raster layers

[Agustin Lobo]

I compare results of moran.test() and moran.mc()
to results of randomizing a raster layer and
calculating the mean and range of the raster::Moran() values
using w=matrix(c(1,1,1,1,0,1,1,1,1), 3, 3) (following advise
https://stat.ethz.ch/pipermail/r-sig-geo/2015-November/023761.html)

While the spdep functions let reject the null hypothesis
of CSR because of regularity (which makes sense
because the test is a check board pattern),
the range of raster::Moran() for the randomized raster
layers include the actual value of the observed test layer.
Could this be caused by the fact that the test raster layer
is too small (8x8)? Because sample() is actually not
randomizing as boot()? An error on my side?

This is what I do:

Using the testGrid example provided by
https://stat.ethz.ch/pipermail/r-sig-geo/2014-December/022106.html

lng <- rep(seq(0, 7, by=1), 8)
counter = 1
subCounter = 0
startNum = 0
lat = NULL
while (counter < 65) {
 if (subCounter == 8) {
   startNum = startNum + 1
   subCounter = 0
 }
 lat = c(lat, startNum)
 subCounter = subCounter + 1
 counter = counter + 1
}
# now add the black/ white chessboard pattern
chess <- rep(c(0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0), 4)
gridNum <- seq(1:64)
testGrid <- data.frame(gridNum, lat, lng, chess)


nb8q <- cell2nb(8, 8, type="queen", torus=FALSE)
moran.test(testGrid$chess, nb2listw(nb8q, style="B"), alternative="less")
moran.mc(testGrid$chess,nb2listw(nb8q, style="B"),nsim=200,alternative="less")
wn=matrix(c(1,1,1,1,0,1,1,1,1), 3, 3)
Moran(rtestGrid, w=wn)
nsim=200
rtest <- rtestGrid
vmoran <- (1:nsim)*0
for (i in 1:nsim){
  rtest[] <- sample(rtestGrid,ncell(rtestGrid),replace=TRUE)
  vmoran[i] <- Moran(rtest, w=wn)
  }
mean(vmoran)
var(vmoran)
range(vmoran)
hist(vmoran)
Moran(rtestGrid, w=wn)

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Problems with mapmisc::insetMap()

[Agustin Lobo]

Hi!

When I make an insetMap for an small country, I get a very
uggly "X" marking its position:

require(rgdal)
require(raster)
require(mapmisc)
nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
map.new(nicabg, 0.7,mar=c(2, 2, 2, 0) + 0.1)
plot(nicabg,axes=TRUE)
plot(nica,add=TRUE)
insetMap(nica,"bottomleft", width=0.3)

is there a way to set another symbol or at least making it smaller?

Also, despite the help page says that
"Additional arguments passed to legend" are accepted, I do not find
any to actually work.
For example, cannot set the exact values of x,y for positioning the
inset map. And cannot use locate() to put it wherever I want to.

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Problems with mapmisc::insetMap()

[Agustin Lobo]

Thanks,

Do you know why your solution does not work on a raster layer?
i.e.

require(rgdal)
require(raster)
require(mapmisc)
require(TeachingDemos)

nicabg2 <- openmap(nica,path="cartodb")
plot(nicabg2)
subplot({
plot(world, col=1, xlim=c(-120, 30), ylim=c(-30, 40))
box()
points(coordinates(nica), pch=20, cex=2, col='red')
}, 'bottomright', inset=c(0.01, 0.01), size=c(1, 1))

Agus

On Thu, Feb 25, 2016 at 12:36 AM, John Baumgartner  wrote:
> On Thu, Feb 25, 2016 at 4:26 AM, Agustin Lobo  wrote:
>>
>> Hi!
>>
>> When I make an insetMap for an small country, I get a very
>> uggly "X" marking its position:
>>
>> require(rgdal)
>> require(raster)
>> require(mapmisc)
>> nica <- getData("GADM", country="NIC", level=0)
>> nicabg <- openmap(nica, path="landscape")
>> map.new(nicabg, 0.7,mar=c(2, 2, 2, 0) + 0.1)
>> plot(nicabg,axes=TRUE)
>> plot(nica,add=TRUE)
>> insetMap(nica,"bottomleft", width=0.3)
>>
>> is there a way to set another symbol or at least making it smaller?
>
> Not quite what you asked for, but you could try TeachingDemos::subplots,
> which you might find to be a little more flexible.
>
> E.g.
>
> library(TeachingDemos)
> library(rgdal)
>
> download.file(file.path('http://www.naturalearthdata.com/http/',
> 'www.naturalearthdata.com/download/110m/physical',
> 'ne_110m_land.zip'), f <- tempfile())
> unzip(f, exdir=tempdir())
> world <- readOGR(tempdir(), 'ne_110m_land')
>
> plot(nica, col='gray90', border='black', axes=TRUE, xlim=c(-88, -81.5))
> subplot({
>   plot(world, col=1, xlim=c(-120, 30), ylim=c(-30, 40))
>   box()
>   points(coordinates(nica), pch=20, cex=2, col='red')
> }, 'bottomright', inset=c(0.01, 0.01), size=c(1, 1))
>
>>
>> Also, despite the help page says that
>> "Additional arguments passed to legend" are accepted, I do not find
>> any to actually work.
>> For example, cannot set the exact values of x,y for positioning the
>> inset map. And cannot use locate() to put it wherever I want to.
>
> The "additional arguments" are for scaleBar, which is also documented
> at ?insetMap.
> Note that there is no ellipsis (...) arg for insetMap.
>
>>
>> Thanks
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] resolution of openmap() raster layers

[Agustin Lobo]

Is there any way to download the raster layers
of openmap() with an increased resolution?
I find the quality of the labels very low,
or am I doing something wrong? i.e.

require(raster)
require(mapmisc)
nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
plot(nicabg)

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Problems with mapmisc::insetMap()

[Agustin Lobo]

Much better now!

nica <- getData("GADM", country="NIC", level=0)
nicabg <- openmap(nica, path="landscape")
map.new(nicabg,1,mar=c(2, 2, 2, 0) + 0.1)
plot(nicabg,axes=TRUE)
plot(nica,add=TRUE)
loc <- insetMap(nica,pos=c(-89.9,8.5), width=0.3, col=NA, lty=0,pch=".")
points(loc)

Just 2 questions:
1. In the example you provide, the resulting mark is too thick. I've
made a thinner version with pch=".", but
it would nice having the option of just drawing the central point and
not the window.
2. Is there a way to calculate the size of the plotting window so that
there is no wasted white space
between the map and the axes? I can do this by hand, but then cannot
read the size.

Thanks

On Thu, Feb 25, 2016 at 6:46 PM, Patrick Brown
 wrote:
> Hello all, my first post here!
>
>
> mapmisc (1.4.8 from r-forge) now allows
>
>
> require(rgdal)
> require(raster)
> require(mapmisc)
> nica <- getData("GADM", country="NIC", level=0)
> nicabg <- openmap(nica, path="landscape")
> map.new(nicabg, 0.7,mar=c(2, 2, 2, 0) + 0.1)
> plot(nicabg,axes=TRUE)
> plot(nica,add=TRUE)
> loc = insetMap(nica,pos=c(-85,12), width=0.3, col=NA, lty=0)
> points(loc)
>
> p
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] odd behaviour of llgridlines() regarding axes labels

[Agustin Lobo]

I see an odd behaviour of llgridlines() regarding axes labels:

require(rgdal)
require(raster)
a <- extent(c(-90,-81.5625,8.407168,16.63619))
a <- as(a, 'SpatialPolygons')
projection(a) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
+towgs84=0,0,0"

No axis labels on x axis:
plot(a,axes=FALSE)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)

#same with:
plot(a,axes=FALSE,mar=c(4, 3, 2, 3)+0.1)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)

#and 2 vertical axes with:
plot(a,axes=TRUE)
llgridlines(a, plotLines=TRUE, plotLabels=TRUE)

Any hint?

Thanks
Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] resolution of openmap() raster layers

[Agustin Lobo]

Stunning!
Can I remove the buttons for saving to a bmp file?
What attribution should be used for publishing?
Agus

On Thu, Feb 25, 2016 at 7:42 PM, Chris Reudenbach
 wrote:
> Hi,
>
> if you just want to map the data, mapview could be an option that among
> others avoid the pixel stretching.
>
> require(mapview)
> require(raster)
> nica <- getData("GADM", country="NIC", level=0)
>
> mapview(nica)
>
> mapview(nica,zcol = "POP2000", color = "#FFA500", lwd= 5, alpha.regions =
> 0.4)
>
>
> cheers Chris
>
>
>
>
> Am 25.02.2016 um 18:49 schrieb Barry Rowlingson:
>>
>> On Thu, Feb 25, 2016 at 5:11 PM, Agustin Lobo 
>> wrote:
>>>
>>> Is there any way to download the raster layers
>>> of openmap() with an increased resolution?
>>> I find the quality of the labels very low,
>>> or am I doing something wrong? i.e.
>>>
>>> require(raster)
>>> require(mapmisc)
>>> nica <- getData("GADM", country="NIC", level=0)
>>> nicabg <- openmap(nica, path="landscape")
>>> plot(nicabg)
>>
>>   Map tiles from OpenStreetMap and other map tile providers are images
>> designed to be shown at a fixed resolution. When you plot them in an R
>> graphics window you could be stretching them so that each pixel in the
>> original maps to 1.273 pixels on your screen. So some kind of
>> interpolation or nearest neighbour replacement has to be done, and
>> this makes text labels look bad. Other line work will look bad too.
>>
>>   If you try and download more map tiles at a higher resolution then
>> you'll find the labels are now way too small, because what you've
>> downloaded are map tiles designed for a higher zoom level on a web
>> browser. Web map browsers have a fixed set of zoom values that
>> correspond to the resolution of the map tiles. With an R window, you
>> are free to choose odd zoom factors that give the ugly behaviour.
>>
>>   If you can resize your R window exactly right then you might get
>> something that looks good!
>>
>>   The alternative is to build a background map yourself from
>> OpenStreetMap *vector* data and some code and some styling. Or use a
>> map tile provider that doesn't have text labels and add them to
>> selected places with R graphics commands. Lines and polygons will
>> still be stretched and a bit "jaggy" but our eyes don't notice this as
>> much as badly scaled text.
>>
>> Barry
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] odd behaviour of llgridlines() regarding axes labels

[Agustin Lobo]

Thanks, that fixes it.
I will try graticules as well.
Agus

On Fri, Feb 26, 2016 at 12:45 PM, Michael Sumner  wrote:
> 1) is ok with xpd clipping turned off
>
> op <- par(xpd = NA)
> plot(a,axes=FALSE)
> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
> par(op)
>
> 2)  same
> op <- par(xpd = NA)
> plot(a,axes=FALSE,mar=c(4, 3, 2, 3)+0.1)
> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
> par(op)
>
> 3)  Consider using graticule package instead? You can control the extents
> independently of the meridian/parallel lines (and optionally provided a
> target projection. )
>
> You have to give it the longs/lats you want, there's no automatic detection
> from an object or plot - but for me that's the price of control over what I
> get.
>
> install.packages("graticule")
> library(graticule)
> plot(a,axes=TRUE)
> library(graticule)
> g <- graticule(seq(-90, -82, by = 2), seq(10, 16, by = 2), xlim =
> par("usr")[1:2], ylim = par("usr")[3:4])
> plot(g, add = TRUE, lty = 2)
>
> Cheers, Mike.
>
>
> On Fri, 26 Feb 2016 at 22:21 Agustin Lobo  wrote:
>
>> I see an odd behaviour of llgridlines() regarding axes labels:
>>
>> require(rgdal)
>> require(raster)
>> a <- extent(c(-90,-81.5625,8.407168,16.63619))
>> a <- as(a, 'SpatialPolygons')
>> projection(a) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
>> +towgs84=0,0,0"
>>
>> No axis labels on x axis:
>> plot(a,axes=FALSE)
>> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
>>
>> #same with:
>> plot(a,axes=FALSE,mar=c(4, 3, 2, 3)+0.1)
>> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
>>
>> #and 2 vertical axes with:
>> plot(a,axes=TRUE)
>> llgridlines(a, plotLines=TRUE, plotLabels=TRUE)
>>
>> Any hint?
>>
>> Thanks
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] resolution of openmap() raster layers

[Agustin Lobo]

I get an error at installing from github on MacOSX10.9.5 that I have reported
to https://github.com/environmentalinformatics-marburg/mapview/issues
and it is being dealt there.

With the cran version mapview_1.0.0, I have problems to select the
type of map, I always get the OSM one. I've tried both
nica <- getData("GADM", country="NIC", level=0)
spplot(mapView(nica["ISO"]),colorkey=FALSE,
   map.type="Thunderforest.Landscape")
as in your example and
spplot(mapView(nica["ISO"],map.types="Thunderforest.Landscape"),
 colorkey=FALSE)

Note that mapView() selects the correct map:
mapView(nica,map.types="Thunderforest.Landscape" )

Also, it seems to me that the mapView() display is always projected
(perhaps Pseudo Mercator epsg:3857 ?) even if the spatial object used
in mapView() has another CRS.
In other words
mapView(nica)
is displayed on Pseudo Mercator even if
projection(nica)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Am I wrong?

Finally, is there a way of including an scale bar? Actually, this
option would be interesting even for interactive display.

Thanks
Agus

On Fri, Feb 26, 2016 at 3:32 PM, Chris Reudenbach
 wrote:
> Agus,
>
> Mapview is using leaflet as engine. Due to this you will have the control
> icons on the map because first of all it is designed for interactive mapping
> within RStudio/R.
>
> I think there are two different approaches to save your maps:
>
> If you want to have a dump of the mapviewobject (but including the graphical
> buttons) you'll find a descriptat stackoverflow
> http://stackoverflow.com/questions/31336898/how-to-save-leaflet-in-rstudio-map-as-png-or-jpg-file
>
>
> You may also use  the spplot function from mapview which is designed for
> basic static mapping and to make usable the adavantages of spplot.
>
> Note even if you are dealing with the mapview map object the spplot function
> uses the Openstreetmap package for retrieving  the background maps  (e.g.
> http://www.inside-r.org/packages/cran/OpenStreetMap/docs/openmap). You can
> use the spplot syntax for designing your maps. Up to now this static
> plotting function is still pretty  basic but you may have a try:
> spplot(mapView(nica["POP2000"]),colorkey=FALSE,  lwd= 15, alpha.regions =
> 0.9,
> map.type="stamen-watercolor" )
>
> I think for using the spplot it is better to install the current stable from
> github:
> library(devtools)
> install_github("environmentalinformatics-marburg/mapview", ref = "master")
>
>
> cheers chris
>
>
>
> Am 26.02.2016 um 12:27 schrieb Agustin Lobo:
>>
>> Stunning!
>> Can I remove the buttons for saving to a bmp file?
>> What attribution should be used for publishing?
>> Agus
>>
>> On Thu, Feb 25, 2016 at 7:42 PM, Chris Reudenbach
>>  wrote:
>>>
>>> Hi,
>>>
>>> if you just want to map the data, mapview could be an option that among
>>> others avoid the pixel stretching.
>>>
>>> require(mapview)
>>> require(raster)
>>> nica <- getData("GADM", country="NIC", level=0)
>>>
>>> mapview(nica)
>>>
>>> mapview(nica,zcol = "POP2000", color = "#FFA500", lwd= 5, alpha.regions =
>>> 0.4)
>>>
>>>
>>> cheers Chris
>>>
>>>
>>>
>>>
>>> Am 25.02.2016 um 18:49 schrieb Barry Rowlingson:
>>>>
>>>> On Thu, Feb 25, 2016 at 5:11 PM, Agustin Lobo 
>>>> wrote:
>>>>>
>>>>> Is there any way to download the raster layers
>>>>> of openmap() with an increased resolution?
>>>>> I find the quality of the labels very low,
>>>>> or am I doing something wrong? i.e.
>>>>>
>>>>> require(raster)
>>>>> require(mapmisc)
>>>>> nica <- getData("GADM", country="NIC", level=0)
>>>>> nicabg <- openmap(nica, path="landscape")
>>>>> plot(nicabg)
>>>>
>>>>Map tiles from OpenStreetMap and other map tile providers are images
>>>> designed to be shown at a fixed resolution. When you plot them in an R
>>>> graphics window you could be stretching them so that each pixel in the
>>>> original maps to 1.273 pixels on your screen. So some kind of
>>>> interpolation or nearest neighbour replacement has to be done, and
>>>> this makes text labels look bad. Other line work will look bad too.
>>>>
>>>>If you try and download more map tiles

Re: [R-sig-Geo] Assign a unique ID number to each patch in a raster

[Agustin Lobo]

I think Nelly points out that package raster lacks a more general
clumping tool (i.e., the equivalent to r.clump in grass) that I've
often missed too.
I'm not sure circumventing the way Ben suggest would work for large
raster layers, for which package raster was designed.
I have a note to myself regarding this
https://gist.github.com/johnbaums/6d155cfca28e02b05ad5
and
www.openforis.org/tools/geospatial-toolkit.html
but have not tried it.

Agus

On Fri, Jan 27, 2017 at 1:19 AM, Ben Tupper  wrote:
> Hi,
>
> I think you are asking for connected component labeling. raster::clump() 
> expects your features to be like islands surrounded by NA or 0 background. 
> There might be more Raster* centric ways to do it, but you could use plain 
> old vanilla image processing using the very good imager package ( 
> https://cran.r-project.org/web/packages/imager/ 
>  ).  The steps below show in 
> color your actual values with the connected component IDs labeled on each 
> cell.  I used the 8-connected labeling but 4-connected is also available.
>
> library(raster)
> library(imager)
>
> set.seed(10)
> nc <- 8
> nr <- 8
> MIN <- 1
> MAX <- 4
> LUT <- c("orange", "green", "darkgoldenrod", "cornflowerblue")
>
> r <- raster(ncols=nc, nrows=nr)
> r[] <- round(runif(ncell(r),MIN, MAX),digits=0)
> cc <- r
>
> img <- imager::as.cimg(as.matrix(r))
> labeled <- imager::label(img, high_connectivity = TRUE)
> cc[] <- as.matrix(labeled)
>
> plot(r, col = LUT)
> txt <- cc[]
> xy <- raster::xyFromCell(cc, 1:raster::ncell(cc))
> text(xy[,1], xy[,2], txt )
>
>
> Cheers,
> Ben
>> On Jan 26, 2017, at 4:48 PM, Nelly Reduan  wrote:
>>
>> Hello,
>> I would like to assign a unique ID number to each patch (a patch is composed 
>> of a set of adjacent cells) as shown in this figure:
>>  <2328.png>
>>
>> I tested the clump function (package raster) by applying an eight adjacent 
>> cells rule but this function mixes all cell values into single patch.
>>
>> Here is a code example to create a raster. The raster contains values ranged 
>> from 1 to 9.
>> r <- raster(ncols=12, nrows=12)
>> r[] <- round(runif(ncell(r),1,9),digits=0)
>> plot(r)
>>
>> Is there a way to assign a unique ID number to grouping cells that form a 
>> patch for each class (i.e., values 1, 2, 3, 4, 5, 6, 7, 8, 9) ?
>>
>> Thanks a lot for your time
>> Nell
>>
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org 
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo 
>> 
> Ben Tupper
> Bigelow Laboratory for Ocean Sciences
> 60 Bigelow Drive, P.O. Box 380
> East Boothbay, Maine 04544
> http://www.bigelow.org
>
>
>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] rasterTopolygons: singleparts often needed

[Agustin Lobo]

Could raster::rasterTopolygons() be modified in future versions
to optionally output singlepart polygons? e.g, in the following
r <-raster(nrow=16, ncol=16)
d <- rep(rep(1:4,rep(4,4)),8)
r[] <- c(d,rev(d))
plot(r)
v <- rasterToPolygons(r,dissolve=TRUE)
plot(v[1,])
dim(v@data)

the new option
v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)
would result on
plot(v[1,])
displaying one single square and
dim(v@data)
would be 8 rows and 1 column

(equivalent to gdal_polygonize)

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] rasterTopolygons: singleparts often needed

[Agustin Lobo]

Yes, this is correct. I refer to integrating disaggregate as an option, as
I think that this is what the user needs in most cases.
Agus

On Sun, Feb 26, 2017 at 1:14 PM, John Baumgartner  wrote:
> Not at my computer, but does sp::disaggregate(v) return your singlepart
> polygons?
>
> On Sun, 26 Feb 2017 at 22:40, Agustin Lobo  wrote:
>
>> Could raster::rasterTopolygons() be modified in future versions
>>
>> to optionally output singlepart polygons? e.g, in the following
>>
>> r <-raster(nrow=16, ncol=16)
>>
>> d <- rep(rep(1:4,rep(4,4)),8)
>>
>> r[] <- c(d,rev(d))
>>
>> plot(r)
>>
>> v <- rasterToPolygons(r,dissolve=TRUE)
>>
>> plot(v[1,])
>>
>> dim(v@data)
>>
>>
>>
>> the new option
>>
>> v <- rasterToPolygons(r,dissolve=TRUE,single=TRUE)
>>
>> would result on
>>
>> plot(v[1,])
>>
>> displaying one single square and
>>
>> dim(v@data)
>>
>> would be 8 rows and 1 column
>>
>>
>>
>> (equivalent to gdal_polygonize)
>>
>>
>>
>> Agus
>>
>>
>>
>> ___
>>
>> R-sig-Geo mailing list
>>
>> R-sig-Geo@r-project.org
>>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Avoid axis labels, ticks and marks in ggRGB() plot

[Agustin Lobo]

In order to save space, I prefer to avoid coordinates in the plot
generated by ggRGB. Is there any way? I've tried:

ggplot() +
ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
  theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank())

but get an empty plot.
I can circumvent the problem as follows, but it is kind of ugly and complicated,
there must be another way:

a <- as(extent(rlogo), 'SpatialPolygons')
y <- fortify(a,region=id)
ggplot(data=y,aes(y=lat, x=long)) +
  ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
  geom_polygon(colour="yellow", fill=NA) +
  guides(fill=FALSE) +
  coord_fixed(ratio=1) +
  theme(
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank())

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Avoid axis labels, ticks and marks in ggRGB() plot

[Agustin Lobo]

Yes! Thanks!
I knew there had to be a better way. For my set of plots I can set
t <- theme(
  axis.text = element_blank(),
  axis.ticks = element_blank(),
  axis.title = element_blank())

and then

ggRGB(rlogo, r=1, g=2, b=3) + t +ggtitle("patatin")

Agus



On Mon, Mar 13, 2017 at 2:05 PM, Benjamin Leutner
 wrote:
> The problem with your first approach is that there is no "data layer" which
> could span the coordinate axes.
> Since ggRGB by default uses an annotation_raster() in order to keep your
> "fill" dimension free for other data,  the coordinate limits are empty,
> hence the plotted raster is lost in the dimensionless void and you get see
> nothing ;-)
>
> I would recommend to start your plot with ggRGB like this:
>
> ggRGB(rlogo, r=1, g=2, b=3) +
>   theme(
> axis.text = element_blank(),
> axis.ticks = element_blank(),
> axis.title = element_blank())
>
>
> If you do want to use it with ggLayer = TRUE, there should be at least one
> "geom" with actual data in it (or you'd have to set ggRGB(...,
> geom_raster=TRUE), but I don't see any benefit in that and it would give you
> a huge legend with one entry per rgb color...)
>
> best,
> Benjamin
>
>
> On 03/13/2017 12:10 PM, Agustin Lobo wrote:
>>
>> In order to save space, I prefer to avoid coordinates in the plot
>> generated by ggRGB. Is there any way? I've tried:
>>
>> ggplot() +
>> ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>>theme(
>>  axis.text = element_blank(),
>>  axis.ticks = element_blank(),
>>  axis.title = element_blank())
>>
>> but get an empty plot.
>> I can circumvent the problem as follows, but it is kind of ugly and
>> complicated,
>> there must be another way:
>>
>> a <- as(extent(rlogo), 'SpatialPolygons')
>> y <- fortify(a,region=id)
>> ggplot(data=y,aes(y=lat, x=long)) +
>>ggRGB(rlogo, r=1, g=2, b=3,ggLayer=TRUE) +
>>geom_polygon(colour="yellow", fill=NA) +
>>guides(fill=FALSE) +
>>coord_fixed(ratio=1) +
>>theme(
>>  axis.text = element_blank(),
>>  axis.ticks = element_blank(),
>>  axis.title = element_blank())
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Benjamin Leutner M.Sc.
> PhD candidate
>
> Department of Remote Sensing
> University of Wuerzburg
> Campus Hubland Nord 86
> 97074 Wuerzburg, Germany
>
> Tel: +49-(0)931-31 89594
> Fax: +49-(0)931-31 89594-0
> Email: benjamin.leut...@uni-wuerzburg.de
> Web: http://www.fernerkundung.uni-wuerzburg.de
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

I get lots of errors at installing rgdal on a recently upgraded machine to
Debian Buster and R3.3.3
There are so many errors that I cannot scroll up enough as to get where
the errors actually start and report.
Is there any way to save the output of install.packages("rgdal") ?

Any other reports on problems installing rgdal on Debian Buster ?

Thanks
Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

Thanks,

As indicated, I do:
$ R CMD INSTALL rgdal_1.2-8.tar.gz 2>&1 | tee rgdal.log

The whole output is here:
https://drive.google.com/open?id=0B_LsZU5NlOztb0ZYQzZZRDQ0WkU

The problem starts at:

configure: PROJ.4 version: 4.8.0
./configure: line 3725:  8390 Segmentation fault  ./proj_conf_test
checking PROJ.4: epsg found and readable... yes
./configure: line 3800:  8399 Segmentation fault  ./proj_conf_test
checking PROJ.4: conus found and readable... yes
configure: Package CPP flags:  -I/usr/include/gdal
configure: Package LIBS:  -L/usr/lib -lgdal -lproj
configure: creating ./config.status
config.status: creating src/Makevars
** libs
g++ -I/usr/share/R/include -DNDEBUG -I/usr/include/gdal
-I"/home/alobo/R/i686-pc-linux-gnu-library/3.3/sp/include"   -fpic  -g
-O2 -fdebug-prefix-map=/build/r-base-yDv1ct/r-base-3.3.3=.
-fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g  -c OGR_write.cpp -o OGR_write.o
In file included from /usr/include/gdal/ogr_feature.h:35:0,
 from /usr/include/gdal/ogrsf_frmts.h:35,
 from OGR_write.cpp:2:
/usr/include/gdal/ogr_featurestyle.h:264:17: warning: override
controls (override/final) only available with -std=c++11 or
-std=gnu++11
 GBool Parse() CPL_OVERRIDE;

I'm totally clueless...

Agus

On Fri, Sep 8, 2017 at 4:00 PM, Roger Bivand  wrote:
> Please use the command line: R CMD INSTALL   and pipe the output 
> (tee) 2>&1 to a file. Please do not post the file, read it to find the root 
> cause. Almost certainly the installed GDAL ha different assumptions about the 
> g++ compiler than R itself, so you may need to install GDAL and R from 
> source. There have been unclear reports of issues with -std=.
>
> Hope this helps,
>
> Roger
>
> Roger Bivand
> Norwegian School of Economics
> Bergen, Norway
>
>
>
> Fra: Agustin Lobo
> Sendt: fredag 8. september, 14.49
> Emne: [R-sig-Geo] Errors at installing rgdal on Debian Buster
> Til: r-sig-geo
>
>
> I get lots of errors at installing rgdal on a recently upgraded machine to 
> Debian Buster and R3.3.3 There are so many errors that I cannot scroll up 
> enough as to get where the errors actually start and report. Is there any way 
> to save the output of install.packages("rgdal") ? Any other reports on 
> problems installing rgdal on Debian Buster ? Thanks Agus 
> ___ R-sig-Geo mailing list 
> R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

Using Synaptic.
I just tested with a live usb with Debian Stretch, and with stretch
rgdal installs with no problems (after installing gdal-bin, gdal-dev
and dependencies) .
Agus

On Tue, Sep 12, 2017 at 10:15 AM, Edzer Pebesma
 wrote:
> How did you install gdal and proj on that machine?
>
> On 12/09/17 09:10, Agustin Lobo wrote:
>> Thanks,
>>
>> As indicated, I do:
>> $ R CMD INSTALL rgdal_1.2-8.tar.gz 2>&1 | tee rgdal.log
>>
>> The whole output is here:
>> https://drive.google.com/open?id=0B_LsZU5NlOztb0ZYQzZZRDQ0WkU
>>
>> The problem starts at:
>>
>> configure: PROJ.4 version: 4.8.0
>> ./configure: line 3725:  8390 Segmentation fault  ./proj_conf_test
>> checking PROJ.4: epsg found and readable... yes
>> ./configure: line 3800:  8399 Segmentation fault  ./proj_conf_test
>> checking PROJ.4: conus found and readable... yes
>> configure: Package CPP flags:  -I/usr/include/gdal
>> configure: Package LIBS:  -L/usr/lib -lgdal -lproj
>> configure: creating ./config.status
>> config.status: creating src/Makevars
>> ** libs
>> g++ -I/usr/share/R/include -DNDEBUG -I/usr/include/gdal
>> -I"/home/alobo/R/i686-pc-linux-gnu-library/3.3/sp/include"   -fpic  -g
>> -O2 -fdebug-prefix-map=/build/r-base-yDv1ct/r-base-3.3.3=.
>> -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
>> -D_FORTIFY_SOURCE=2 -g  -c OGR_write.cpp -o OGR_write.o
>> In file included from /usr/include/gdal/ogr_feature.h:35:0,
>>  from /usr/include/gdal/ogrsf_frmts.h:35,
>>  from OGR_write.cpp:2:
>> /usr/include/gdal/ogr_featurestyle.h:264:17: warning: override
>> controls (override/final) only available with -std=c++11 or
>> -std=gnu++11
>>  GBool Parse() CPL_OVERRIDE;
>>
>> I'm totally clueless...
>>
>> Agus
>>
>> On Fri, Sep 8, 2017 at 4:00 PM, Roger Bivand  wrote:
>>> Please use the command line: R CMD INSTALL   and pipe the output 
>>> (tee) 2>&1 to a file. Please do not post the file, read it to find the root 
>>> cause. Almost certainly the installed GDAL ha different assumptions about 
>>> the g++ compiler than R itself, so you may need to install GDAL and R from 
>>> source. There have been unclear reports of issues with -std=.
>>>
>>> Hope this helps,
>>>
>>> Roger
>>>
>>> Roger Bivand
>>> Norwegian School of Economics
>>> Bergen, Norway
>>>
>>>
>>>
>>> Fra: Agustin Lobo
>>> Sendt: fredag 8. september, 14.49
>>> Emne: [R-sig-Geo] Errors at installing rgdal on Debian Buster
>>> Til: r-sig-geo
>>>
>>>
>>> I get lots of errors at installing rgdal on a recently upgraded machine to 
>>> Debian Buster and R3.3.3 There are so many errors that I cannot scroll up 
>>> enough as to get where the errors actually start and report. Is there any 
>>> way to save the output of install.packages("rgdal") ? Any other reports on 
>>> problems installing rgdal on Debian Buster ? Thanks Agus 
>>> ___ R-sig-Geo mailing list 
>>> R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> ___
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Edzer Pebesma
> Institute for Geoinformatics  (ifgi),  University of Münster
> Heisenbergstraße 2, 48149 Münster, Germany; +49 251 83 33081
> Journal of Statistical Software:   http://www.jstatsoft.org/
> Computers & Geosciences:   http://elsevier.com/locate/cageo/
>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

I assume that from the repos I have defined, this is an standard installation:
deb http://ftp.de.debian.org/debian/ testing main
deb-src http://ftp.de.debian.org/debian/ testing main
deb http://ftp.de.debian.org/debian/ testing-updates main
deb-src http://ftp.de.debian.org/debian/ testing-updates main
deb ftp://security.debian.org/debian-security testing/updates main
deb-src ftp://security.debian.org/debian-security testing/updates main

deb http://deb.debian.org/debian/ testing non-free contrib main
deb-src http://deb.debian.org/debian/ testing non-free contrib main
deb http://deb.debian.org/debian/ testing-updates non-free contrib main
deb-src http://deb.debian.org/debian/ testing-updates non-free contrib main

Also:
alobo@debi:~$ dpkg -s gdal-bin
Package: gdal-bin
Status: install ok installed
Priority: optional
Section: science
Installed-Size: 724
Maintainer: Debian GIS Project 
Architecture: i386
Source: gdal (2.2.1+dfsg-2)
Version: 2.2.1+dfsg-2+b2
Depends: gdal-abi-2-2-1, libc6 (>= 2.15), libgcc1 (>= 1:4.2),
libgdal20 (>= 2.2.0), libstdc++6 (>= 5.2)
Suggests: libgdal-grass, python-gdal
Breaks: gdal-bin (<< 1.10.0-0~)

alobo@debi:~$ dpkg -s libgdal-dev
Package: libgdal-dev
Status: install ok installed
Priority: optional
Section: libdevel
Installed-Size: 32633
Maintainer: Debian GIS Project 
Architecture: i386
Source: gdal (2.2.1+dfsg-2)
Version: 2.2.1+dfsg-2+b2
Depends: libgdal20 (= 2.2.1+dfsg-2+b2), libc6-dev, libarmadillo-dev,
libcurl4-gnutls-dev | libcurl-ssl-dev, libdap-dev, libepsilon-dev,
libexpat1-dev, libfreexl-dev, libfyba-dev, libgeos-dev,
libgeotiff-dev, libgif-dev, libhdf4-alt-dev, libhdf5-dev, libjpeg-dev,
libjson-c-dev, libkml-dev, libltdl-dev, liblzma-dev,
default-libmysqlclient-dev, libnetcdf-dev, libogdi3.2-dev,
libopenjp2-7-dev, libpcre3-dev, libpng-dev, libpoppler-private-dev,
libpq-dev, libproj-dev, libqhull-dev, libspatialite-dev,
libsqlite3-dev, libtiff-dev, liburiparser-dev, libwebp-dev,
libxerces-c-dev, libxml2-dev, unixodbc-dev
Suggests: libgdal-doc

I've contacted pkg-grass-de...@lists.alioth.debian.org

Agus

On Tue, Sep 12, 2017 at 11:37 AM, Roger Bivand  wrote:
> My guess is that your g++ is being used differently to build R and GDAL
> (seen before off-list in August). If you install PROJ, GDAL and R from
> source on your platform, they will (most likely) match. It looks as though R
> was not built with -std=c++11 or -std=gnu++11. On your current platform, you
> almost certainly have a more modern g++, and your R should be built with
> this for packages using C++. Do you know where synaptic is getting its
> binaries?
>
> Roger
>
>
> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>
>> Using Synaptic.
>> I just tested with a live usb with Debian Stretch, and with stretch
>> rgdal installs with no problems (after installing gdal-bin, gdal-dev
>> and dependencies) .
>> Agus
>>
>> On Tue, Sep 12, 2017 at 10:15 AM, Edzer Pebesma
>>  wrote:
>>>
>>> How did you install gdal and proj on that machine?
>>>
>>> On 12/09/17 09:10, Agustin Lobo wrote:
>>>>
>>>> Thanks,
>>>>
>>>> As indicated, I do:
>>>> $ R CMD INSTALL rgdal_1.2-8.tar.gz 2>&1 | tee rgdal.log
>>>>
>>>> The whole output is here:
>>>> https://drive.google.com/open?id=0B_LsZU5NlOztb0ZYQzZZRDQ0WkU
>>>>
>>>> The problem starts at:
>>>>
>>>> configure: PROJ.4 version: 4.8.0
>>>> ./configure: line 3725:  8390 Segmentation fault  ./proj_conf_test
>>>> checking PROJ.4: epsg found and readable... yes
>>>> ./configure: line 3800:  8399 Segmentation fault  ./proj_conf_test
>>>> checking PROJ.4: conus found and readable... yes
>>>> configure: Package CPP flags:  -I/usr/include/gdal
>>>> configure: Package LIBS:  -L/usr/lib -lgdal -lproj
>>>> configure: creating ./config.status
>>>> config.status: creating src/Makevars
>>>> ** libs
>>>> g++ -I/usr/share/R/include -DNDEBUG -I/usr/include/gdal
>>>> -I"/home/alobo/R/i686-pc-linux-gnu-library/3.3/sp/include"   -fpic  -g
>>>> -O2 -fdebug-prefix-map=/build/r-base-yDv1ct/r-base-3.3.3=.
>>>> -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
>>>> -D_FORTIFY_SOURCE=2 -g  -c OGR_write.cpp -o OGR_write.o
>>>> In file included from /usr/include/gdal/ogr_feature.h:35:0,
>>>>  from /usr/include/gdal/ogrsf_frmts.h:35,
>>>>  from OGR_write.cpp:2:
>>>> /usr/include/gdal/ogr_featurestyle.h:264:17: warning: override
>>>> controls (override/final) only available with -std=c++11 or
>>>> -std=gnu++11
>>

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

I just installed r-base through synaptic.
alobo@debi:~$ dpkg -s r-base
Package: r-base
Status: install ok installed
Priority: optional
Section: gnu-r
Installed-Size: 60
Maintainer: Dirk Eddelbuettel 
Architecture: all
Version: 3.3.3-1
Depends: r-base-core (>= 3.3.3-1), r-recommended (= 3.3.3-1)
Recommends: r-base-html, r-doc-html
Suggests: ess, r-doc-info | r-doc-pdf

I think I'll try to downgrade to stretch...

Agus


On Tue, Sep 12, 2017 at 12:28 PM, Roger Bivand  wrote:
> From the entrails, GDAL was built with a CXX11 compiler, but R was very
> possibly not, hence the messages. Where did R come from, and was it itself
> built for the same Debian version?
>
>
> Roger
>
>
> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>
>> I assume that from the repos I have defined, this is an standard
>> installation:
>> deb http://ftp.de.debian.org/debian/ testing main
>> deb-src http://ftp.de.debian.org/debian/ testing main
>> deb http://ftp.de.debian.org/debian/ testing-updates main
>> deb-src http://ftp.de.debian.org/debian/ testing-updates main
>> deb ftp://security.debian.org/debian-security testing/updates main
>> deb-src ftp://security.debian.org/debian-security testing/updates main
>>
>> deb http://deb.debian.org/debian/ testing non-free contrib main
>> deb-src http://deb.debian.org/debian/ testing non-free contrib main
>> deb http://deb.debian.org/debian/ testing-updates non-free contrib main
>> deb-src http://deb.debian.org/debian/ testing-updates non-free contrib
>> main
>>
>> Also:
>> alobo@debi:~$ dpkg -s gdal-bin
>> Package: gdal-bin
>> Status: install ok installed
>> Priority: optional
>> Section: science
>> Installed-Size: 724
>> Maintainer: Debian GIS Project 
>> Architecture: i386
>> Source: gdal (2.2.1+dfsg-2)
>> Version: 2.2.1+dfsg-2+b2
>> Depends: gdal-abi-2-2-1, libc6 (>= 2.15), libgcc1 (>= 1:4.2),
>> libgdal20 (>= 2.2.0), libstdc++6 (>= 5.2)
>> Suggests: libgdal-grass, python-gdal
>> Breaks: gdal-bin (<< 1.10.0-0~)
>>
>> alobo@debi:~$ dpkg -s libgdal-dev
>> Package: libgdal-dev
>> Status: install ok installed
>> Priority: optional
>> Section: libdevel
>> Installed-Size: 32633
>> Maintainer: Debian GIS Project 
>> Architecture: i386
>> Source: gdal (2.2.1+dfsg-2)
>> Version: 2.2.1+dfsg-2+b2
>> Depends: libgdal20 (= 2.2.1+dfsg-2+b2), libc6-dev, libarmadillo-dev,
>> libcurl4-gnutls-dev | libcurl-ssl-dev, libdap-dev, libepsilon-dev,
>> libexpat1-dev, libfreexl-dev, libfyba-dev, libgeos-dev,
>> libgeotiff-dev, libgif-dev, libhdf4-alt-dev, libhdf5-dev, libjpeg-dev,
>> libjson-c-dev, libkml-dev, libltdl-dev, liblzma-dev,
>> default-libmysqlclient-dev, libnetcdf-dev, libogdi3.2-dev,
>> libopenjp2-7-dev, libpcre3-dev, libpng-dev, libpoppler-private-dev,
>> libpq-dev, libproj-dev, libqhull-dev, libspatialite-dev,
>> libsqlite3-dev, libtiff-dev, liburiparser-dev, libwebp-dev,
>> libxerces-c-dev, libxml2-dev, unixodbc-dev
>> Suggests: libgdal-doc
>>
>> I've contacted pkg-grass-de...@lists.alioth.debian.org
>>
>> Agus
>>
>> On Tue, Sep 12, 2017 at 11:37 AM, Roger Bivand 
>> wrote:
>>>
>>> My guess is that your g++ is being used differently to build R and GDAL
>>> (seen before off-list in August). If you install PROJ, GDAL and R from
>>> source on your platform, they will (most likely) match. It looks as
>>> though R
>>> was not built with -std=c++11 or -std=gnu++11. On your current platform,
>>> you
>>> almost certainly have a more modern g++, and your R should be built with
>>> this for packages using C++. Do you know where synaptic is getting its
>>> binaries?
>>>
>>> Roger
>>>
>>>
>>> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>>>
>>>> Using Synaptic.
>>>> I just tested with a live usb with Debian Stretch, and with stretch
>>>> rgdal installs with no problems (after installing gdal-bin, gdal-dev
>>>> and dependencies) .
>>>> Agus
>>>>
>>>> On Tue, Sep 12, 2017 at 10:15 AM, Edzer Pebesma
>>>>  wrote:
>>>>>
>>>>>
>>>>> How did you install gdal and proj on that machine?
>>>>>
>>>>> On 12/09/17 09:10, Agustin Lobo wrote:
>>>>>>
>>>>>>
>>>>>> Thanks,
>>>>>>
>>>>>> As indicated, I do:
>>>>>> $ R CMD INSTALL rgdal_1.2-8.tar.gz 2>&1 | tee rgdal.log
>>>>>>
>>>

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

Yes, r-base in Buster is 3.3.3.1:
https://packages.debian.org/buster/r-base

I'll wait to see if the maintainers answer...

Thanks


On Tue, Sep 12, 2017 at 1:09 PM, Roger Bivand  wrote:
> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>
>> I just installed r-base through synaptic.
>> alobo@debi:~$ dpkg -s r-base
>> Package: r-base
>> Status: install ok installed
>> Priority: optional
>> Section: gnu-r
>> Installed-Size: 60
>> Maintainer: Dirk Eddelbuettel 
>> Architecture: all
>> Version: 3.3.3-1
>> Depends: r-base-core (>= 3.3.3-1), r-recommended (= 3.3.3-1)
>> Recommends: r-base-html, r-doc-html
>> Suggests: ess, r-doc-info | r-doc-pdf
>
>
> This is not current R. An offlist solution for someone else was installing
> PROJ, GDAL, and R from source. Is the R 3.3.3-1 binary for your Debian
> version?
>
> Roger
>
>
>>
>> I think I'll try to downgrade to stretch...
>>
>> Agus
>>
>>
>> On Tue, Sep 12, 2017 at 12:28 PM, Roger Bivand 
>> wrote:
>>>
>>> From the entrails, GDAL was built with a CXX11 compiler, but R was very
>>> possibly not, hence the messages. Where did R come from, and was it
>>> itself
>>> built for the same Debian version?
>>>
>>>
>>> Roger
>>>
>>>
>>> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>>>
>>>> I assume that from the repos I have defined, this is an standard
>>>> installation:
>>>> deb http://ftp.de.debian.org/debian/ testing main
>>>> deb-src http://ftp.de.debian.org/debian/ testing main
>>>> deb http://ftp.de.debian.org/debian/ testing-updates main
>>>> deb-src http://ftp.de.debian.org/debian/ testing-updates main
>>>> deb ftp://security.debian.org/debian-security testing/updates main
>>>> deb-src ftp://security.debian.org/debian-security testing/updates main
>>>>
>>>> deb http://deb.debian.org/debian/ testing non-free contrib main
>>>> deb-src http://deb.debian.org/debian/ testing non-free contrib main
>>>> deb http://deb.debian.org/debian/ testing-updates non-free contrib main
>>>> deb-src http://deb.debian.org/debian/ testing-updates non-free contrib
>>>> main
>>>>
>>>> Also:
>>>> alobo@debi:~$ dpkg -s gdal-bin
>>>> Package: gdal-bin
>>>> Status: install ok installed
>>>> Priority: optional
>>>> Section: science
>>>> Installed-Size: 724
>>>> Maintainer: Debian GIS Project 
>>>> Architecture: i386
>>>> Source: gdal (2.2.1+dfsg-2)
>>>> Version: 2.2.1+dfsg-2+b2
>>>> Depends: gdal-abi-2-2-1, libc6 (>= 2.15), libgcc1 (>= 1:4.2),
>>>> libgdal20 (>= 2.2.0), libstdc++6 (>= 5.2)
>>>> Suggests: libgdal-grass, python-gdal
>>>> Breaks: gdal-bin (<< 1.10.0-0~)
>>>>
>>>> alobo@debi:~$ dpkg -s libgdal-dev
>>>> Package: libgdal-dev
>>>> Status: install ok installed
>>>> Priority: optional
>>>> Section: libdevel
>>>> Installed-Size: 32633
>>>> Maintainer: Debian GIS Project 
>>>> Architecture: i386
>>>> Source: gdal (2.2.1+dfsg-2)
>>>> Version: 2.2.1+dfsg-2+b2
>>>> Depends: libgdal20 (= 2.2.1+dfsg-2+b2), libc6-dev, libarmadillo-dev,
>>>> libcurl4-gnutls-dev | libcurl-ssl-dev, libdap-dev, libepsilon-dev,
>>>> libexpat1-dev, libfreexl-dev, libfyba-dev, libgeos-dev,
>>>> libgeotiff-dev, libgif-dev, libhdf4-alt-dev, libhdf5-dev, libjpeg-dev,
>>>> libjson-c-dev, libkml-dev, libltdl-dev, liblzma-dev,
>>>> default-libmysqlclient-dev, libnetcdf-dev, libogdi3.2-dev,
>>>> libopenjp2-7-dev, libpcre3-dev, libpng-dev, libpoppler-private-dev,
>>>> libpq-dev, libproj-dev, libqhull-dev, libspatialite-dev,
>>>> libsqlite3-dev, libtiff-dev, liburiparser-dev, libwebp-dev,
>>>> libxerces-c-dev, libxml2-dev, unixodbc-dev
>>>> Suggests: libgdal-doc
>>>>
>>>> I've contacted pkg-grass-de...@lists.alioth.debian.org
>>>>
>>>> Agus
>>>>
>>>> On Tue, Sep 12, 2017 at 11:37 AM, Roger Bivand 
>>>> wrote:
>>>>>
>>>>>
>>>>> My guess is that your g++ is being used differently to build R and GDAL
>>>>> (seen before off-list in August). If you install PROJ, GDAL and R from
>>>>> source on your platform, they will (most likely) match. It looks as
>>>&g

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

How do I find out that?
Agus

On Tue, Sep 12, 2017 at 1:39 PM, Roger Bivand  wrote:
> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>
>> Yes, r-base in Buster is 3.3.3.1:
>> https://packages.debian.org/buster/r-base
>>
>> I'll wait to see if the maintainers answer...
>
>
> What is the version of g++ on your buster system - the buster r-base does
> not seem to require a version that supports CXX11 or later?
>
> Roger
>
>
>>
>> Thanks
>>
>>
>> On Tue, Sep 12, 2017 at 1:09 PM, Roger Bivand  wrote:
>>>
>>> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>>>
>>>> I just installed r-base through synaptic.
>>>> alobo@debi:~$ dpkg -s r-base
>>>> Package: r-base
>>>> Status: install ok installed
>>>> Priority: optional
>>>> Section: gnu-r
>>>> Installed-Size: 60
>>>> Maintainer: Dirk Eddelbuettel 
>>>> Architecture: all
>>>> Version: 3.3.3-1
>>>> Depends: r-base-core (>= 3.3.3-1), r-recommended (= 3.3.3-1)
>>>> Recommends: r-base-html, r-doc-html
>>>> Suggests: ess, r-doc-info | r-doc-pdf
>>>
>>>
>>>
>>> This is not current R. An offlist solution for someone else was
>>> installing
>>> PROJ, GDAL, and R from source. Is the R 3.3.3-1 binary for your Debian
>>> version?
>>>
>>> Roger
>>>
>>>
>>>>
>>>> I think I'll try to downgrade to stretch...
>>>>
>>>> Agus
>>>>
>>>>
>>>> On Tue, Sep 12, 2017 at 12:28 PM, Roger Bivand 
>>>> wrote:
>>>>>
>>>>>
>>>>> From the entrails, GDAL was built with a CXX11 compiler, but R was very
>>>>> possibly not, hence the messages. Where did R come from, and was it
>>>>> itself
>>>>> built for the same Debian version?
>>>>>
>>>>>
>>>>> Roger
>>>>>
>>>>>
>>>>> On Tue, 12 Sep 2017, Agustin Lobo wrote:
>>>>>
>>>>>> I assume that from the repos I have defined, this is an standard
>>>>>> installation:
>>>>>> deb http://ftp.de.debian.org/debian/ testing main
>>>>>> deb-src http://ftp.de.debian.org/debian/ testing main
>>>>>> deb http://ftp.de.debian.org/debian/ testing-updates main
>>>>>> deb-src http://ftp.de.debian.org/debian/ testing-updates main
>>>>>> deb ftp://security.debian.org/debian-security testing/updates main
>>>>>> deb-src ftp://security.debian.org/debian-security testing/updates main
>>>>>>
>>>>>> deb http://deb.debian.org/debian/ testing non-free contrib main
>>>>>> deb-src http://deb.debian.org/debian/ testing non-free contrib main
>>>>>> deb http://deb.debian.org/debian/ testing-updates non-free contrib
>>>>>> main
>>>>>> deb-src http://deb.debian.org/debian/ testing-updates non-free contrib
>>>>>> main
>>>>>>
>>>>>> Also:
>>>>>> alobo@debi:~$ dpkg -s gdal-bin
>>>>>> Package: gdal-bin
>>>>>> Status: install ok installed
>>>>>> Priority: optional
>>>>>> Section: science
>>>>>> Installed-Size: 724
>>>>>> Maintainer: Debian GIS Project
>>>>>> 
>>>>>> Architecture: i386
>>>>>> Source: gdal (2.2.1+dfsg-2)
>>>>>> Version: 2.2.1+dfsg-2+b2
>>>>>> Depends: gdal-abi-2-2-1, libc6 (>= 2.15), libgcc1 (>= 1:4.2),
>>>>>> libgdal20 (>= 2.2.0), libstdc++6 (>= 5.2)
>>>>>> Suggests: libgdal-grass, python-gdal
>>>>>> Breaks: gdal-bin (<< 1.10.0-0~)
>>>>>>
>>>>>> alobo@debi:~$ dpkg -s libgdal-dev
>>>>>> Package: libgdal-dev
>>>>>> Status: install ok installed
>>>>>> Priority: optional
>>>>>> Section: libdevel
>>>>>> Installed-Size: 32633
>>>>>> Maintainer: Debian GIS Project
>>>>>> 
>>>>>> Architecture: i386
>>>>>> Source: gdal (2.2.1+dfsg-2)
>>>>>> Version: 2.2.1+dfsg-2+b2
>>>>>> Depends: libgdal20 (= 2.2.1+dfsg-2+b2), libc6-dev, libarmadillo-dev,
>>>>>> libcurl4-gnutls-dev | libcurl-ssl-dev, libdap-dev, libepsilon-dev

[R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

The problem seems to be related to the compiler. I had:
gcc version 5.3.1 20160101 (Debian 5.3.1-5)

I have now:
gcc version 7.2.0 (Debian 7.2.0-4)

All weird previus errors are gone but I still have:
configure: PROJ.4 version: 4.8.0
./configure: line 3725:  8390 Segmentation fault  ./proj_conf_test
checking PROJ.4: epsg found and readable... yes
./configure: line 3800:  8399 Segmentation fault  ./proj_conf_test

and finally

** testing if installed package can be loaded
Fatal error: glibc detected an invalid stdio handle
Aborted
ERROR: loading failed
* removing ‘/home/alobo/R/i686-pc-linux-gnu-library/3.3/rgdal’

Surprisingly,  I do have 4.9 installed, not 4.8:

$ dpkg -s proj-bin
Package: proj-bin
Status: install ok installed
Priority: optional
Section: science
Installed-Size: 125
Maintainer: Debian GIS Project 
Architecture: i386
Source: proj
Version: 4.9.3-2

and
$ proj
Rel. 4.9.3, 15 August 2016
usage: proj [ -bCeEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

and I cannot find proj_conf_test anywhere:
$ sudo find / -name 'proj_conf_test'
gives no result

How is the test ./proj_conf_test done so that I can investigate why it
results into PROJ 4.8 instead of 4.9?

Thanks

Agus



On Tue, Sep 12, 2017 at 1:49 PM, Roger Bivand  wrote:
> $ g++ -v
> ...
> gcc version 7.1.1 20170622 (Red Hat 7.1.1-3) (GCC)
>
> in my case on Fedora 26.
>
> Recent versions should be OK, but < 5 may be problematic. If you have
> upgraded in place but not upgraded the compile trains, installing r-base
> will not force that.
>
>
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

I'm now very confused. It seems to me, in line with what you say, that
I had 2 different
kinds of problems:
- Those derived from upgrading from Debian Stretch to Buster (a
consequence of keeping the system as Debian testing: it is just the
consequence of testing passing from Stretch to Buster as
Stretch becomes Stable) and keeping some older versions of different
components (probably installed by packages outside Debian repos). I
guess these problems have been fixed (although I still do not
understand why the rgdal installation claims I have PROJ4.8, which I
cannot find in my system).
- Those derived from Debian packages having been compiled with
different and incompatible options. This is something I cannot fix
myself.

Therefore:

- I will go on compiling myself as Roger suggests to regain an
operational system.
- I will try to test on a fresh Buster system (or a docker, which I
will have to find out what it is...)
- I will report to r-sig-debian and gdal-dev (please suggest an
alternative gdal linst if
gdal-dev is not appropriate).

I'll keep this list informed.

Thanks!

Agus







On Thu, Sep 14, 2017 at 9:26 AM, Roger Bivand  wrote:
> As Edzer said earlier in this thread (and Agus was quite right to post
> rather than contact me directly), reproducing the problem(s) is hard without
> access to the specific platform (thread started here):
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2017-September/025974.html
>
> My guesses so far are that the gcc build trains used for R, and for the
> GDAL/PROJ.4 components are not the same, and may not be the same as the
> installed gcc build train used to install rgdal from source. This probably
> also affects sf on debian buster (and ubuntu?). It may also affect
> GEOS/rgeos.
>
> The gcc versions do change as platforms are upgraded, so for anything using
> C++, the ABI changes will bite hard, and keeping the compiler versions in
> harmony is crucial. The baseline comparison is to install GDAL and its
> dependencies (including PROJ.4) from source (download and unpack source
> tarball, ./configure && make && make install), then probably R from source,
> finally rgdal. This has to work, but unpicks the debian/ubuntu packaging
> system on which many rely.
>
> Please also consider taking this up on r-sig-debian, with a reproducible
> example, best in a docker container.
>
> A further thought is that the most frequent cause of trouble on
> debian/ubuntu previously were multiple installs of different versions of
> GDAL, pulled in by different downstream packages, but the current reports do
> not suggest this as the cause. Still worth checking, though ...
>
> Roger
>
>
> On Wed, 13 Sep 2017, Roger Bivand wrote:
>
>> On Wed, 13 Sep 2017, Agustin Lobo wrote:
>>
>>> The problem seems to be related to the compiler. I had:
>>> gcc version 5.3.1 20160101 (Debian 5.3.1-5)
>>>
>>> I have now:
>>> gcc version 7.2.0 (Debian 7.2.0-4)
>>>
>>> All weird previus errors are gone but I still have:
>>> configure: PROJ.4 version: 4.8.0
>>> ./configure: line 3725:  8390 Segmentation fault  ./proj_conf_test
>>> checking PROJ.4: epsg found and readable... yes
>>> ./configure: line 3800:  8399 Segmentation fault  ./proj_conf_test
>>>
>>> and finally
>>>
>>> ** testing if installed package can be loaded
>>> Fatal error: glibc detected an invalid stdio handle
>>> Aborted
>>> ERROR: loading failed
>>> * removing ‘/home/alobo/R/i686-pc-linux-gnu-library/3.3/rgdal’
>>>
>>> Surprisingly,  I do have 4.9 installed, not 4.8:
>>>
>>> $ dpkg -s proj-bin
>>> Package: proj-bin
>>> Status: install ok installed
>>> Priority: optional
>>> Section: science
>>> Installed-Size: 125
>>> Maintainer: Debian GIS Project 
>>> Architecture: i386
>>> Source: proj
>>> Version: 4.9.3-2
>>>
>>> and
>>> $ proj
>>> Rel. 4.9.3, 15 August 2016
>>> usage: proj [ -bCeEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]
>>>
>>> and I cannot find proj_conf_test anywhere:
>>> $ sudo find / -name 'proj_conf_test'
>>> gives no result
>>>
>>> How is the test ./proj_conf_test done so that I can investigate why it
>>> results into PROJ 4.8 instead of 4.9?
>>
>>
>> It is defined in rgdal/configure.ac (line 273 ff.):
>>
>> if test ${proj_version} -ge 480; then
>> [cat > proj_conf_test.c <<_EOCONF
>> #include 
>> #include 
>> #if PJ_VERSION == 480
>> FILE *pj_open_lib(projCtx, const char *, const char *);
>> #endif
&g

Re: [R-sig-Geo] Errors at installing rgdal on Debian Buster

[Agustin Lobo]

rgdal compiles fine on my Debian Testing (Buster) system by removing
an older pro_api.h file:

1. As suggested by Roger, I located and removed an older version of
proj_api.h taht had not been removed when I removed the older proj4
version:
$ sudo find / -name proj_api.h
/usr/local/include/proj_api.h
/usr/include/proj_api.h
alobo@debi:~$ sudo rm /usr/local/include/proj_api.h

2. Then the rgdal installation goes smoothly;

$ R CMD INSTALL rgdal_1.2-8.tar.gz
.../...
configure: PROJ.4 version: > 4.8.0
.../...

(see https://www.dropbox.com/s/9i8rk74k0vk95uc/rgdal2.log?dl=0 for the
entire output)

3. Then, in R, rgdal loads with no problems

> require(rgdal)
Loading required package: rgdal
Loading required package: sp
rgdal: version: 1.2-8, (SVN revision 663)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.1, released 2017/06/23
 Path to GDAL shared files: /usr/share/gdal/2.2
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-5
> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: i686-pc-linux-gnu (32-bit)
Running under: Debian GNU/Linux buster/sid

locale:
 [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8
LC_PAPER=en_US.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C
LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_1.2-8 sp_1.2-5

loaded via a namespace (and not attached):
[1] tools_3.3.3 grid_3.3.3  lattice_0.20-35

4. Thus I guess all problems were caused by an old gcc compiler +
leftovers in /usr/local/... directories and that I should
report to gdal-dev and r-sig-debian that there is actually no problem
with versions of R. gfal and proj4 in the Debian testing (Buster)
repos.

Correct?

Agus

On Fri, Sep 15, 2017 at 9:32 AM, Roger Bivand  wrote:
> On Fri, 15 Sep 2017, Agustin Lobo wrote:
>
>> I'm now very confused. It seems to me, in line with what you say, that
>> I had 2 different
>> kinds of problems:
>> - Those derived from upgrading from Debian Stretch to Buster (a
>> consequence of keeping the system as Debian testing: it is just the
>> consequence of testing passing from Stretch to Buster as
>> Stretch becomes Stable) and keeping some older versions of different
>> components (probably installed by packages outside Debian repos). I
>> guess these problems have been fixed (although I still do not
>> understand why the rgdal installation claims I have PROJ4.8, which I
>> cannot find in my system).
>
>
> Please see whether there are multiple proj_api.h files on your system - if
> so, rgdal is picking up the first in your search path, but may link against
> a later *.so, explaining the crashes during the configure run.
>
>> - Those derived from Debian packages having been compiled with
>> different and incompatible options. This is something I cannot fix
>> myself.
>>
>> Therefore:
>>
>> - I will go on compiling myself as Roger suggests to regain an
>> operational system.
>> - I will try to test on a fresh Buster system (or a docker, which I
>> will have to find out what it is...)
>
>
> https://github.com/rocker-org/rocker
> http://dirk.eddelbuettel.com/papers/useR2015_docker.pdf
>
> for example.
>
>> - I will report to r-sig-debian and gdal-dev (please suggest an
>> alternative gdal linst if
>> gdal-dev is not appropriate).
>>
>
> (please keep your heads covered as GEOS and GDAL transition to requiring
>>
>> = C++11 ...) Yes, at this stage sharing information is valuable so that
>
> others can learn from hard-won experience.
>
> Roger
>
>
>> I'll keep this list informed.
>>
>> Thanks!
>>
>> Agus
>>
>>
>>
>>
>>
>>
>>
>> On Thu, Sep 14, 2017 at 9:26 AM, Roger Bivand  wrote:
>>>
>>> As Edzer said earlier in this thread (and Agus was quite right to post
>>> rather than contact me directly), reproducing the problem(s) is hard
>>> without
>>> access to the specific platform (thread started here):
>>>
>>> https://stat.ethz.ch/pipermail/r-sig-geo/2017-September/025974.html
>>>
>>> My guesses so far are that the gcc build trains used for R, and for the
>>> GDAL/PROJ.4 components are not the same, and may not be the same as the
>>> installed gcc build train used to install rgdal from source. This
>>> probably
>>> also affects sf on debian buster (and ubuntu?). It may also af

[R-sig-Geo] temp directory not actually reset by rasterOptions()

[Agustin Lobo]

I have problems with an small /tmp partition
and thought using
rasterOptions(tmpdir="/home/alobo/provi")

Nevertheless, this has no effect, and temp raster files keep growing
in the default directory within /tmp

Is this a bug?

Thanks
Agus

> sessionInfo()
R version 3.3.3 (2017-03-06)
Platform: i686-pc-linux-gnu (32-bit)
Running under: Debian GNU/Linux buster/sid

locale:
 [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
LC_TIME=en_US.UTF-8
 [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8   LC_NAME=C
LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] raster_2.5-8 sp_1.2-5

loaded via a namespace (and not attached):
[1] rgdal_1.2-11parallel_3.3.3  tools_3.3.3 Rcpp_0.12.12
grid_3.3.3
[6] lattice_0.20-35

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] temp directory not actually reset by rasterOptions()

[Agustin Lobo]

Thanks. Yes, I also had found using Simon Urbanek’s unixtools to be an
alternative solution.
But this does not address the core of my question:
raster::rasterOptions(tmpdir="/home/alobo/provi")
does not work.
Agus

On Thu, Sep 21, 2017 at 12:34 AM, Michael Sumner  wrote:
> FYI It's also important to set the max size for raster,  default is way too
> low:
>
> https://discuss.ropensci.org/t/how-to-avoid-space-hogging-raster-tempfiles/864
>
> On Thu, 21 Sep 2017, 01:25 Agustin Lobo  wrote:
>
>> I have problems with an small /tmp partition
>> and thought using
>> rasterOptions(tmpdir="/home/alobo/provi")
>>
>> Nevertheless, this has no effect, and temp raster files keep growing
>> in the default directory within /tmp
>>
>> Is this a bug?
>>
>> Thanks
>> Agus
>>
>> > sessionInfo()
>> R version 3.3.3 (2017-03-06)
>> Platform: i686-pc-linux-gnu (32-bit)
>> Running under: Debian GNU/Linux buster/sid
>>
>> locale:
>>  [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
>> LC_TIME=en_US.UTF-8
>>  [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
>> LC_MESSAGES=en_US.UTF-8
>>  [7] LC_PAPER=en_US.UTF-8   LC_NAME=C
>> LC_ADDRESS=C
>> [10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8
>> LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] stats graphics  grDevices utils datasets  methods   base
>>
>> other attached packages:
>> [1] raster_2.5-8 sp_1.2-5
>>
>> loaded via a namespace (and not attached):
>> [1] rgdal_1.2-11parallel_3.3.3  tools_3.3.3 Rcpp_0.12.12
>> grid_3.3.3
>> [6] lattice_0.20-35
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> Dr. Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> 203 Channel Highway
> Kingston Tasmania 7050 Australia
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Change in rasterToPolygons() behaviour

[Agustin Lobo]

Recycling an old script, I've found that rasterToPolygons()
produces an error if not all stack bands are of the same type:

file <- system.file("external/test.grd", package="raster") s <-
stack(file, file, file)
s2 <- (s[[3]]-s[[2]])/(s[[3]]-s[[2]])
s1 <- stack(file, file, file)
s2 <- (s1[[3]]-s1[[2]])/(s1[[3]]-s1[[2]])
s <- addLayer(s1,s2)
sg <- rasterToPolygons(s,fun=NULL)
Error in `colnames<-`(`*tmp*`, value = c("x", "y", names(x))) :
  length of 'dimnames' [2] not equal to array extent

This used to work. The following works now:
 s <- addLayer(s1*1.0,s2)
 sg <- rasterToPolygons(s,fun=NULL)

Just reporting for other inattentive users like me (if any).

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] intersections to individual polygons

[Agustin Lobo]

Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and
,B,  is there a function
that would split the original polygons into "onlyA", "onlyB" and
"AintersectingB" polygons?

Thanks
Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] intersections to individual polygons

[Agustin Lobo]

Mimicking your vignette:

require(rgdal)
require(raster)
require(spatstat)
require(rgeos)
require(maptools)

dzfile 
<-"https://www.dropbox.com/s/xxujcwqy3ec21sa/TDM1_DEM__04_S63W062_epsg32720.zip?dl=0";
download.file(dzfile,"TDM1_DEM__04_S63W062_epsg32720.zip",method="wget")
unzip("TDM1_DEM__04_S63W062_epsg32720.zip")
v <- readOGR(dsn="TDM1_DEM__04_S63W062_epsg32720",
 layer="TDM1_DEM__04_S63W062_epsg32720", stringsAsFactors = FALSE)
plot(v)
p <- slot(v, "polygons")
p2 <- lapply(p, function(x) { SpatialPolygons(list(x)) })
w <- lapply(p2, as.owin) #maptools required

te <- tess(tiles=w)
class(te)
plot.tess(te,do.labels=TRUE)

But this is still the original polygons,
not the mosaic of polygon parts I'm looking for.
Would I need to go doing all possible intersections and then
adding the "A not B" and "B not A" parts for all possible pairs?

Or is there a function in spatstat to convert "te" into a tesselation
of adjacent polygons?

Thanks



On Thu, Jun 28, 2018 at 2:02 PM, Rolf Turner  wrote:
>
> On 29/06/18 00:51, Agustin Lobo wrote:
>
>> Given an SpatialPolygons DF with i.e. 2 intersecting polygons, A and
>> ,B,  is there a function
>> that would split the original polygons into "onlyA", "onlyB" and
>> "AintersectingB" polygons?
>
>
> You can do this easily in spatstat by converting your polygons to "owin"
> objects.
>
> cheers,
>
> Rolf Turner
>
> --
> Technical Editor ANZJS
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] gdalUtils::gdal_translate()

[Agustin Lobo]

I find the following problem when trying gdal_translate():
>  gdal_translate(src_dataset=file.path(dirMAJAIma,"imain20.vrt"),
+ dst_dataset="test.tif",
+ ot="Int16",projwin=ext@xmin,ext@ymax,ext@xmax,ext@ymin)
Error in gdal_chooseInstallation(hasDrivers = of) :
  No installations match.

Other functions in gdalUtils do work and gdal_translate works on the
Debian terminal.
Does igdal_translate() work for other people?
Thanks
Agus

sessionInfo()
R version 3.5.0 (2018-04-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux buster/sid

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8   LC_NUMERIC=C
LC_TIME=en_GB.UTF-8
 [4] LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8
LC_MESSAGES=en_GB.UTF-8
 [7] LC_PAPER=en_GB.UTF-8   LC_NAME=C
LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_GB.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] gdalUtils_2.0.1.14 rgdal_1.3-6R.utils_2.7.0
R.oo_1.22.0R.methodsS3_1.7.1
[6] raster_2.7-15  sp_1.3-1

loaded via a namespace (and not attached):
 [1] compiler_3.5.0   tools_3.5.0  yaml_2.2.0   Rcpp_0.12.19
  codetools_0.2-15 grid_3.5.0
 [7] iterators_1.0.10 foreach_1.4.4knitr_1.20   lattice_0.20-35

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] gdalUtils::gdal_translate()

[Agustin Lobo]

Alex,
It was not a problem with the gdal installation, I had already checked:
> gdal_setInstallation()
> getOption("gdalUtils_gdalPath")[[1]]$path
[1] "/usr/bin/"
> getOption("gdalUtils_gdalPath")[[1]]$version
version
"2.3.2"
> getOption("gdalUtils_gdalPath")[[1]]$date
date
"2018-09-21"

Also note that the rest of gdalUtilities do work.

But thanks to your spotting of the missing c(), the error message
vanishes out. So it was a problem of my mistake with the missing c()
+ a wrong error message that was totally misleading me.

Many thanks!

Agus

On Wed, Nov 7, 2018 at 7:18 PM Alex M  wrote:
>
> gdal_translate does work fine for me on R 3.4
>
> That error sounds like your install is having issues picking what gdal
> installation to use.
> ?gdal_chooseInstallation
>
> perhaps your install was not found correctly amd you need to run
> gdal_setInstallation?
>
> But I also see something in your call that looks off to me.
>
> Based on the docs projwin takes a numeric vector. Which to me should
> look like projwin=c(ext@xmin,ext@ymax,ext@xmax,ext@ymin)
>
> Thanks,
> Alex
>
> On 11/7/18 01:02, Agustin Lobo wrote:
> > I find the following problem when trying gdal_translate():
> >>  gdal_translate(src_dataset=file.path(dirMAJAIma,"imain20.vrt"),
> > + dst_dataset="test.tif",
> > + ot="Int16",projwin=ext@xmin,ext@ymax,ext@xmax,ext@ymin)
> > Error in gdal_chooseInstallation(hasDrivers = of) :
> >   No installations match.
> >
> > Other functions in gdalUtils do work and gdal_translate works on the
> > Debian terminal.
> > Does igdal_translate() work for other people?
> > Thanks
> > Agus
> >
> > sessionInfo()
> > R version 3.5.0 (2018-04-23)
> > Platform: x86_64-pc-linux-gnu (64-bit)
> > Running under: Debian GNU/Linux buster/sid
> >
> > Matrix products: default
> > BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
> > LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
> >
> > locale:
> >  [1] LC_CTYPE=en_GB.UTF-8   LC_NUMERIC=C
> > LC_TIME=en_GB.UTF-8
> >  [4] LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8
> > LC_MESSAGES=en_GB.UTF-8
> >  [7] LC_PAPER=en_GB.UTF-8   LC_NAME=C
> > LC_ADDRESS=C
> > [10] LC_TELEPHONE=C LC_MEASUREMENT=en_GB.UTF-8
> > LC_IDENTIFICATION=C
> >
> > attached base packages:
> > [1] stats graphics  grDevices utils datasets  methods   base
> >
> > other attached packages:
> > [1] gdalUtils_2.0.1.14 rgdal_1.3-6R.utils_2.7.0
> > R.oo_1.22.0R.methodsS3_1.7.1
> > [6] raster_2.7-15  sp_1.3-1
> >
> > loaded via a namespace (and not attached):
> >  [1] compiler_3.5.0   tools_3.5.0  yaml_2.2.0   Rcpp_0.12.19
> >   codetools_0.2-15 grid_3.5.0
> >  [7] iterators_1.0.10 foreach_1.4.4knitr_1.20   lattice_0.20-35
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] writeRaster: BSQ despite BIL is requested

[Agustin Lobo]

I do

writeRaster(testin,
file="testin",format="ENVI",datatype="INT2U",bandorder="BIL",overwrite=TRUE)

but get an image with BSQ order
(rda object: https://www.dropbox.com/s/or5h7yb9nd52e94/testin.rda?dl=0)
Thanks
Agus

R version 3.5.0 (2018-04-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux buster/sid

Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0

locale:
 [1] LC_CTYPE=en_GB.UTF-8   LC_NUMERIC=C
LC_TIME=en_GB.UTF-8LC_COLLATE=en_GB.UTF-8
LC_MONETARY=en_GB.UTF-8
 [6] LC_MESSAGES=en_GB.UTF-8LC_PAPER=en_GB.UTF-8   LC_NAME=C
   LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] raster_2.7-15 rgdal_1.3-6   sp_1.3-1


Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] writeRaster: BSQ despite BIL is requested

[Agustin Lobo]

Thanks Roger,

Just to summarize, not only the hdr files are correct, the hdr files
are consistent with the ordering in the data files:

require(rgdal)
require(raster)
m <- matrix(rep(10,20),ncol=5)
b <- brick(raster(m),raster(10*m),raster(100*m))
b
writeRaster(b, file="bwrast", format="ENVI", datatype="INT2U",
bandorder="BIL", overwrite=TRUE)

writes a correct BSQ file, despite having requested BIL order

writeRaster(b, file="bwrast2", format="ENVI", datatype="INT2U",
options="INTERLEAVE=BIL", overwrite=TRUE)

writes a correct BIL file (and this is a useful solution, thanks).

I tend to think that the fact that writeRaster() ignores
bandorder="BIL" is a bug.

Agus
(mi R is using Loaded GDAL runtime: GDAL 2.3.2, released 2018/09/21)

On Wed, Dec 19, 2018 at 11:52 AM Roger Bivand  wrote:
>
> On Wed, 19 Dec 2018, Agustin Lobo wrote:
>
> > I do
> >
> > writeRaster(testin,
> > file="testin",format="ENVI",datatype="INT2U",bandorder="BIL",overwrite=TRUE)
> >
> > but get an image with BSQ order
> > (rda object: https://www.dropbox.com/s/or5h7yb9nd52e94/testin.rda?dl=0)
>
> Which version of GDAL is under all of this? It doesn't matter here, but
> might have mattered.
>
> raster::writeRaster() puts interleave = bsq in the *.hdr files. "BIP" and
> does not seem to do anything. However:
>
> writeGDAL(as(testin, "SpatialGridDataFrame"), fname="testinsp.envi",
> drivername="ENVI", type="UInt16", options="INTERLEAVE=BIL",
> setStatistics=TRUE)
>
> yields interleave = bil in the *.hdr file. Using the same route, I get:
>
> writeRaster(testin, file="testin", format="ENVI", datatype="INT2U",
> options="INTERLEAVE=BIL", overwrite=TRUE)
>
> with interleave = bil in the *.hdr file. I'm not sure where
> raster::writeRaster() drops the options definitions given in bandorder=,
> maybe ENVI is not a 'native' file format?
>
> Hope this helps,
>
> Roger
>
>
> > Thanks
> > Agus
> >
> > R version 3.5.0 (2018-04-23)
> > Platform: x86_64-pc-linux-gnu (64-bit)
> > Running under: Debian GNU/Linux buster/sid
> >
> > Matrix products: default
> > BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
> > LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
> >
> > locale:
> > [1] LC_CTYPE=en_GB.UTF-8   LC_NUMERIC=C
> > LC_TIME=en_GB.UTF-8LC_COLLATE=en_GB.UTF-8
> > LC_MONETARY=en_GB.UTF-8
> > [6] LC_MESSAGES=en_GB.UTF-8LC_PAPER=en_GB.UTF-8   LC_NAME=C
> >   LC_ADDRESS=C   LC_TELEPHONE=C
> > [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
> >
> > attached base packages:
> > [1] stats graphics  grDevices utils datasets  methods   base
> >
> > other attached packages:
> > [1] raster_2.7-15 rgdal_1.3-6   sp_1.3-1
> >
> >
> > Agus
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J&hl=en
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Brick and Stack in package raster

[Agustin Lobo]

Thanks
Just noting that you also have BIL for envi format and that
BIP seems to be the default for geotiff (at least that's what gdalinfo
reports for
tif files written by writeRaster() :
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL

Also, I understand that the order is BIL also when brick() writes a
file by default (i.e. /tmp/R_raster_tmp/raster_tmp_789131022710.grd)

Agus

2010/11/26 Robert J. Hijmans :
> On Fri, Nov 26, 2010 at 2:25 AM, Agustin Lobo  wrote:
>> Robert,
>>
>> At some point you made this comment:
>>
>>>>> (as you are dealing with large
>>>>> files, I think it is more efficient to first make a stack and then
>>>>> write the results to a brick)
>>
>> But I have some doubts:
>>
>> bN = brick(N1,N2,N3,...,N72)
>> where
>> N1 = raster(paste(midir,"NDV_19990101__Extract.img",sep="")) etc
>> takes a very long time. Is a file on disk being written somewhere?
>>
>
> That is why I am saying that it is better to do
>
> sN = stack(N1,N2,N3,...,N72)
>
> and then, perhaps,
>
> bN <- brick(sN )
>
> But that is probably not worth it. I may change a few things in the
> code and manual to make this go smoother.
> I think the approach is to always makes a stack if you have separate
> files and a brick if it is a single file. Adding a layer to a brick
> should perhaps result in a stack. Computations on a stack/brick always
> result in a brick (single file if on disk)
>
>
>> Is overlay() going to be any faster with brick than with stack objects? ie:
>> bNv2 <- overlay(bN, bSM,
>> fun=fun1,filename="bNv2.tif",overwrite=TRUE,NAflag=0) faster than
>> sNv2 <- overlay(sN, sSM, 
>> fun=fun1,filename="sNv2.tif",overwrite=TRUE,NAflag=0) ?
>> (bN, bSM, brick objects; sN, sSM, stack objects)
>
> Yes, I would think that would be faster with a brick, because the data
> need to come out of a single object/file (i.e. 72 times fewer files
> need to be opened and closed). But I am not sure if it is worth it to
> first make a brick from the stack, because that also takes time.
>
>> When you write a brick to a file, can you control if the resulting
>> file is going to be BSQ or BIL?
>
> Yes, but only for the native format (not when writing via rgdal, e.g.
> GeoTIFF). BIL is default and should work. BIP is rather inefficient
> except for extracting point values. BSQ best for extracting single
> layers.
>
>> Are the files written by writeRaster() identical for stack and brick objects?
>> writeRaster(bNv2,
>> filename="bNv2b.tif",overwrite=TRUE,format="GTiff",datatype='INT2U',NAflag=0)
>> writeRaster(sNv2,
>> filename="sNv2b.tif",overwrite=TRUE,format="GTiff",datatype='INT2U',NAflag=0)
>>
>
> Yes, always writes a single file (brick object).
>
>
>> Thanks
>>
>> Agus
>>
>
> Cheers, Robert
>
> ___
> R-sig-Geo mailing list
> r-sig-...@stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster::crop

[Agustin Lobo]

Is it possible to crop a raster by selecting rows and cols? i.e. the
equivalent to

subs = bN3[1000:1100,1000:1100]

but subs should remain a raster (or a brick)
(crop seems to require geographic coordinates even if the raster is projected)

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] warning at creating SpatialPolygonsDataFrame and error at saving through rgdal

[Agustin Lobo]

Hi!

I create a SPDF with a function
HypNAr = meta2sppol()

that ends up with:

sitios <- 1:length(noms)
row.names(datos) <- sitios
SpatialPolygonsDataFrame(dummy, data=datos, match.ID = T)

but get:

Warning message:
In data.row.names(row.names, rowsi, i) :
  some row.names duplicated: 2 --> row.names NOT used

I do not understand, as I define the row names as a sequence of
integers, how can they be duplicated?
I can plot the SPDF and the data slot is correct, but when I try to
save as shapefile:

> writeOGR(HypNAr,dsn="HypNAr",layer="HypNAr",driver="ESRI Shapefile",verbose=T)
Error in writeOGR(HypNAr, dsn = "HypNAr", layer = "HypNAr", driver =
"ESRI Shapefile",  :

GDAL Error 1: Invalid index : -1
Calls: writeOGR -> .Call
In addition: There were 11 warnings (use warnings() to see them)

Is this error related to the previous warning?

Thanks!

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster::crop

[Agustin Lobo]

Thanks,

and summarizing messages exchanged with Robert,

subs = bN3[1000:1100,1000:1100, drop=FALSE]

is the easiest, while

e<- extent(c (xFromCol(bN3, 1000), xFromCol(bN3, 1100), yFromRow(bN3,
1100), yFromRow(bN3, 1000)))
x<- crop(bN3, e)

is the fastest. And Robert plans a simplifies syntax for using extent
with row, cols in the future, i.e. like:
extent(x, 1000,1100,1000,1100)

Agus

2010/11/29 steven mosher :
>  you should be able to use xFromCol, yFromCol  etc  to get the coordinate
> for
> the extent you want to crop. or xyFromCell.
>
>
>
>
>
> On Sun, Nov 28, 2010 at 9:47 AM, Agustin Lobo  wrote:
>
>> Is it possible to crop a raster by selecting rows and cols? i.e. the
>> equivalent to
>>
>> subs = bN3[1000:1100,1000:1100]
>>
>> but subs should remain a raster (or a brick)
>> (crop seems to require geographic coordinates even if the raster is
>> projected)
>>
>> Thanks
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>        [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] warning at creating SpatialPolygonsDataFrame and error at saving through rgdal

[Agustin Lobo]

The error in wirteOGR() was caused by the long names of the variables
in the table
of the SPDF object that were listed in the warnings. Fixed by renaming
some columns in the data slot to shorter names.

The error was not related to the warnings issued at creating the SPDF
object by my function,
those were easily fixed and are of no general interest (and of particular shame)

I wonder if the error  message in writeOGR() could be more
informative. Actually, I initially discarded the long names
as the cause of the error because, according to the warning messages,
those names had been automatically truncated.

Thanks!

Agus

2010/11/29 Agustin Lobo :
> Roger,
>
> Please find enclosed the function and 3 files such as the ones that
> are processed. Please extract the zip
> into a folder named GLIS_scenes and indicate the path in the call to the
> function.
>
> The function is run just by:
>
> source("meta2sppol.R")
> HypAr <- meta2sppol(midir="../GLIS_scenes")
>
> What I try to do is read each file, extract the 4 pairs of coordinates
> and make a SPDF with the polygons and the data
> from the files.
>
> Thanks
>
> Agus
>
> 2010/11/29 Roger Bivand :
>> On Mon, 29 Nov 2010, Agustin Lobo wrote:
>>
>>> Hi!
>>>
>>> I create a SPDF with a function
>>> HypNAr = meta2sppol()
>>>
>>> that ends up with:
>>>
>>> sitios <- 1:length(noms)
>>> row.names(datos) <- sitios
>>> SpatialPolygonsDataFrame(dummy, data=datos, match.ID = T)
>>>
>>> but get:
>>>
>>> Warning message:
>>> In data.row.names(row.names, rowsi, i) :
>>>  some row.names duplicated: 2 --> row.names NOT used
>>
>> The message is coming from data.frame() in base package. Do
>>
>> debug(meta2sppol)
>> debug(data.frame)
>> HypNAr = meta2sppol()
>>
>> and see what happens inside the data.frame() call that creates datos, which
>> is where the warning is issued (it may be inside a coercion method call,
>> such as as.data.frame()); open:
>>
>> https://svn.r-project.org/R/trunk/src/library/base/R/dataframe.R
>>
>> and search for "some row.names duplicated".
>>
>> That is, the problem is in your own function. Please also note that
>> providing a reproducible example is crucial to getting precise assistance,
>> and usually the exercise of writing such an example reveals the cause, thus
>> solving the problem more quickly.
>>
>> Roger
>>
>>>
>>> I do not understand, as I define the row names as a sequence of
>>> integers, how can they be duplicated?
>>> I can plot the SPDF and the data slot is correct, but when I try to
>>> save as shapefile:
>>>
>>>> writeOGR(HypNAr,dsn="HypNAr",layer="HypNAr",driver="ESRI
>>>> Shapefile",verbose=T)
>>>
>>> Error in writeOGR(HypNAr, dsn = "HypNAr", layer = "HypNAr", driver =
>>> "ESRI Shapefile",  :
>>>
>>>        GDAL Error 1: Invalid index : -1
>>> Calls: writeOGR -> .Call
>>> In addition: There were 11 warnings (use warnings() to see them)
>>>
>>> Is this error related to the previous warning?
>>>
>>> Thanks!
>>>
>>> Agus
>>>
>>> ___
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>> --
>> Roger Bivand
>> Economic Geography Section, Department of Economics, Norwegian School of
>> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
>> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
>> e-mail: roger.biv...@nhh.no
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

Hi!

Could an option overwrite=T (default overwrite=F) be included in writeOGR()?
Thanks!

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

Roger,

It does not for me. The error is solved by deleting files xygcatLO.* ,
while directory xygcatLO may continue to exist.

Also, I've found that with this version the below command does not
save the files within
the folder indicated by dsn, but outside. This was not the case in the
past, but do not know if this is
because of a change in rgdal, gdal or ubuntu 10.04 (I've changed all 3
since the last time I used writeOGR()).

> writeOGR(xygcatLO,dsn="xygcatLO",layer="xygcatLO", driver="ESRI Shapefile")
Error in writeOGR(xygcatLO, dsn = "xygcatLO", layer = "xygcatLO",
driver = "ESRI Shapefile") :

GDAL Error 1: Layer 'xygcatLO' already exists
Calls: writeOGR -> .Call

> sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.utf8  LC_NUMERIC=C
 [3] LC_TIME=en_US.utf8   LC_COLLATE=en_US.utf8
 [5] LC_MONETARY=en_US.utf8   LC_MESSAGES=en_US.utf8
 [7] LC_PAPER=en_US.utf8  LC_NAME=en_US.utf8
 [9] LC_ADDRESS=en_US.utf8LC_TELEPHONE=en_US.utf8
[11] LC_MEASUREMENT=en_US.utf8LC_IDENTIFICATION=en_US.utf8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] RANN_2.1.2  outliers_0.13-3 rgdal_0.6-30raster_1.7-2
[5] sp_0.9-73   rkward_0.5.4

loaded via a namespace (and not attached):
[1] grid_2.12.0 lattice_0.19-13 tools_2.12.0
2010/12/10 Roger Bivand :
> On Thu, 9 Dec 2010, Barry Rowlingson wrote:
>
>> On Thu, Dec 9, 2010 at 8:56 AM, Agustin Lobo 
>> wrote:
>>>
>>> Hi!
>>>
>>> Could an option overwrite=T (default overwrite=F) be included in
>>> writeOGR()?
>>> Thanks!
>>
>> Submit a Feature Request here:
>>
>> http://sourceforge.net/tracker/?group_id=84716&atid=573625
>
> No, rgdal is on R-Forge, and a list called rgdal-devel will be available
> from tomorrow.
>
>>
>> It looks like it just needs to be a wrapper around OGR's delete layer
>> function.
>
> Contributions welcome. Note that the facility is present for rasters, which
> I fould out by accidentally deleting GRASS database objects (for which of
> course I did not have backups). Is this really needed - which drivers do not
> already overwrite without complaint? These do not complain for me:
> driver="ESRI Shapefile"; driver="KML"; driver="MapInfo File";
> driver="MapInfo File", dataset_options="FORMAT=MIF".
>
> Hope this helps,
>
> Roger
>
>>
>> Barry
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: roger.biv...@nhh.no
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

Roger,

is this all you need? I thought sessionInfo was providing everything
needed. The driver is
"ESRI Shapefile" as stated in the command, have not tested any other
(it does work for
raster though, as you mention):

> require(rgdal)
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.8dev, released 2010/01/19
Path to GDAL shared files: /usr/local/share/gdal
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)

Thanks

Agus

2010/12/10 Roger Bivand :
> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>
>> Roger,
>>
>> It does not for me. The error is solved by deleting files xygcatLO.* ,
>> while directory xygcatLO may continue to exist.
>
> Agus,
>
> With which drivers, OGR version (in header when rgdal loads)? If I cannot
> reproduce this, I cannot readily help.
>
> Roger
>
>>
>> Also, I've found that with this version the below command does not
>> save the files within
>> the folder indicated by dsn, but outside. This was not the case in the
>> past, but do not know if this is
>> because of a change in rgdal, gdal or ubuntu 10.04 (I've changed all 3
>> since the last time I used writeOGR()).
>>
>>> writeOGR(xygcatLO,dsn="xygcatLO",layer="xygcatLO", driver="ESRI
>>> Shapefile")
>>
>> Error in writeOGR(xygcatLO, dsn = "xygcatLO", layer = "xygcatLO",
>> driver = "ESRI Shapefile") :
>>
>>        GDAL Error 1: Layer 'xygcatLO' already exists
>> Calls: writeOGR -> .Call
>>
>>> sessionInfo()
>>
>> R version 2.12.0 (2010-10-15)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>>
>> locale:
>> [1] LC_CTYPE=en_US.utf8          LC_NUMERIC=C
>> [3] LC_TIME=en_US.utf8           LC_COLLATE=en_US.utf8
>> [5] LC_MONETARY=en_US.utf8       LC_MESSAGES=en_US.utf8
>> [7] LC_PAPER=en_US.utf8          LC_NAME=en_US.utf8
>> [9] LC_ADDRESS=en_US.utf8        LC_TELEPHONE=en_US.utf8
>> [11] LC_MEASUREMENT=en_US.utf8    LC_IDENTIFICATION=en_US.utf8
>>
>> attached base packages:
>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>
>> other attached packages:
>> [1] RANN_2.1.2      outliers_0.13-3 rgdal_0.6-30    raster_1.7-2
>> [5] sp_0.9-73       rkward_0.5.4
>>
>> loaded via a namespace (and not attached):
>> [1] grid_2.12.0     lattice_0.19-13 tools_2.12.0
>> 2010/12/10 Roger Bivand :
>>>
>>> On Thu, 9 Dec 2010, Barry Rowlingson wrote:
>>>
>>>> On Thu, Dec 9, 2010 at 8:56 AM, Agustin Lobo 
>>>> wrote:
>>>>>
>>>>> Hi!
>>>>>
>>>>> Could an option overwrite=T (default overwrite=F) be included in
>>>>> writeOGR()?
>>>>> Thanks!
>>>>
>>>> Submit a Feature Request here:
>>>>
>>>> http://sourceforge.net/tracker/?group_id=84716&atid=573625
>>>
>>> No, rgdal is on R-Forge, and a list called rgdal-devel will be available
>>> from tomorrow.
>>>
>>>>
>>>> It looks like it just needs to be a wrapper around OGR's delete layer
>>>> function.
>>>
>>> Contributions welcome. Note that the facility is present for rasters,
>>> which
>>> I fould out by accidentally deleting GRASS database objects (for which of
>>> course I did not have backups). Is this really needed - which drivers do
>>> not
>>> already overwrite without complaint? These do not complain for me:
>>> driver="ESRI Shapefile"; driver="KML"; driver="MapInfo File";
>>> driver="MapInfo File", dataset_options="FORMAT=MIF".
>>>
>>> Hope this helps,
>>>
>>> Roger
>>>
>>>>
>>>> Barry
>>>>
>>>> ___
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo@r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>> --
>>> Roger Bivand
>>> Economic Geography Section, Department of Economics, Norwegian School of
>>> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
>>> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
>>> e-mail: roger.biv...@nhh.no
>>>
>>> ___
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>
>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: roger.biv...@nhh.no
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

Roger,

it works fine on another computer:

> require(rgdal)
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.2, released 2010/04/23
Path to GDAL shared files: /usr/share/gdal17
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)

> writeOGR(xygcatLO,dsn="xygcatLO",layer="xygcatLO", driver="ESRI Shapefile")
> sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: i486-pc-linux-gnu (32-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C
LC_TIME=en_US.UTF-8
 [4] LC_COLLATE=en_US.UTF-8LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8  LC_NAME=en_US.UTF-8
LC_ADDRESS=en_US.UTF-8
[10] LC_TELEPHONE=en_US.UTF-8  LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.6-27  raster_1.6-22 sp_0.9-72 rkward_0.5.4

loaded via a namespace (and not attached):
[1] grid_2.12.0 lattice_0.19-13 tools_2.12.0

Agus


2010/12/10 Agustin Lobo :
> Roger,
>
> is this all you need? I thought sessionInfo was providing everything
> needed. The driver is
> "ESRI Shapefile" as stated in the command, have not tested any other
> (it does work for
> raster though, as you mention):
>
>> require(rgdal)
> Loading required package: rgdal
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.8dev, released 2010/01/19
> Path to GDAL shared files: /usr/local/share/gdal
> Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
> Path to PROJ.4 shared files: (autodetected)
>
> Thanks
>
> Agus
>
> 2010/12/10 Roger Bivand :
>> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>>
>>> Roger,
>>>
>>> It does not for me. The error is solved by deleting files xygcatLO.* ,
>>> while directory xygcatLO may continue to exist.
>>
>> Agus,
>>
>> With which drivers, OGR version (in header when rgdal loads)? If I cannot
>> reproduce this, I cannot readily help.
>>
>> Roger
>>
>>>
>>> Also, I've found that with this version the below command does not
>>> save the files within
>>> the folder indicated by dsn, but outside. This was not the case in the
>>> past, but do not know if this is
>>> because of a change in rgdal, gdal or ubuntu 10.04 (I've changed all 3
>>> since the last time I used writeOGR()).
>>>
>>>> writeOGR(xygcatLO,dsn="xygcatLO",layer="xygcatLO", driver="ESRI
>>>> Shapefile")
>>>
>>> Error in writeOGR(xygcatLO, dsn = "xygcatLO", layer = "xygcatLO",
>>> driver = "ESRI Shapefile") :
>>>
>>>        GDAL Error 1: Layer 'xygcatLO' already exists
>>> Calls: writeOGR -> .Call
>>>
>>>> sessionInfo()
>>>
>>> R version 2.12.0 (2010-10-15)
>>> Platform: x86_64-pc-linux-gnu (64-bit)
>>>
>>> locale:
>>> [1] LC_CTYPE=en_US.utf8          LC_NUMERIC=C
>>> [3] LC_TIME=en_US.utf8           LC_COLLATE=en_US.utf8
>>> [5] LC_MONETARY=en_US.utf8       LC_MESSAGES=en_US.utf8
>>> [7] LC_PAPER=en_US.utf8          LC_NAME=en_US.utf8
>>> [9] LC_ADDRESS=en_US.utf8        LC_TELEPHONE=en_US.utf8
>>> [11] LC_MEASUREMENT=en_US.utf8    LC_IDENTIFICATION=en_US.utf8
>>>
>>> attached base packages:
>>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>>
>>> other attached packages:
>>> [1] RANN_2.1.2      outliers_0.13-3 rgdal_0.6-30    raster_1.7-2
>>> [5] sp_0.9-73       rkward_0.5.4
>>>
>>> loaded via a namespace (and not attached):
>>> [1] grid_2.12.0     lattice_0.19-13 tools_2.12.0
>>> 2010/12/10 Roger Bivand :
>>>>
>>>> On Thu, 9 Dec 2010, Barry Rowlingson wrote:
>>>>
>>>>> On Thu, Dec 9, 2010 at 8:56 AM, Agustin Lobo 
>>>>> wrote:
>>>>>>
>>>>>> Hi!
>>>>>>
>>>>>> Could an option overwrite=T (default overwrite=F) be included in
>>>>>> writeOGR()?
>>>>>> Thanks!
>>>>>
>>>>> Submit a Feature Request here:
>>>>>
>>>>> http://sourceforge.net/tracker/?group_id=84716&atid=573625
>>>>
>>>> No, rgdal is on R-Forge, and a list called rgdal-devel will be available
>>>> from tomorrow.
&

Re: [R-sig-Geo] Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

Roger,

I was not aware that 1.8dev is actually older than 1.7.2
I'll check why I have that gdal on the other machine and will try
to fix that on my side on Monday. But perhaps this problem will come up
again for future stable versions of 1.8.* ? Anyway, 1.8dev is clearly not a good
choice and I have to change it.

Thanks

Agus

2010/12/10 Roger Bivand :
> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>
>> Roger,
>>
>> is this all you need? I thought sessionInfo was providing everything
>> needed.
>
> You had given sessionInfo() but that does not report which GDAL you are
> using. The report unfortunately does not give an svn revision, maybe
> gdal-config --version does? I develop using 1.7.* GDAL. Your 1.8dev is
> pretty old, certainly older that released 1.7.3. I can try to build a
> current GDAL trunk from source next week, if I get time, but this looks like
> something in your setup (development GDAL that is now stale?).
>
>> The driver is "ESRI Shapefile" as stated in the command, have not tested
>> any other (it does work for raster though, as you mention):
>
> I mentioned four variants of vector drivers, not raster.
>
> Roger
>
>>
>>> require(rgdal)
>>
>> Loading required package: rgdal
>> Geospatial Data Abstraction Library extensions to R successfully loaded
>> Loaded GDAL runtime: GDAL 1.8dev, released 2010/01/19
>> Path to GDAL shared files: /usr/local/share/gdal
>> Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
>> Path to PROJ.4 shared files: (autodetected)
>>
>> Thanks
>>
>> Agus
>>
>> 2010/12/10 Roger Bivand :
>>>
>>> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>>>
>>>> Roger,
>>>>
>>>> It does not for me. The error is solved by deleting files xygcatLO.* ,
>>>> while directory xygcatLO may continue to exist.
>>>
>>> Agus,
>>>
>>> With which drivers, OGR version (in header when rgdal loads)? If I cannot
>>> reproduce this, I cannot readily help.
>>>
>>> Roger
>>>
>>>>
>>>> Also, I've found that with this version the below command does not
>>>> save the files within
>>>> the folder indicated by dsn, but outside. This was not the case in the
>>>> past, but do not know if this is
>>>> because of a change in rgdal, gdal or ubuntu 10.04 (I've changed all 3
>>>> since the last time I used writeOGR()).
>>>>
>>>>> writeOGR(xygcatLO,dsn="xygcatLO",layer="xygcatLO", driver="ESRI
>>>>> Shapefile")
>>>>
>>>> Error in writeOGR(xygcatLO, dsn = "xygcatLO", layer = "xygcatLO",
>>>> driver = "ESRI Shapefile") :
>>>>
>>>>        GDAL Error 1: Layer 'xygcatLO' already exists
>>>> Calls: writeOGR -> .Call
>>>>
>>>>> sessionInfo()
>>>>
>>>> R version 2.12.0 (2010-10-15)
>>>> Platform: x86_64-pc-linux-gnu (64-bit)
>>>>
>>>> locale:
>>>> [1] LC_CTYPE=en_US.utf8          LC_NUMERIC=C
>>>> [3] LC_TIME=en_US.utf8           LC_COLLATE=en_US.utf8
>>>> [5] LC_MONETARY=en_US.utf8       LC_MESSAGES=en_US.utf8
>>>> [7] LC_PAPER=en_US.utf8          LC_NAME=en_US.utf8
>>>> [9] LC_ADDRESS=en_US.utf8        LC_TELEPHONE=en_US.utf8
>>>> [11] LC_MEASUREMENT=en_US.utf8    LC_IDENTIFICATION=en_US.utf8
>>>>
>>>> attached base packages:
>>>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>>>
>>>> other attached packages:
>>>> [1] RANN_2.1.2      outliers_0.13-3 rgdal_0.6-30    raster_1.7-2
>>>> [5] sp_0.9-73       rkward_0.5.4
>>>>
>>>> loaded via a namespace (and not attached):
>>>> [1] grid_2.12.0     lattice_0.19-13 tools_2.12.0
>>>> 2010/12/10 Roger Bivand :
>>>>>
>>>>> On Thu, 9 Dec 2010, Barry Rowlingson wrote:
>>>>>
>>>>>> On Thu, Dec 9, 2010 at 8:56 AM, Agustin Lobo 
>>>>>> wrote:
>>>>>>>
>>>>>>> Hi!
>>>>>>>
>>>>>>> Could an option overwrite=T (default overwrite=F) be included in
>>>>>>> writeOGR()?
>>>>>>> Thanks!
>>>>>>
>>>>>> Submit a Feature Request here:
>>>>>>
>>>>>> http://sourceforge.net/tracker/?group_id=84716&atid=573625
>&g

[R-sig-Geo] Fwd: Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

-- Forwarded message --
From: Agustin Lobo 
Date: 2010/12/13
Subject: Re: [R-sig-Geo] Suggestion: option overwrite=T for writeOGR()
To: roger.biv...@nhh.no


Roger,

According to Synaptic, my gdal is 1.7.3,
I do not understand where the reported
GDAL runtime: GDAL 1.8dev
comes from. Could it have been installed by rgdal?

I've searched my disk and find gdal16 and gdal17 but not gdal18
Do you know if it is possible to tell rgdal which gdal should be used?

Thanks

Agus



2010/12/11 Roger Bivand :
> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>
>> Roger,
>>
>> I was not aware that 1.8dev is actually older than 1.7.2 I'll check why I
>> have that gdal on the other machine and will try to fix that on my side on
>> Monday. But perhaps this problem will come up again for future stable
>> versions of 1.8.* ? Anyway, 1.8dev is clearly not a good choice and I have
>> to change it.
>
> As of GDAL 1.8.0 SVN revision 21234, the problem exists. I'll fix this in
> rgdal for the drivers I am familiar with before 1.8.0 is released.
>
> Roger
>
>>
>> Thanks
>>
>> Agus
>>
>> 2010/12/10 Roger Bivand :
>>>
>>> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>>>
>>>> Roger,
>>>>
>>>> is this all you need? I thought sessionInfo was providing everything
>>>> needed.
>>>
>>> You had given sessionInfo() but that does not report which GDAL you are
>>> using. The report unfortunately does not give an svn revision, maybe
>>> gdal-config --version does? I develop using 1.7.* GDAL. Your 1.8dev is
>>> pretty old, certainly older that released 1.7.3. I can try to build a
>>> current GDAL trunk from source next week, if I get time, but this looks
>>> like
>>> something in your setup (development GDAL that is now stale?).
>>>
>>>> The driver is "ESRI Shapefile" as stated in the command, have not tested
>>>> any other (it does work for raster though, as you mention):
>>>
>>> I mentioned four variants of vector drivers, not raster.
>>>
>>> Roger
>>>
>>>>
>>>>> require(rgdal)
>>>>
>>>> Loading required package: rgdal
>>>> Geospatial Data Abstraction Library extensions to R successfully loaded
>>>> Loaded GDAL runtime: GDAL 1.8dev, released 2010/01/19
>>>> Path to GDAL shared files: /usr/local/share/gdal
>>>> Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
>>>> Path to PROJ.4 shared files: (autodetected)
>>>>
>>>> Thanks
>>>>
>>>> Agus
>>>>
>>>> 2010/12/10 Roger Bivand :
>>>>>
>>>>> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>>>>>
>>>>>> Roger,
>>>>>>
>>>>>> It does not for me. The error is solved by deleting files xygcatLO.* ,
>>>>>> while directory xygcatLO may continue to exist.
>>>>>
>>>>> Agus,
>>>>>
>>>>> With which drivers, OGR version (in header when rgdal loads)? If I
>>>>> cannot
>>>>> reproduce this, I cannot readily help.
>>>>>
>>>>> Roger
>>>>>
>>>>>>
>>>>>> Also, I've found that with this version the below command does not
>>>>>> save the files within
>>>>>> the folder indicated by dsn, but outside. This was not the case in the
>>>>>> past, but do not know if this is
>>>>>> because of a change in rgdal, gdal or ubuntu 10.04 (I've changed all 3
>>>>>> since the last time I used writeOGR()).
>>>>>>
>>>>>>> writeOGR(xygcatLO,dsn="xygcatLO",layer="xygcatLO", driver="ESRI
>>>>>>> Shapefile")
>>>>>>
>>>>>> Error in writeOGR(xygcatLO, dsn = "xygcatLO", layer = "xygcatLO",
>>>>>> driver = "ESRI Shapefile") :
>>>>>>
>>>>>>        GDAL Error 1: Layer 'xygcatLO' already exists
>>>>>> Calls: writeOGR -> .Call
>>>>>>
>>>>>>> sessionInfo()
>>>>>>
>>>>>> R version 2.12.0 (2010-10-15)
>>>>>> Platform: x86_64-pc-linux-gnu (64-bit)
>>>>>>
>>>>>> locale:
>>>>>> [1] LC_CTYPE=en_US.utf8          LC_NUMERIC=C
>>>>>> [3] LC_TIME=en_U

Re: [R-sig-Geo] Fwd: Suggestion: option overwrite=T for writeOGR()

[Agustin Lobo]

Roger,

Your new rgdal version from rgdal_0.6-32.tar.gz works fine once I've
removed gdal1.8dev (using gdal1.7.3 now):
overwrite_layer=T works and the files go to the directory indicated by dsn

Thanks
Agus

> setwd ("/media/FREECOM_HDD/GRAVI")
> require(rgdal)
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.3, released 2010/11/10
Path to GDAL shared files: /usr/share/gdal17
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)

> writeOGR(xygcatLO,dsn="/media/FREECOM_HDD/GRAVI/xygcatLO",layer="xygcatLO", 
> driver="ESRI Shapefile",overwrite_layer=T)
> writeOGR(xygcatLO,dsn="/media/FREECOM_HDD/GRAVI/xygcatLO",layer="xygcatLO", 
> driver="ESRI Shapefile",overwrite_layer=T)
> writeOGR(xygcatLO,dsn="/media/FREECOM_HDD/GRAVI/xygcatLO",layer="xygcatLO", 
> driver="ESRI Shapefile",overwrite_layer=T)
> sessionInfo()
R version 2.12.0 (2010-10-15)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.utf8  LC_NUMERIC=C
 [3] LC_TIME=en_US.utf8   LC_COLLATE=en_US.utf8
 [5] LC_MONETARY=en_US.utf8   LC_MESSAGES=en_US.utf8
 [7] LC_PAPER=en_US.utf8  LC_NAME=en_US.utf8
 [9] LC_ADDRESS=en_US.utf8LC_TELEPHONE=en_US.utf8
[11] LC_MEASUREMENT=en_US.utf8LC_IDENTIFICATION=en_US.utf8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.6-32 raster_1.7-2 sp_0.9-74rkward_0.5.4

loaded via a namespace (and not attached):
[1] grid_2.12.0 lattice_0.19-13 tools_2.12.0

2010/12/14 Roger Bivand :
> On Mon, 13 Dec 2010, Roger Bivand wrote:
>
>> On Mon, 13 Dec 2010, Agustin Lobo wrote:
>>
>>> -- Forwarded message --
>>> From: Agustin Lobo 
>>> Date: 2010/12/13
>>> Subject: Re: [R-sig-Geo] Suggestion: option overwrite=T for writeOGR()
>>> To: roger.biv...@nhh.no
>>>
>>>
>>> Roger,
>>>
>>> According to Synaptic, my gdal is 1.7.3,
>>> I do not understand where the reported
>>> GDAL runtime: GDAL 1.8dev
>>> comes from. Could it have been installed by rgdal?
>>
>> rgdal package binaries from CRAN for Windows are built static against
>> GDAL, so include GDAL but not as a shared object. For OSX on CRAN extras,
>> I'm not sure whether the binary rgdal package is built static or not, but
>> GDAL is included.
>>
>> You are on Linux, so you provide GDAL yourself. If it declares on query
>> (asking the shared object what version it is, which is what rgdal does on
>> startup) that it is 1.8.0dev, then that is what it is. Be aware that
>> downstream packagers may also mistake the version numbers, so asking the
>> shared object is authoritative (unless the releasers of the original source
>> enter the wrong string, as with the current PROJ.4, which declares that it
>> is an as yet unreleased version, these things happen).
>>
>> Either revert properly to GDAL 1.7.*, or wait for me to fix this for
>> 1.8.*, please. I'll post in this thread when I'm ready, and then you can
>> check the rgdal source from R-Forge, OK?
>
> Could you, or others who may be interested, please try the draft version by
> anonymous checkout from the R-forge rgdal project, or using the source
> package at:
>
> http://spatial.nhh.no/R/Devel/rgdal_0.6-32.tar.gz
>
> If you report back that it suits your needs, I'll submit it to CRAN; if you
> need further changes, please let me know.
>
> Roger
>
>
>>
>> Roger
>>
>>>
>>> I've searched my disk and find gdal16 and gdal17 but not gdal18
>>> Do you know if it is possible to tell rgdal which gdal should be used?
>>>
>>> Thanks
>>>
>>> Agus
>>>
>>>
>>>
>>> 2010/12/11 Roger Bivand :
>>>>
>>>> On Fri, 10 Dec 2010, Agustin Lobo wrote:
>>>>
>>>>> Roger,
>>>>>
>>>>> I was not aware that 1.8dev is actually older than 1.7.2 I'll check why
>>>>> I
>>>>> have that gdal on the other machine and will try to fix that on my side
>>>>> on
>>>>> Monday. But perhaps this problem will come up again for future stable
>>>>> versions of 1.8.* ? Anyway, 1.8dev is clearly not a good choice and I
>>>>> have
>>>>> to change it.
>>>>
>>>> As of GDAL 1.8.0 SVN revision 21234, the problem exists. I'll fix this
>>>> in
&

[R-sig-Geo] Fwd: Problem with writing parameters in overlay()

[Agustin Lobo]

-- Forwarded message --
From: Agustin Lobo 
Date: 2010/12/16
Subject: Re: Problem with writing parameters in overlay()
To: "Robert J. Hijmans" , roger.biv...@nhh.no,
alobolis...@gmail.com


Robert,

The weird problem at writing raster layers in some cases is also
solved once I remove gdal1.8dev:
> require(rgdal)
Loading required package: rgdal
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.3, released 2010/11/10
Path to GDAL shared files: /usr/share/gdal17
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)
 require(raster)

> bN3v2 <- overlay(bN3, bSM3, 
> fun=fun1,filename="bN3v2",overwrite=TRUE,format="EHdr", NAflag=0)

works fine now.

Considering the bug fixed in gdal1.8 that I forwarded you, the problem
should also be solved even using the newest gdal1.8,
but I rather do not test until gdal1.8 is available in ubuntugis. You
might want to keep an eye on this, though.

Agus

On 11/28/2010 01:33 AM, Robert J. Hijmans wrote:
>
> Is it only with overlay, or do other functions have the same problem?
>
>  a<- calc(r2, fun=function(x)x*2, filename="test2", 
> format='GTiff',overwrite=T)
>  a<- aggregate(r2, 10, filename="test2", format='GTiff',overwrite=T)
>
> Robert
>
> On Sat, Nov 27, 2010 at 12:34 PM, Robert J. Hijmans  
> wrote:
>>
>> Hi Agus,
>>
>> This is very strange. Somehow the filetype and filename get mixed up.
>> I cannot reproduce this.
>> Does the below all work for you? (it does for me)
>>
>> library(raster)
>> r<- raster(ncol=10, nrow=10)
>> r<- writeStart(r, filename="test", format="GTiff", overwrite=TRUE)
>> filename(r)
>> r<- writeValues(r, 1:100, 1)
>> filename(r)
>> r<- writeStop(r)
>> filename(r)
>>
>> r<- raster(ncol=10, nrow=10)
>> r[]<- 1:100
>> r2<- r
>>
>> a<- overlay(r, r2, fun=sum, filename="test.tif", overwrite=TRUE, NAflag=0)
>> b<- overlay(r, r2, fun=sum, filename="test.img",overwrite=TRUE, NAflag=0)
>>
>> setOptions(todisk=TRUE)
>> a<- overlay(r, r2, fun=sum, filename="test1", format='GTiff',
>> overwrite=TRUE, NAflag=0)
>>
>> r<- writeRaster(r, filename='r1', overwrite=TRUE)
>> r2<- writeRaster(r2, filename='r1', overwrite=TRUE)
>> a<- overlay(r, r2, fun=sum, filename="test2", format='GTiff',
>> overwrite=TRUE, NAflag=0)
>> setOptions(todisk=FALSE)
>>
>>
>> Thanks for your help, Robert
>>
>>
>> On Sat, Nov 27, 2010 at 3:37 AM, Agustin Lobo  
>> wrote:
>>>
>>> Robert,
>>>
>>> There is still this problem:
>>> These work:
>>>>
>>>> bN3v2<- overlay(bN3, bSM3, fun=fun1,filename="bN3v2",NAflag=0)
>>>> bN3v2<- overlay(bN3, bSM3,
>>>> fun=fun1,filename="bN3v2",overwrite=TRUE,NAflag=0)
>>>
>>> These do not work (is not a parsing problem?)
>>> (Not critical for me, just trying to help out testing):
>>>>
>>>> bN3v2<- overlay(bN3, bSM3, fun=fun1,filename="bN3v2",overwrite=TRUE,
>>>> format="EHdr", NAflag=0)
>>>
>>> Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",
>>>  :
>>>  file: /home/alobo/EHdr does not exist
>>> Calls: overlay ... raster ->  raster ->  .local ->  .rasterObjectFromFile
>>>
>>>> bN3v2<- overlay(bN3, bSM3, fun=fun1,filename="bN3v2", format="EHdr",
>>>> NAflag=0)
>>>
>>> Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",
>>>  :
>>>  file: /home/alobo/EHdr does not exist
>>> Calls: overlay ... raster ->  raster ->  .local ->  .rasterObjectFromFile
>>>
>>>> bN3v2<- overlay(bN3, bSM3, fun=fun1,filename="bN3v2",overwrite=TRUE,
>>>> format="EHdr", NAflag=0)
>>>
>>> Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",
>>>  :
>>>  file: /home/alobo/EHdr does not exist
>>> Calls: overlay ... raster ->  raster ->  .local ->  .rasterObjectFromFile
>>>
>>>> bN3v2<- overlay(bN3, bSM3, fun=fun1,filename="bN3v2", NAflag=0,
>>>> format="GTiff")
>>>
>>> Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer"

[R-sig-Geo] sampleRandom

[Agustin Lobo]

Robert,

Could you make the output of sampleRandom() to optionally be an
spatial points data frame with the raster values across bands in the
table?
(I know it can be done currently through cells=T etc, but it would be
handy having the sppdf directly)
BTW, having row,col in addition to cell number would be handy as well

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] function parameters in overlay()

[Agustin Lobo]

Robert,

How can I define parameters in the functions within overlay() ?
For example:

> fun3 = function(x,y)
{a=-1.1867422
 s=1.0287373
 return((y-(a+s*x))/sqrt(s^2+1))}
> NmaxDev = 
> overlay(bN1999max,bN2009max,fun=fun2,filename="NmaxDev",overwrite=TRUE,format="EHdr",datatype='FLT4S',
>  NAflag=0)

want it to be something like:

> fun3 = function(x,y,a,b)
{return((y-(a+s*x))/sqrt(s^2+1))}
> NmaxDev = 
> overlay(bN1999max,bN2009max,fun=fun2(a=-1.1867422,s=1.0287373),filename="NmaxDev",overwrite=TRUE,format="EHdr",datatype='FLT4S',
>  NAflag=0)

Is it possible or should I redefine the function every time?

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] function parameters in overlay()

[Agustin Lobo]

Robert,
if possible, a syntax more similar to apply() would be more familiar
to regular R users, i.e.,
x1 = overlay(r1, r2, fun=fun2,a=-1.1867422 s=1.0287373)

I personally would find more intuitive the following:
x1 = overlay(x=r1, y=r2, fun=fun2(x,y,a=-1.1867422 s=1.0287373))
but I recognize that would be quite not a regular R style.


Agus



2010/12/17 Robert J. Hijmans :
> Agus,
>
> # self contained example...
> library(raster)
> r1 <- raster(nrow=10, ncol=10)
> r1[] = 1:ncell(r1)
> r2 <- r1
>
> # base case
> fun1 = function(x,y){ a=-1.1867422; s=1.0287373;
> return((y-(a+s*x))/sqrt(s^2+1))}
> x = overlay(r1, r2, fun=fun1 )
>
> # changing a & s without changing fun2, perhaps clumsy
> fun2 = function(x,y,a,s) {return((y-(a+s*x))/sqrt(s^2+1))}
> x1 = overlay(r1, r2, fun=function(x,y){fun2 (x, y, a=-1.1867422,s=1.0287373)})
> x2 = overlay(r1, r2, fun=function(...){fun2 (..., a=-1.1867422,s=1.0287373)})
>
> # changing a & s without changing fun3, perhaps cleaner
> # I believe this works because a and s are defined in the same
> environment as fun3
> fun3 = function(x,y) { return((y-(a+s*x))/sqrt(s^2+1)) }
> a = -1.1867422
> s = 1.0287373
> x3 = overlay(r1, r2, fun=fun3)
>
> Best, Robert
>
> On Fri, Dec 17, 2010 at 1:04 PM, Agustin Lobo  wrote:
>> Robert,
>>
>> How can I define parameters in the functions within overlay() ?
>> For example:
>>
>>> fun3 = function(x,y)
>> {a=-1.1867422
>>  s=1.0287373
>>  return((y-(a+s*x))/sqrt(s^2+1))}
>>> NmaxDev = 
>>> overlay(bN1999max,bN2009max,fun=fun2,filename="NmaxDev",overwrite=TRUE,format="EHdr",datatype='FLT4S',
>>>  NAflag=0)
>>
>> want it to be something like:
>>
>>> fun3 = function(x,y,a,b)
>> {return((y-(a+s*x))/sqrt(s^2+1))}
>>> NmaxDev = 
>>> overlay(bN1999max,bN2009max,fun=fun2(a=-1.1867422,s=1.0287373),filename="NmaxDev",overwrite=TRUE,format="EHdr",datatype='FLT4S',
>>>  NAflag=0)
>>
>> Is it possible or should I redefine the function every time?
>>
>> Agus
>>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Converting array into raster brick

[Agustin Lobo]

Robert,

Converting from a regular 3D array to brick is probably going to be
quite a  common operation,
perhaps you could consider a simpler syntax in the future, i,e:

r = as.brick(x)

where x would be a 3D array?

Agus

2010/12/16 Robert J. Hijmans :
> Christian,
>
> Could be quicker if you make a matrix out of the array (or, if
> possible, avoid the array in the first place). See below
>
> library(raster)
>
> # create an array
> m <- matrix(1:100, ncol=10, byrow=TRUE)
> a <- array(list(m, m, m, m, m))
>
> # empty RasterBrick
> x <- brick(ncol=10, nrow=10)
>
> # make matrix from array
> v <- sapply(a, function(x) as.vector(t(x)))
>
> x <- setValues(x, v)
> plot(x)
>
>
> Best, Robert
>
> On Thu, Dec 16, 2010 at 12:26 PM,   wrote:
>> Dear all,
>>
>> I needed to convert an array into a RasterBrick object. The output did not
>> need to be written to memory, as assessed by canProcessInMemory(). The only
>> way to put the data into the RasterBrick that worked for me was layer by
>> layer:
>>
>> for (l in 1:nlayers(x)) x
>> <-setValues(x,values=as.vector(t(xx[,,l])),layer=l)
>>
>> This was quite complicated and slow. So I am wondering if there is a better
>> way of doing this.
>>
>> Many thanks, Christian
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Help with rgdal on MACOS

[Agustin Lobo]

Dea all,

Quite a number of students will be using MAC for a course in which we
will make use
of rgdal and raster. As I do not have an easy access to MAC machines
and I've found quite
a number of different options, I'm seeking some clarifications.

Students will have gdal and other dependencies installed from
kyngchaos as they should have QGIS installed.
Under these circumstances, which alternative is better:

1. As suggested by Roger in
http://r-sig-geo.2731867.n2.nabble.com/rgdal-on-MAC-td5073766.html
> setRepositories(ind=1:2)
> install.packages('rgdal')

2. Install rgdal from source.
2.1. As stated at the end of
http://spatialanalysis.co.uk/2010/11/02/installing-rgdal-on-mac-os-x/
> R CMD INSTALL
 
–configure-args=”–with-gdal-config=/Library/Frameworks/GDAL.framework/unix/bin/gdal-config
 –with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib
 –with-proj-include=/Library/Frameworks/PROJ.framework/unix/include”
rgdal_0.6-20.tar.gz
2.2. Using homebrew:
brew install gdal
R CMD INSTALL –configure-args=’–with-proj-include=/usr/local/lib’
rgdal_0.6-28.tar.gz
(in this case, is the brew command using the kyngchaos gdal?)

3. Other alternatives

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Help with rgdal on MACOS

[Agustin Lobo]

Actually, more alternatives:

3. Suggested by Ralf Strasse:
Install Xcode tools from Apple (is recommended anyway for many
applications). This will install several compilers.Afterwards
installing from source should work fine.
http://developer.apple.com/technologies/tools/xcode.html
And also select on CRAN under Mac OS X installation "Tools" and
installe gfortran (and the other two - you may need them another
time).
install.packages("rgdal",
configure.args='--with-gdal-config=/Library/Frameworks/GDAL.Framework/Unix/bin/gdal-config
--with-proj-include=/Library/Frameworks/PROJ.framework/unix/include-with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib/',dependencies=TRUE,
type="source")

4. Downloading binary from kyngchaos
http://www.kyngchaos.com/software/frameworks?s[]=rgdal
and installing with
install.packages("~/Downloads/rgdal-0.6.29-1.zip")
(I'm guessing the directory, is this the downloading directory by
default in Macs?)

Any advice on which one is the best alternative?

Thanks

Agus

2010/12/23 Agustin Lobo :
> Dea all,
>
> Quite a number of students will be using MAC for a course in which we
> will make use
> of rgdal and raster. As I do not have an easy access to MAC machines
> and I've found quite
> a number of different options, I'm seeking some clarifications.
>
> Students will have gdal and other dependencies installed from
> kyngchaos as they should have QGIS installed.
> Under these circumstances, which alternative is better:
>
> 1. As suggested by Roger in
> http://r-sig-geo.2731867.n2.nabble.com/rgdal-on-MAC-td5073766.html
>> setRepositories(ind=1:2)
>> install.packages('rgdal')
>
> 2. Install rgdal from source.
> 2.1. As stated at the end of
> http://spatialanalysis.co.uk/2010/11/02/installing-rgdal-on-mac-os-x/
>> R CMD INSTALL
>  –configure-args=”–with-gdal-config=/Library/Frameworks/GDAL.framework/unix/bin/gdal-config
>  –with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib
>  –with-proj-include=/Library/Frameworks/PROJ.framework/unix/include”
> rgdal_0.6-20.tar.gz
> 2.2. Using homebrew:
> brew install gdal
> R CMD INSTALL –configure-args=’–with-proj-include=/usr/local/lib’
> rgdal_0.6-28.tar.gz
> (in this case, is the brew command using the kyngchaos gdal?)
>
> 3. Other alternatives
>
> Thanks
>
> Agus
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] R-sig-Geo Digest, Vol 88, Issue 23

[Agustin Lobo]

Many thanks for your clarifications, Ralf.
But if
> Actually the 4 options more or less boil down to two if you check what is 
> done in detail:
.../...

then  we are just discarding alternative 1 in my first post:
"> setRepositories(ind=1:2)
> install.packages('rgdal')
"
Also, what about homebrew? is this not a good way to make this type of
things easier on Macs?

Finally, where you say
> BTW: "Linux Techies"  should have no major probs with Mac, given that the Mac 
> system commands are more or less the same as for Linux systems
The problem is that most mac users I've had as students are not
Techies on anything. Otherwise, Macs are wonderful machines.

Agus

2010/12/23 Ralf Schäfer :
> Just a quick reply:
>
> Actually the 4 options more or less boil down to two if you check what is 
> done in detail:
> 1) install from source and install external dependencies before.
> This option should always be mentioned in case binaries do not work on a 
> specific system (Tiger, Panther, Leopard etc) - and apart from that it can 
> take a bit after rgdal binaries are available on kyngchaos after R updates 
> (and some people may not want to wait after future updates). In addition, the 
> installed tools may also be necessary for the R commander GUI, GIMP and other 
> programs or if you need to install other packages from source.
> So the way to go is:
> 1. Install GDAL Framework on kyngchaos GDAL 1.7 Complete
> 2. Install R AND all tools for Mac: Gfortran gfortran-4.2.3.dmg TclTk 
> tcltk-8.5.5-x11.dmg Developer pack devpack4-darwin8-bin4.tar.gz
> 3. Install Xcode tools - either from DVD that ships with every mac or after 
> free registration here http://developer.apple.com/programs/register/ (the 
> Xcode tools are also recommended in the Mac OS X FAQ for R 
> http://cran.r-project.org/bin/macosx/RMacOSX-FAQ.html)
> 4. Install rgdal as already outlined
> install.packages("rgdal",
> configure.args='--with-gdal-config=/Library/Frameworks/GDAL.Framework/Unix/bin/gdal-config
> --with-proj-include=/Library/Frameworks/PROJ.framework/unix/include-with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib/',dependencies=TRUE,
> type="source")
>
> 2) as outlined by Barry
>> 1. Download the kyngchaos GDAL Complete 'Framework':
>> http://www.kyngchaos.com/files/software/frameworks/GDAL_Complete-1.7.dmg
>> 2. Double click it, I think that opened it up and then somehow things
>> from that get installed. Wish I'd kept screenshots
>> 2.5 Wonder where in the .dmg file the rgdal binary is.
>> 3. Download the kyngchaos rgdal binary:
>> http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
>> and did install.packages on that.
>
>
>
> In general option 2  is easiest - but on the long run option 1) will help you 
> also with other packages (Can`t remember exactly but just installed a new Mac 
> Pro and had to install a couple of R packages from source).
> BTW: "Linux Techies"  should have no major probs with Mac, given that the Mac 
> system commands are more or less the same as for Linux systems
>
> Have a nice Xmas
> Ralf
>
>
>> Date: Thu, 23 Dec 2010 08:35:13 +
>> From: Barry Rowlingson 
>> To: agustin.l...@ija.csic.es
>> Cc: r-sig-geo 
>> Subject: Re: [R-sig-Geo] Help with rgdal on MACOS
>> Message-ID:
>>       
>> Content-Type: text/plain; charset=ISO-8859-1
>>
>> On Thu, Dec 23, 2010 at 7:36 AM, Agustin Lobo  wrote:
>>
>>> 4. Downloading binary from kyngchaos
>>> http://www.kyngchaos.com/software/frameworks?s[]=rgdal
>>> and installing with
>>> install.packages("~/Downloads/rgdal-0.6.29-1.zip")
>>> (I'm guessing the directory, is this the downloading directory by
>>> default in Macs?)
>>>
>>> Any advice on which one is the best alternative?
>>>
>>
>> I was recently asked by a new student with a Mac to sort out how to
>> read shapefiles into R for her. I know nothing about Mac packaging,
>> binaries, compilers etc, and she knew even less.
>>
>> After reading far too many confusing 'how to get rgdal working on a
>> Mac' mailing list messages, I got it working. What I think I did was:
>>
>> 1. Download the kyngchaos GDAL Complete 'Framework':
>> http://www.kyngchaos.com/files/software/frameworks/GDAL_Complete-1.7.dmg
>> 2. Double click it, I think that opened it up and then somehow things
>> from that get installed. Wish I'd kept screenshots
>> 2.5 Wonder where in the .dmg file the rgdal binary is.
>> 3. Download the kyngchaos rgdal binary:
>> http://www.kyngchaos.com/

[R-sig-Geo] MACOSX: rgdal cannot find GDAL

[Agustin Lobo]

After having installed GDAL Complete from
http://www.kyngchaos.com/files/software/frameworks/GSL_Framework-1.14-1.dmg,

I try gdalinfo but:

imac-de-agustin-lobo:~ agustinlobo$ gdalinfo --version
-bash: gdalinfo: command not found

Also Spotlight cannot find but the uncompressed folder GDAL Complete,
or any of the installed
programs (i.e. gdalinfo or gdal-config)

imac-de-agustin-lobo:~ agustinlobo$ gdal-config-bash: gdal-config:
command not found

Obviously, rgdal cannot find GDAL

> setRepositories(ind=1:2)
> install.packages("rgdal")
--- Please select a CRAN mirror for use in this session ---
probando la URL
'http://www.stats.ox.ac.uk/pub/RWin/bin/macosx/leopard/contrib/2.12/rgdal_0.6-31.tgz'
Content type 'application/x-gzip' length 10203689 bytes (9.7 Mb)
URL abierta
==
downloaded 9.7 Mb


The downloaded packages are in

/var/folders/iR/iRxOvaG9F10viNoa7WceNU+++TI/-Tmp-//RtmpLu8kR8/downloaded_packages
> require(rgdal)
Loading required package: rgdal
Loading required package: sp
Error in dyn.load(file, DLLpath = DLLpath, ...) :
  unable to load shared object
'/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so':
  
dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so,
6): Library not loaded: libR.dylib
  Referenced from:
/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so
  Reason: image not found

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] rgdal from kyngchaos: Error in get(Info[i, 1], envir = env) :

[Agustin Lobo]

After having installed GDAL from
http://www.kyngchaos.com/files/software/frameworks/GDAL_Complete-1.7.dmg
and downloading
http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
I install rgdal and get:

> install.packages("/Users/agustinlobo/Downloads/rgdal/rgdal_0.6-29.tgz")
inferring 'repos = NULL' from the file name
> require(rgdal)
Loading required package: rgdal
Error in get(Info[i, 1], envir = env) :
  no se puede asignar un bloque de memoria de tamaño 2.2 Gb

(== cannot assign a memory block of size 2.2 Gb)

Is this a particular problem in this computer? Any fix?

(BTW, cannot find how to get these messages in English)

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] MACOSX: rgdal cannot find GDAL

[Agustin Lobo]

Ralf,

2010/12/24 Ralf Schäfer :
> Hello,
> your download link referred to  another framework. Check whether you have
> really installed this one:
> GDAL 1.7 Complete
Sorry, my mistake, but  just at writing the link, I actually downloaded
the correct one:
http://www.kyngchaos.com/files/software/frameworks/GDAL_Complete-1.7.dmg

> Can you allocate the folder:
> /Library/Frameworks/GDAL.framework 
Yes:
imac-de-agustin-lobo:Frameworks agustinlobo$ pwd
/Library/Frameworks
imac-de-agustin-lobo:Frameworks agustinlobo$ ls -l
total 0
drwxrwxrwx  5 root  wheel  170 10 sep  2008 Adobe AIR.framework
drwxrwxr-x  8 root  admin  272 21 dic 17:44 FreeType.framework
drwxrwxr-x  8 root  admin  272 24 dic 09:13 GDAL.framework
drwxrwxr-x  8 root  admin  272 24 dic 09:13 GEOS.framework
drwxrwxr-x  8 root  admin  272 21 dic 17:58 GSL.framework
drwxrwxr-x  8 root  admin  272 12 dic  2009 HPDeviceModel.framework
drwxrwxr-x  6 root  admin  204 12 dic  2009 HPPml.framework
drwxrwxr-x  7 root  admin  238 12 dic  2009 HPServicesInterface.framework
drwxrwxr-x  7 root  admin  238 12 dic  2009 HPSmartPrint.framework
drwxr-xr-x  5 1422  wheel  170 12 dic  2009 HPSmartX.framework
drwxr-xr-x  6 root  admin  204 27 ene  2010 MacFUSE.framework
drwxrwxr-x  8 root  admin  272 24 dic 09:13 PROJ.framework
drwxrwxr-x  8 root  admin  272 24 dic 08:59 R.framework
drwxrwxr-x  8 root  admin  272 24 dic 09:13 SQLite3.framework
drwxrwxr-x  8 root  admin  272 24 dic 09:13 UnixImageIO.framework
drwxrwxr-x  8 root  admin  272 21 dic 17:56 cairo.framework

This problem persists:
imac-de-agustin-lobo:GDAL.framework agustinlobo$ gdalinfo --version
-bash: gdalinfo: command not found

But actually, the R installation of rgdal works now. Could it be that
I needed to reboot after
installing GDAL? I'd like to clarify as I'm trying to write up directions:

> install.packages("/Users/agustinlobo/Downloads/rgdal/rgdal_0.6-29.tgz")
inferring 'repos = NULL' from the file name
> require(rgdal)
Loading required package: rgdal
Loading required package: sp
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.7.3, released 2010/11/10
Path to GDAL shared files:
/Library/Frameworks/GDAL.framework/Versions/1.7/Resources/gdal
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009
Path to PROJ.4 shared files: (autodetected)

>
> But you obviously fixed that already as I can read from the second post?!

No, those are different problems. The one I report here is the installation with
the rgdal provided in kyngchaos
http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip

The other one is the problem installing rgdal from CRAN repositories 1 and 2

I still do not know which one is to be preferred...

Agus



> Regards
> Ralf
>
>
>
> ------
>
> Message: 11
> Date: Fri, 24 Dec 2010 10:21:13 +0100
> From: Agustin Lobo 
> To: r-sig-geo 
> Subject: [R-sig-Geo] MACOSX: rgdal cannot find GDAL
> Message-ID:
> 
> Content-Type: text/plain; charset=ISO-8859-1
>
> After having installed GDAL Complete from
> http://www.kyngchaos.com/files/software/frameworks/GSL_Framework-1.14-1.dmg,
>
> I try gdalinfo but:
>
> imac-de-agustin-lobo:~ agustinlobo$ gdalinfo --version
> -bash: gdalinfo: command not found
>
> Also Spotlight cannot find but the uncompressed folder GDAL Complete,
> or any of the installed
> programs (i.e. gdalinfo or gdal-config)
>
> imac-de-agustin-lobo:~ agustinlobo$ gdal-config-bash: gdal-config:
> command not found
>
> Obviously, rgdal cannot find GDAL
>
> setRepositories(ind=1:2)
>
> install.packages("rgdal")
>
> --- Please select a CRAN mirror for use in this session ---
> probando la URL
> 'http://www.stats.ox.ac.uk/pub/RWin/bin/macosx/leopard/contrib/2.12/rgdal_0.6-31.tgz'
> Content type 'application/x-gzip' length 10203689 bytes (9.7 Mb)
> URL abierta
> ==
> downloaded 9.7 Mb
>
>
> The downloaded packages are in
> /var/folders/iR/iRxOvaG9F10viNoa7WceNU+++TI/-Tmp-//RtmpLu8kR8/downloaded_packages
>
> require(rgdal)
>
> Loading required package: rgdal
> Loading required package: sp
> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>  unable to load shared object
> '/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so':
>  dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so,
> 6): Library not loaded: libR.dylib
>  Referenced from:
> /Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so
>  Reason: image not found
>
>
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] rgdal from kyngchaos: Error in get(Info[i, 1], envir = env) :

[Agustin Lobo]

Sorry, I was mixing installations. The problem I report here is the
one with the
rgdal from CRAN repositories. This is the correct story:

After having installed GDAL from
http://www.kyngchaos.com/files/software/frameworks/GDAL_Complete-1.7.dmg

> remove.packages("rgdal")
Removing package(s) from ‘/Library/Frameworks/R.framework/Resources/library’
(as ‘lib’ is unspecified)
Error en .find.package(pkgs, lib) : there is no package called 'rgdal'
> setRepositories(ind=1:2)
> install.packages("rgdal")
probando la URL
'http://www.stats.ox.ac.uk/pub/RWin/bin/macosx/leopard/contrib/2.12/rgdal_0.6-31.tgz'
Content type 'application/x-gzip' length 10203689 bytes (9.7 Mb)
URL abierta
==
downloaded 9.7 Mb


The downloaded packages are in

/var/folders/iR/iRxOvaG9F10viNoa7WceNU+++TI/-Tmp-//Rtmpx5JXhi/downloaded_packages
> require(rgdal)
Loading required package: rgdal
Error in dyn.load(file, DLLpath = DLLpath, ...) :
  unable to load shared object
'/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so':
  
dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so,
6): Library not loaded: libR.dylib
  Referenced from:
/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so
  Reason: image not found

I do have rgdal.so in that directory:
imac-de-agustin-lobo:i386 agustinlobo$ pwd
/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386
imac-de-agustin-lobo:i386 agustinlobo$ ls -l
total 15088
-rwxr-xr-x  1 agustinlobo  admin  7722168  4 dic 15:57 rgdal.so


I can get rgdal running using the version from kyngchaos, but not with this one.

Thanks,

Agus


2010/12/24 Agustin Lobo :
> After having installed GDAL from
> http://www.kyngchaos.com/files/software/frameworks/GDAL_Complete-1.7.dmg
> and downloading
> http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
> I install rgdal and get:
>
>> install.packages("/Users/agustinlobo/Downloads/rgdal/rgdal_0.6-29.tgz")
> inferring 'repos = NULL' from the file name
>> require(rgdal)
> Loading required package: rgdal
> Error in get(Info[i, 1], envir = env) :
>  no se puede asignar un bloque de memoria de tamaño 2.2 Gb
>
> (== cannot assign a memory block of size 2.2 Gb)
>
> Is this a particular problem in this computer? Any fix?
>
> (BTW, cannot find how to get these messages in English)
>
> Agus
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] MACOSX: rgdal cannot find GDAL

[Agustin Lobo]

Ralf,

Setting the path as you indicate solves the problem of running gdalinfo:
imac-de-agustin-lobo:~ agustinlobo$ export
PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH
imac-de-agustin-lobo:~ agustinlobo$ echo
$PATH/Library/Frameworks/GDAL.framework/Programs:/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin:/usr/X11/bin
imac-de-agustin-lobo:~ agustinlobo$ gdalinfo --version
GDAL 1.7.3, released 2010/11/10

Therefore we can close this thread on installing GDAL and rgdal on
MACOSX from the kyngchaos server!
I'll be sending the link with directions for dummies shortly

Thank you all

Agus

2010/12/27 Ralf Schäfer :
> Hi Agustin,
>
>>> Can you allocate the folder:
>>> /Library/Frameworks/GDAL.framework 
>> Yes:
>> imac-de-agustin-lobo:Frameworks agustinlobo$ pwd
>> /Library/Frameworks
>> This problem persists:
>> imac-de-agustin-lobo:GDAL.framework agustinlobo$ gdalinfo --version
>> -bash: gdalinfo: command not found
>> But actually, the R installation of rgdal works now. Could it be that
>> I needed to reboot after
>> installing GDAL? I'd like to clarify as I'm trying to write up directions:
>
> I am not a Shell expert - but my guess would be that the GDAL commands are 
> not in your PATH. So when you type the commands Shell does not find them.
> Try to include them in your Path:
> export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH
> I actually think that I had to do this as well...
>
> In case that works you can write this path to your profile so you have it for 
> every terminal session - here is a short guide:
> http://www.macgasm.net/2008/04/10/adding-a-new-location-to-your-path-variable-within-terminal/
>
>>
>> The other one is the problem installing rgdal from CRAN repositories 1 and 2
>
> Unfortunately, I have no idea why this fails for you. In fact, I checked it 
> as well, and did not work for me either. I get the same error even when I 
> specify GDAL and Proj path:
> install.packages("rgdal", 
> configure.args='--with-gdal-config=/Library/Frameworks/GDAL.Framework/unix/bin/gdal-config
>  
> --with-proj-include=/Library/Frameworks/PROJ.framework/unix/include-with-proj-lib=/Library/Frameworks/PROJ.framework/unix/lib/')
>
> So that is why I also installed from source as outlined previously- and I 
> noticed that on the CRAN source there is a later version of rgdal 0.6-33 
> instead of 0.6-31, may be this is the problem (I updated to R 2.12.1).
> I reckon other people can be of more help here.
>
>> I still do not know which one is to be preferred...
>
> Well, as Roger Bivand pointed out (24.12)  Installing from CRAN/CRAN extra is 
> certainly the best option. May be there is a fix to install the binary from 
> CRAN/CRAN extra
> Otherwise certainly the kyngchaos binaries are the most convenient way (yet 
> they are rgdal 0.6-29)
>
> Cheers
> Ralf
>        [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Installation of rgdal from CRAN on MACOSX does not work

[Agustin Lobo]

While the installation from
http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
works fine
the installation of rgdal from CRAN on MACOSX does not work, despite
my previous message that this issue was solved as well
(my confusion was a consequence of mixing both installations in the
same session):

> setRepositories(ind=1:2)
> install.packages("rgdal")
--- Please select a CRAN mirror for use in this session ---
probando la URL
'http://www.stats.ox.ac.uk/pub/RWin/bin/macosx/leopard/contrib/2.12/rgdal_0.6-31.tgz'
Content type 'application/x-gzip' length 10203689 bytes (9.7 Mb)
URL abierta
==
downloaded 9.7 Mb


The downloaded packages are in

/var/folders/iR/iRxOvaG9F10viNoa7WceNU+++TI/-Tmp-//RtmpqIhqn9/downloaded_packages
> require(rgdal)
Loading required package: rgdal
Loading required package: sp
Error in dyn.load(file, DLLpath = DLLpath, ...) :
  unable to load shared object
'/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so':
  
dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so,
6): Library not loaded: libR.dylib
  Referenced from:
/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so
  Reason: image not found
> sessionInfo()
R version 2.12.1 (2010-12-16)
Platform: i386-apple-darwin9.8.0/i386 (32-bit)

locale:
[1] es_ES.UTF-8/es_ES.UTF-8/C/C/es_ES.UTF-8/es_ES.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] sp_0.9-76

loaded via a namespace (and not attached):
[1] grid_2.12.1 lattice_0.19-13 tools_2.12.1

Note that rgdal.so is where it should be:

imac-de-agustin-lobo:~ agustinlobo$ ls -l
/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/
total 15088
-rwxr-xr-x  1 agustinlobo  admin  7722168  4 dic 15:57 rgdal.so

This error is not solved by setting the GDAL path.

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Installation of rgdal from CRAN on MACOSX does not work

[Agustin Lobo]

Roger,

But
http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
installs with no problems on the same machine and on the same R, would that not
discard the possibility of an
> issue of your R installation,
> possibly a multi-architecture issue

as you suggest?

Please note I'm not personally in a hurry any more as we can work with
the rgdal from kyngchaos. I just want to report the error to help
having an rgdal installation for MACOSX as
simple as for Windows.
I understand installing from source will always be a bit more
difficult on MACOSX than on linux as there are some basic tools and
libraries that are missing on regular MAC systems that must be
previously installed, but this is not my concern now as I'm focusing
on installations from binaries.

Thanks!

Agus



2010/12/29 Roger Bivand :
> On Wed, 29 Dec 2010, Agustin Lobo wrote:
>
>> While the installation from
>> http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
>> works fine
>> the installation of rgdal from CRAN on MACOSX does not work, despite
>> my previous message that this issue was solved as well
>
> Please do read the error message:
>
>>  unable to load shared object
>>
> '/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i$
>>
> dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal$
>>
>> 6): Library not loaded: libR.dylib
>>  Referenced from:
>>
> /Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i3$
>>
>>  Reason: image not found
>
> The image that is not found is libR.dylib, the main R dylib against which
> rgdal.so is linked. So this is still an issue of your R installation,
> possibly a multi-architecture issue (is libR.dylib the i386 version,
> guessing - is libR64.dylib the 64-bit Intel version?).
>
> Could any OSX users who are using CRAN extra rgdal please help Agus to
> resolve this - it could be that local path modifications are hiding
> libR.dylib, or that R was installed 64-bit only, and the INSTALL_opts=
> argument to install.packages() was not set correctly? I assume that the CRAN
> extras OSX rgdal binary is built against a standard R OSX Intel install,
> which is multi-architecture.
>
> For a fully updated rgdal, the best combination remains installing the full
> GDAL framework from Kyngchaos, the full tools collection:
>
> http://r.research.att.com/tools/
> http://cran.r-project.org/bin/macosx/tools
>
> and install from CRAN source. For a lab, the instructor would rather build a
> local OSX rgdal binary package using the tools collection on one development
> machine with all the correct path variables set (R CMD check rgdal; R CMD
> INSTALL --build rgdal), paying attention to multi-architectures, and
> distribute to lab machines from there. This means that the lab instructor
> can keep rgdal updated if need be, as both Kyngchaos and CRAN extras rgdal
> OSX binary packages typically lag CRAN source.
>
> Please also note that there are very many very old GDAL binaries out there
> that users pick up and apply, missing all the subsequent bug-fixes and
> improvements - GDAL 1.8.0 release candidate 1 is now out, so we are looking
> at an upgrade to 1.8.0 during January, which Kyngchaos will most likely pick
> up quickly, probably with a current rgdal binary.
>
> Roger
>
>
>> (my confusion was a consequence of mixing both installations in the
>> same session):
>>
>>> setRepositories(ind=1:2)
>>> install.packages("rgdal")
>>
>> --- Please select a CRAN mirror for use in this session ---
>> probando la URL
>>
>> 'http://www.stats.ox.ac.uk/pub/RWin/bin/macosx/leopard/contrib/2.12/rgdal_0.6-31.tgz'
>> Content type 'application/x-gzip' length 10203689 bytes (9.7 Mb)
>> URL abierta
>> ==
>> downloaded 9.7 Mb
>>
>>
>> The downloaded packages are in
>>
>>  /var/folders/iR/iRxOvaG9F10viNoa7WceNU+++TI/-Tmp-//RtmpqIhqn9/downloaded_packages
>>>
>>> require(rgdal)
>>
>> Loading required package: rgdal
>> Loading required package: sp
>> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>>  unable to load shared object
>>
>> '/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so':
>>
>>  dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so,
>> 6): Library not loaded: libR.dylib
>>  Referenced from:
>>
>> /Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so
>>  Reason: image not found
&g

Re: [R-sig-Geo] Installation of rgdal from CRAN on MACOSX does not work

[Agustin Lobo]

Thanks,

it has worked with no problems on 4 Macs, both Snow Leopard and Leopard.

Agus


2011/1/10 Roger Bivand :
> On Sat, 8 Jan 2011, Roger Bivand wrote:
>
>> On Thu, 30 Dec 2010, Roger Bivand wrote:
>>
>>> On Wed, 29 Dec 2010, Agustin Lobo wrote:
>>>
>>>> Roger,
>>>>
>>>> But
>>>> http://www.kyngchaos.com/files/software/frameworks/rgdal-0.6.29-1.zip
>>>> installs with no problems on the same machine and on the same R, would
>>>> that not
>>>> discard the possibility of an
>>>>>
>>>>> issue of your R installation,
>>>>> possibly a multi-architecture issue
>>>>
>>>> as you suggest?
>>>>
>>>> Please note I'm not personally in a hurry any more as we can work with
>>>> the rgdal from kyngchaos.
>>>
>>> I will have access to an OSX laptop in a little over a week, and will try
>>> out the CRAN rgdal binaries. Just one report of installation failure isn't
>>> enough. I think that maybe someone on R-sig-mac might know why rgdal.so is
>>> looking for libR.dylib; I don't use OSX, and don't have easy access to such
>>> a machine, which is why I asked for others to comment.
>>
>> On a Leopard 10.5.8, I see the problem reported by Agus for the CRAN-extra
>> OSX binary rgdal package and freshly installed R 2.12.1 (libR.dylib image
>> not found) in the 32-bit R-GUI, but no problem when running R from the
>> command line. A different problem (mach-o, but wrong architecture) occurs in
>> 64-bit systems, both run from the command line and from R-GUI. Could anyone
>> who wants to follow this up please post to R-sig-Mac, which is the correct
>> list for OSX issues?
>>
>> Conclusion: on Leopard, the CRAN extras rgdal binary can be used in the
>> 32-bit R application run from the command line, but not otherwise. R 2.12.1
>> installs the 64-bit GUI into the OSX applications toolbar - one needs to
>> open a terminal and type "R", checking that the length of a pointer is 4, to
>> run the version that works.
>
> Prof. Ripley has updated the OSX CRAN extras rgdal binaries to rgdal 0.6.33,
> and this new version is confirmed to work under Leopard both 32 and 64 bit,
> both GUI and terminal. So the simplest way to install rgdal for OSX for the
> basic set of GDAL and OGR drivers is once again:
>
> setRepositories(ind=1:2)
> install.packages("rgdal")
>
> Roger
>
>>
>> Roger
>>
>> PS. The symptoms are:
>>
>> R64, R-GUI 64
>>
>>> library(rgdal)
>>
>> Loading required package: sp
>> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>>  unable to load shared object
>> '/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/x86_64/rgdal.so':
>>
>>
>> dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/x86_64/rgdal.so,
>> 6): no suitable image found.  Did find:
>>
>>
>> /Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/x86_64/rgdal.so:
>> mach-o, but wrong architecture
>> Error: package/namespace load failed for 'rgdal'
>>
>>
>>
>> R-GUI 32
>>
>>> library(rgdal)
>>
>> Loading required package: sp
>> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>>  unable to load shared object
>> '/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so':
>>
>>
>> dlopen(/Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so,
>> 6): Library not loaded: libR.dylib
>>  Referenced from:
>> /Library/Frameworks/R.framework/Versions/2.12/Resources/library/rgdal/libs/i386/rgdal.so
>>  Reason: image not found
>> Error: package/namespace load failed for 'rgdal'
>>
>>
>>>
>>> Roger
>>>
>>>> I just want to report the error to help
>>>> having an rgdal installation for MACOSX as
>>>> simple as for Windows.
>>>> I understand installing from source will always be a bit more
>>>> difficult on MACOSX than on linux as there are some basic tools and
>>>> libraries that are missing on regular MAC systems that must be
>>>> previously installed, but this is not my concern now as I'm focusing
>>>> on installations from binaries.
>>>>
>>>> Thanks!
>>>>
>>>> Agus
>>>>
>>>>
>>>>
>>>> 2010

[R-sig-Geo] listw for local Moran in spdep

[Agustin Lobo]

Hi!

I'd like to run localmoran() from spdep on a small grid (matrix). Is
there a function to calculate the spatial weights (listw) from
the grid  (defining 3x3 neighborhoods)?

(BTW: the definition of LISA in the spdep manual is unformatted)

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] listw for local Moran in spdep

[Agustin Lobo]

Thanks. I have problems to understand the order in which I have to
input the data from the matrix.
This is what I'm doing:

require(spdep)
gpclibPermit()
require(ncf)
#data
x <- expand.grid(1:20, 1:5)[,1]
y <- expand.grid(1:20, 1:5)[,2]
zori <- rmvn.spa(x=x, y=y, p=2, method="exp")
z = scale(zori)
dim(z)=c(20,5)
z = t(z)[5:1,]
rz = raster(z)

#weights
wrz = cell2nb(nrow=nrow(rz), ncol=ncol(rz), type="queen", torus=FALSE)
#lisa
lisarz = localmoran(x=rz@data@values, listw=nb2listw(wrz))

(but I do not know if the ordering of the data x and the weights listw
are matching)

l3r = lisarz[,1]
dim(l3r)=c(20,5)
l3r = t(l3r)[5:1,]
l3r = raster(l3r)
plot(l3r)



Agus

2011/2/7 Roger Bivand :
> On Mon, 7 Feb 2011, Agustin Lobo wrote:
>
>> Hi!
>>
>> I'd like to run localmoran() from spdep on a small grid (matrix). Is
>> there a function to calculate the spatial weights (listw) from
>> the grid  (defining 3x3 neighborhoods)?
>
> For a 3x3 neighbourhood on a grid simply defined by numbers of rows and
> columnd, see ?cell2nb. Otherwise, use dnearneigh() with the threshold set to
> the length of the diagonal distance between cell centres giving the cell
> centre coordinates as input.
>
>>
>> (BTW: the definition of LISA in the spdep manual is unformatted)
>>
>
> The LaTeX version looks fine! I'll try to repair the regular page.
>
> Roger
>
>> Thanks
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: roger.biv...@nhh.no
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster: extract is empty with polygon

[Agustin Lobo]

Hi!
I'm trying to calculate the mean values for a set of polygons.
This is what I do:

require(raster)
require(rgdal)
#SGRGBWBPS125F40
#"read" tif file
SGRGBF40 = 
brick("/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
> SGRGBF40
class   : RasterBrick
filename: 
/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif
nlayers : 3
nrow: 1760
ncol: 2640
ncell   : 4646400
projection  : NA
min value   : 0 0 0
max value   : 65535 65535 65535
extent  : 0, 2640, 0, 1760  (xmin, xmax, ymin, ymax)
resolution  : 1, 1  (x, y)
> projection(SGRGBF40)
[1] "NA"

#read the shape file
calibf2 = 
readOGR(dsn="/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
> projection(calibf2)
[1] "+proj=utm +zone=31 +ellps=intl +units=m +no_defs"

> projection(SGRGBF40) = projection(calibf2)

But they do not overlap as shown by
> plot(subset(SGRGBF40,1))
> plot(calibf2,add=T)

and the extracted object is empty
#Calculate mean values for each polygon:
> v <- extract(SGRGBF40, calibf2,fun=mean,nl=3)
> summary(v)
  Length Class  Mode
 [1,] 0  -none- NULL
 [2,] 0  -none- NULL
 [3,] 0  -none- NULL
 [4,] 0  -none- NULL
 [5,] 0  -none- NULL
etc

> sessionInfo()
R version 2.12.2 (2011-02-25)
Platform: i486-pc-linux-gnu (32-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8   LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=C  LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8   LC_NAME=C
 [9] LC_ADDRESS=C   LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.6-31 raster_1.7-8 sp_0.9-72

loaded via a namespace (and not attached):
[1] grid_2.12.2 lattice_0.19-17 tools_2.12.2

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster: extract is empty with polygon

[Agustin Lobo]

I forgot to add that raster and vector do overlap in qgis (the vector
was actually digitized on top of the raster)
Agus

2011/3/24 Agustin Lobo :
> Hi!
> I'm trying to calculate the mean values for a set of polygons.
> This is what I do:
>
> require(raster)
> require(rgdal)
> #SGRGBWBPS125F40
> #"read" tif file
> SGRGBF40 = 
> brick("/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
>> SGRGBF40
> class       : RasterBrick
> filename    : 
> /media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif
> nlayers     : 3
> nrow        : 1760
> ncol        : 2640
> ncell       : 4646400
> projection  : NA
> min value   : 0 0 0
> max value   : 65535 65535 65535
> extent      : 0, 2640, 0, 1760  (xmin, xmax, ymin, ymax)
> resolution  : 1, 1  (x, y)
>> projection(SGRGBF40)
> [1] "NA"
>
> #read the shape file
> calibf2 = 
> readOGR(dsn="/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
>> projection(calibf2)
> [1] "+proj=utm +zone=31 +ellps=intl +units=m +no_defs"
>
>> projection(SGRGBF40) = projection(calibf2)
>
> But they do not overlap as shown by
>> plot(subset(SGRGBF40,1))
>> plot(calibf2,add=T)
>
> and the extracted object is empty
> #Calculate mean values for each polygon:
>> v <- extract(SGRGBF40, calibf2,fun=mean,nl=3)
>> summary(v)
>      Length Class  Mode
>  [1,] 0      -none- NULL
>  [2,] 0      -none- NULL
>  [3,] 0      -none- NULL
>  [4,] 0      -none- NULL
>  [5,] 0      -none- NULL
> etc
>
>> sessionInfo()
> R version 2.12.2 (2011-02-25)
> Platform: i486-pc-linux-gnu (32-bit)
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] rgdal_0.6-31 raster_1.7-8 sp_0.9-72
>
> loaded via a namespace (and not attached):
> [1] grid_2.12.2     lattice_0.19-17 tools_2.12.2
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster: extract is empty with polygon

[Agustin Lobo]

Robert,

I have "reprojection on the fly" not activated. The raster
is an image from an standard camera, hence no projection.
The shapefile takes the projection from the project, but the
coordinates should be within the limits of the image, as it has been
digitized on top of the image. Thus I thought it was just a matter of
setting the CRS for the raster in R.

> intersectExtent(calibf2, SGRGBF40)
Error in intersectExtent(calibf2, SGRGBF40) : Invalid extent

which is in agreement with the rest of results, but still in
contradiction with what we see in QGIS

Agus

2011/3/24 Robert J. Hijmans :
> Hi Agus,
>
> It would be useful to see
>
> extent(calibf2)
>
> Apparently it does not overlap with SGRGBF40  (see
> intersectExtent(calibf2, SGRGBF40)  )
> because the polys did not plot on top of the raster.
>
> Perhaps QGIS does some on-the-fly projection trick or other such that
> the coordinates from your on-screen digitization do match those of the
> raster?
>
> This will not fix that! :
>> projection(SGRGBF40) = projection(calibf2)
>
> You can perhaps compare with a polygon you draw in R on top of SGRGBF40
>
> plot(SGRGBF40,1)
> p <- drawPoly()
> extract(SGRGBF40, p)
>
> And save p and load into QGIS?
>
> Best, Robert
>
> ( I would also update the raster package (but I do not think that
> matters for this problem).
>
>
> On Thu, Mar 24, 2011 at 11:27 AM, Agustin Lobo  wrote:
>> I forgot to add that raster and vector do overlap in qgis (the vector
>> was actually digitized on top of the raster)
>> Agus
>>
>> 2011/3/24 Agustin Lobo :
>>> Hi!
>>> I'm trying to calculate the mean values for a set of polygons.
>>> This is what I do:
>>>
>>> require(raster)
>>> require(rgdal)
>>> #SGRGBWBPS125F40
>>> #"read" tif file
>>> SGRGBF40 = 
>>> brick("/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
>>>> SGRGBF40
>>> class       : RasterBrick
>>> filename    : 
>>> /media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif
>>> nlayers     : 3
>>> nrow        : 1760
>>> ncol        : 2640
>>> ncell       : 4646400
>>> projection  : NA
>>> min value   : 0 0 0
>>> max value   : 65535 65535 65535
>>> extent      : 0, 2640, 0, 1760  (xmin, xmax, ymin, ymax)
>>> resolution  : 1, 1  (x, y)
>>>> projection(SGRGBF40)
>>> [1] "NA"
>>>
>>> #read the shape file
>>> calibf2 = 
>>> readOGR(dsn="/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
>>>> projection(calibf2)
>>> [1] "+proj=utm +zone=31 +ellps=intl +units=m +no_defs"
>>>
>>>> projection(SGRGBF40) = projection(calibf2)
>>>
>>> But they do not overlap as shown by
>>>> plot(subset(SGRGBF40,1))
>>>> plot(calibf2,add=T)
>>>
>>> and the extracted object is empty
>>> #Calculate mean values for each polygon:
>>>> v <- extract(SGRGBF40, calibf2,fun=mean,nl=3)
>>>> summary(v)
>>>      Length Class  Mode
>>>  [1,] 0      -none- NULL
>>>  [2,] 0      -none- NULL
>>>  [3,] 0      -none- NULL
>>>  [4,] 0      -none- NULL
>>>  [5,] 0      -none- NULL
>>> etc
>>>
>>>> sessionInfo()
>>> R version 2.12.2 (2011-02-25)
>>> Platform: i486-pc-linux-gnu (32-bit)
>>>
>>> locale:
>>>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>>>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>>>  [5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8
>>>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
>>>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
>>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>>
>>> attached base packages:
>>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>>
>>> other attached packages:
>>> [1] rgdal_0.6-31 raster_1.7-8 sp_0.9-72
>>>
>>> loaded via a namespace (and not attached):
>>> [1] grid_2.12.2     lattice_0.19-17 tools_2.12.2
>>>
>>
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Fwd: raster: extract is empty with polygon

[Agustin Lobo]

But then this is not as in gdal. Gdalinfo states that the bottom left is 0,1760
while R states 0,0
Should not raster do as gdal?

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1760.0)
Upper Right ( 2640.0,    0.0)
Lower Right ( 2640.0, 1760.0)
Center      ( 1320.0,  880.0)
Band 1 Block=2640x1 Type=UInt16, ColorInterp=Red
Band 2 Block=2640x1 Type=UInt16, ColorInterp=Green
Band 3 Block=2640x1 Type=UInt16, ColorInterp=Blue

Agus

2011/3/24 Robert J. Hijmans :
> Agus,
>
> If this is from a standard camera, perhaps the problem is that there
> is no associated georeference, and different programs may come up with
> a different solution to put them somewhere on the map.
> R (via gdal) does this: 0, 2640, 0, 1760 (xmin, xmax, ymin, ymax),
> with resolution 1 (cell). What does QGIS do?
> Perhaps you need to first save it to a geotiff, or add a (fake
> coordinates) "world file" so that the treatment becomes consistent.
>
> What is extent(calibf2) ?
>
> Robert
>
> On Thu, Mar 24, 2011 at 12:56 PM, Agustin Lobo  wrote:
>> Robert,
>>
>> I have "reprojection on the fly" not activated. The raster
>> is an image from an standard camera, hence no projection.
>> The shapefile takes the projection from the project, but the
>> coordinates should be within the limits of the image, as it has been
>> digitized on top of the image. Thus I thought it was just a matter of
>> setting the CRS for the raster in R.
>>
>>> intersectExtent(calibf2, SGRGBF40)
>> Error in intersectExtent(calibf2, SGRGBF40) : Invalid extent
>>
>> which is in agreement with the rest of results, but still in
>> contradiction with what we see in QGIS
>>
>> Agus
>>
>> 2011/3/24 Robert J. Hijmans :
>>> Hi Agus,
>>>
>>> It would be useful to see
>>>
>>> extent(calibf2)
>>>
>>> Apparently it does not overlap with SGRGBF40  (see
>>> intersectExtent(calibf2, SGRGBF40)  )
>>> because the polys did not plot on top of the raster.
>>>
>>> Perhaps QGIS does some on-the-fly projection trick or other such that
>>> the coordinates from your on-screen digitization do match those of the
>>> raster?
>>>
>>> This will not fix that! :
>>>> projection(SGRGBF40) = projection(calibf2)
>>>
>>> You can perhaps compare with a polygon you draw in R on top of SGRGBF40
>>>
>>> plot(SGRGBF40,1)
>>> p <- drawPoly()
>>> extract(SGRGBF40, p)
>>>
>>> And save p and load into QGIS?
>>>
>>> Best, Robert
>>>
>>> ( I would also update the raster package (but I do not think that
>>> matters for this problem).
>>>
>>>
>>> On Thu, Mar 24, 2011 at 11:27 AM, Agustin Lobo  
>>> wrote:
>>>> I forgot to add that raster and vector do overlap in qgis (the vector
>>>> was actually digitized on top of the raster)
>>>> Agus
>>>>
>>>> 2011/3/24 Agustin Lobo :
>>>>> Hi!
>>>>> I'm trying to calculate the mean values for a set of polygons.
>>>>> This is what I do:
>>>>>
>>>>> require(raster)
>>>>> require(rgdal)
>>>>> #SGRGBWBPS125F40
>>>>> #"read" tif file
>>>>> SGRGBF40 = 
>>>>> brick("/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
>>>>>> SGRGBF40
>>>>> class       : RasterBrick
>>>>> filename    : 
>>>>> /media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif
>>>>> nlayers     : 3
>>>>> nrow        : 1760
>>>>> ncol        : 2640
>>>>> ncell       : 4646400
>>>>> projection  : NA
>>>>> min value   : 0 0 0
>>>>> max value   : 65535 65535 65535
>>>>> extent      : 0, 2640, 0, 1760  (xmin, xmax, ymin, ymax)
>>>>> resolution  : 1, 1  (x, y)
>>>>>> projection(SGRGBF40)
>>>>> [1] "NA"
>>>>>
>>>>> #read the shape file
>>>>> calibf2 = 
>>>>> readOGR(dsn="/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
>>>>>> projection(calibf2)
>>>>> [1] "+proj=utm +zone=31 +ellps=intl +units=m +no_defs"
>>>>>
>>>>>> projection(SGRGBF40) = proj

Re: [R-sig-Geo] raster: extract is empty with polygon

[Agustin Lobo]

Anyway, the problem is solved by including the appropriate wld file:
1.0
0
0
-1.0
0.5
0.5

In this case, both QGIS and R (raster) keep the same coordinates and
the vector layer matches the raster layer.

Agus

2011/3/25 Agustin Lobo :
> But then this is not as in gdal. Gdalinfo states that the bottom left is 
> 0,1760
> while R states 0,0
> Should not raster do as gdal?
>
> Corner Coordinates:
> Upper Left  (    0.0,    0.0)
> Lower Left  (    0.0, 1760.0)
> Upper Right ( 2640.0,    0.0)
> Lower Right ( 2640.0, 1760.0)
> Center      ( 1320.0,  880.0)
> Band 1 Block=2640x1 Type=UInt16, ColorInterp=Red
> Band 2 Block=2640x1 Type=UInt16, ColorInterp=Green
> Band 3 Block=2640x1 Type=UInt16, ColorInterp=Blue
>
> Agus
>
> 2011/3/24 Robert J. Hijmans :
>> Agus,
>>
>> If this is from a standard camera, perhaps the problem is that there
>> is no associated georeference, and different programs may come up with
>> a different solution to put them somewhere on the map.
>> R (via gdal) does this: 0, 2640, 0, 1760 (xmin, xmax, ymin, ymax),
>> with resolution 1 (cell). What does QGIS do?
>> Perhaps you need to first save it to a geotiff, or add a (fake
>> coordinates) "world file" so that the treatment becomes consistent.
>>
>> What is extent(calibf2) ?
>>
>> Robert
>>
>> On Thu, Mar 24, 2011 at 12:56 PM, Agustin Lobo  wrote:
>>> Robert,
>>>
>>> I have "reprojection on the fly" not activated. The raster
>>> is an image from an standard camera, hence no projection.
>>> The shapefile takes the projection from the project, but the
>>> coordinates should be within the limits of the image, as it has been
>>> digitized on top of the image. Thus I thought it was just a matter of
>>> setting the CRS for the raster in R.
>>>
>>>> intersectExtent(calibf2, SGRGBF40)
>>> Error in intersectExtent(calibf2, SGRGBF40) : Invalid extent
>>>
>>> which is in agreement with the rest of results, but still in
>>> contradiction with what we see in QGIS
>>>
>>> Agus
>>>
>>> 2011/3/24 Robert J. Hijmans :
>>>> Hi Agus,
>>>>
>>>> It would be useful to see
>>>>
>>>> extent(calibf2)
>>>>
>>>> Apparently it does not overlap with SGRGBF40  (see
>>>> intersectExtent(calibf2, SGRGBF40)  )
>>>> because the polys did not plot on top of the raster.
>>>>
>>>> Perhaps QGIS does some on-the-fly projection trick or other such that
>>>> the coordinates from your on-screen digitization do match those of the
>>>> raster?
>>>>
>>>> This will not fix that! :
>>>>> projection(SGRGBF40) = projection(calibf2)
>>>>
>>>> You can perhaps compare with a polygon you draw in R on top of SGRGBF40
>>>>
>>>> plot(SGRGBF40,1)
>>>> p <- drawPoly()
>>>> extract(SGRGBF40, p)
>>>>
>>>> And save p and load into QGIS?
>>>>
>>>> Best, Robert
>>>>
>>>> ( I would also update the raster package (but I do not think that
>>>> matters for this problem).
>>>>
>>>>
>>>> On Thu, Mar 24, 2011 at 11:27 AM, Agustin Lobo  
>>>> wrote:
>>>>> I forgot to add that raster and vector do overlap in qgis (the vector
>>>>> was actually digitized on top of the raster)
>>>>> Agus
>>>>>
>>>>> 2011/3/24 Agustin Lobo :
>>>>>> Hi!
>>>>>> I'm trying to calculate the mean values for a set of polygons.
>>>>>> This is what I do:
>>>>>>
>>>>>> require(raster)
>>>>>> require(rgdal)
>>>>>> #SGRGBWBPS125F40
>>>>>> #"read" tif file
>>>>>> SGRGBF40 = 
>>>>>> brick("/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
>>>>>>> SGRGBF40
>>>>>> class       : RasterBrick
>>>>>> filename    : 
>>>>>> /media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif
>>>>>> nlayers     : 3
>>>>>> nrow        : 1760
>>>>>> ncol        : 2640
>>>>>> ncell       : 4646400
>>>>>> projection  : NA
>>>>>> min value   : 0 0 0
>>>>>> max value   : 65535 65535

[R-sig-Geo] raster::extract NAs introduced?

[Agustin Lobo]

Hi

I've found that extract() introduces NAs:
SGRGBF40 = 
brick("/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
calibf2 = 
readOGR(dsn="/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
projection(SGRGBF40) = projection(calibf2)
plot(subset(SGRGBF40,1))
plot(calibf2,add=T)
(overlay is ok now)

> summary(SGRGBF40)
Cells:  4646400
NAs  :  0 0 0

1 2 3
Min. 2629 3 3
1st Qu. 34100 22200 18290
Median  42470 27480 24680
Mean40900 28510 25360
3rd Qu. 48360 33100 29010
Max.65510 65510 64300
summary based on a sample of 5000 cells, which is 0.107610192837466 %
of all cells
To make sure (NAs could outside the sample):
> any(is.na(SGRGBF40))
class   : RasterLayer
dimensions  : 1760, 2640, 1  (nrow, ncol, nlayers)
resolution  : 1, 1  (x, y)
extent  : 0, 2640, -1759, 1  (xmin, xmax, ymin, ymax)
projection  : +proj=utm +zone=31 +ellps=intl +units=m +no_defs
values  : in memory
min value   : 0
max value   : 1
So no NAs
But:
> summary(v)
 SGRGBWBPS125F40_1 SGRGBWBPS125F40_2 SGRGBWBPS125F40_3
 Min.   :18444 Min.   :  119.4   Min.   :3
 1st Qu.:33227 1st Qu.:25516.4   1st Qu.:15633
 Median :37752 Median :36264.1   Median :22374
 Mean   :40208 Mean   :36275.1   Mean   :29737
 3rd Qu.:46899 3rd Qu.:48052.5   3rd Qu.:48116
 Max.   :64992 Max.   :63667.1   Max.   :63826
 NA's   :5 NA's   :1.0


Any explanation to this? Not a big problem (can use na.rm) but I'm
concerned on this problem implying
the means not being correct.

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Problem with raster::extract

[Agustin Lobo]

Hi!

I'm doing

SGRGBF40 = 
brick("/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
calibf2 = 
readOGR(dsn="/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
projection(SGRGBF40) = projection(calibf2)
#I make sure vector calibf2 and raster SGRGBF40 do overlap:
plot(subset(SGRGBF40,1))
plot(calibf2,add=T)

#Calculate mean values for each polygon:
v <- extract(SGRGBF40, calibf2,fun=mean,na.rm=T,nl=3)
summary(v)
SGRGBF40data = cbind(calibf2@data,v)

But then, despite na.rm=T:
> SGRGBF40data[7,]
  ID code codetxt SGRGBWBPS125F40_1 SGRGBWBPS125F40_2 SGRGBWBPS125F40_3
6  77  C7   NaN  39781.02 3

> summary(v)
 SGRGBWBPS125F40_1 SGRGBWBPS125F40_2 SGRGBWBPS125F40_3
 Min.   :18444 Min.   :  119.4   Min.   :3
 1st Qu.:34098 1st Qu.:25875.7   1st Qu.:15633
 Median :37776 Median :36668.2   Median :22374
 Mean   :41314 Mean   :36882.2   Mean   :29737
 3rd Qu.:52236 3rd Qu.:49452.6   3rd Qu.:48116
 Max.   :65396 Max.   :65417.5   Max.   :63826
 NA's   :3

While there are no NA in the raster brick:
> summary(SGRGBF40)
Cells:  4646400
NAs  :  0 0 0

1 2 3
Min. 1777 3 3
1st Qu. 34910 22460 18640
Median  42660 27540 24420
Mean41070 28570 25170
3rd Qu. 48470 33320 29060
Max.65520 65530 64220


A possibility would be a division by 0 (weird, because all polygons
are large), so I've tried:
> vl <- extract(SGRGBF40, calibf2,fun=length,nl=3)
but...
Error in t(sapply(res[!i], function(x) apply(x, 2, FUN = fun, na.rm =
na.rm))) :
  error in evaluating the argument 'x' in selecting a method for
function 't': Error in FUN(newX[, i], ...) :
  2 arguments passed to 'length' which requires 1
Calls: sapply -> lapply -> FUN -> apply
Calls: extract -> extract -> .polygonValues -> t

Any hint?

Data:
http://dl.dropbox.com/u/3180464/SGRGBWBPS125F40.tif
http://dl.dropbox.com/u/3180464/SGRGBWBPS125F40.tfw
http://dl.dropbox.com/u/3180464/calibf2.zip

> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8   LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8   LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8  LC_NAME=en_US.UTF-8
 [9] LC_ADDRESS=en_US.UTF-8LC_TELEPHONE=en_US.UTF-8
[11] LC_MEASUREMENT=en_US.UTF-8LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] grid  stats graphics  grDevices utils datasets  methods
[8] base

other attached packages:
[1] gridExtra_0.7 ggplot2_0.8.9 rgdal_0.6-33  raster_1.7-29 sp_0.9-76
[6] reshape_0.8.4 plyr_1.4  proto_0.3-8   rkward_0.5.6

loaded via a namespace (and not attached):
[1] lattice_0.19-26 tools_2.13.0
>

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Problem with raster::calc

[Agustin Lobo]

I've calculated a linear model using a subset of data from raster files
> sgcirRlm = lm(SGRGBF40data[,4] ~ 
> SGCIRF56data[,4]+SGCIRF56data[,5]+SGCIRF56data[,6]
> coefficients(sgcirRlm)
  (Intercept) SGCIRF56data[, 4] SGCIRF56data[, 5] SGCIRF56data[, 6]
-2075.0432230 0.9229117 0.9963710-0.8310794

and want to apply the coefficients to the multiband raster
> show(SGCIRF56)
class   : RasterBrick
dimensions  : 1760, 2640, 3  (nrow, ncol, nlayers)
resolution  : 1, 1  (x, y)
extent  : 0, 2640, -1759, 1  (xmin, xmax, ymin, ymax)
projection  : +proj=utm +zone=31 +ellps=intl +units=m +no_defs
values  : 
/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL/SGCIR/SGCIRWBPS125F56v2.tif
min values  : 0 0 0
max values  : 65535 65535 65535

to calculate the predicted raster, for which I'm using subset() within calc():
> sgcirRpred = 
> calc(SGCIRF56,fun=function(x){coefficients(sgcirRlm)[1]+coefficients(sgcirRlm)[2]*subset(x,1)
>  + coefficients(sgcirRlm)[3]*subset(x,2) + 
> coefficients(sgcirRlm)[4]*subset(x,3)})

but this results into an error:
Error in .local(x, fun, ...) : cannot use this function
Calls: calc -> calc -> .local

is this because subset() cannot be used within the function to be
provided to calc()?

Thanks
> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8   LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8   LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8  LC_NAME=en_US.UTF-8
 [9] LC_ADDRESS=en_US.UTF-8LC_TELEPHONE=en_US.UTF-8
[11] LC_MEASUREMENT=en_US.UTF-8LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] grid  stats graphics  grDevices utils datasets  methods
[8] base

other attached packages:
[1] gridExtra_0.7 ggplot2_0.8.9 rgdal_0.6-33  raster_1.8-16 sp_0.9-76
[6] reshape_0.8.4 plyr_1.4  proto_0.3-8   rkward_0.5.6

loaded via a namespace (and not attached):
[1] lattice_0.19-26 tools_2.13.0

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Problem with raster::extract

[Agustin Lobo]

I confirm the problem is solved in  raster_1.8-16

Thanks!

Agus

2011/4/28 Robert J. Hijmans :
> This was caused by raster incorrectly interpreting the highest
> possible INT2U value (65530) as NA. Now fixed (version 1.8-16), at
> least for this case. Thanks for reporting, Robert
>
> On Wed, Apr 27, 2011 at 8:52 AM, Robert J. Hijmans  
> wrote:
>> Hello Agus,
>>
>> The NaN suggests that there are only NA values for that polygon/layer.
>> To investigate this I would do
>>
>> v <- extract(SGRGBF40, calibf2, nl=3)
>>
>> # i.e. do not use a function and then inspect
>>
>> v[[7]]
>>
>>
>> All functions must accept an na.rm argument (even if they choose to
>> ignore it. So for length you can do something like:
>>
>> vl <- extract(SGRGBF40, calibf2, fun=function(x,...)length(x), nl=3)
>>
>>
>> Hope this helps, Robert
>>
>>> v
>> [1] 387.8158 329.3913
>>> v <- extract(r, polys, length)
>> Error in FUN(X[[1L]], ...) :
>>  2 arguments passed to 'length' which requires 1
>>> v <- extract(r, polys, function(x,...)length(x))
>>> v
>> [1] 38 23
>>>
>>
>>
>> On Wed, Apr 27, 2011 at 3:14 AM, Agustin Lobo  wrote:
>>> Hi!
>>>
>>> I'm doing
>>>
>>> SGRGBF40 = 
>>> brick("/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
>>> calibf2 = 
>>> readOGR(dsn="/media/Iomega_HDD/UAVetal/CALIBRACIONRADIOM/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
>>> projection(SGRGBF40) = projection(calibf2)
>>> #I make sure vector calibf2 and raster SGRGBF40 do overlap:
>>> plot(subset(SGRGBF40,1))
>>> plot(calibf2,add=T)
>>>
>>> #Calculate mean values for each polygon:
>>> v <- extract(SGRGBF40, calibf2,fun=mean,na.rm=T,nl=3)
>>> summary(v)
>>> SGRGBF40data = cbind(calibf2@data,v)
>>>
>>> But then, despite na.rm=T:
>>>> SGRGBF40data[7,]
>>>  ID code codetxt SGRGBWBPS125F40_1 SGRGBWBPS125F40_2 SGRGBWBPS125F40_3
>>> 6  7    7      C7               NaN          39781.02                 3
>>>
>>>> summary(v)
>>>  SGRGBWBPS125F40_1 SGRGBWBPS125F40_2 SGRGBWBPS125F40_3
>>>  Min.   :18444     Min.   :  119.4   Min.   :    3
>>>  1st Qu.:34098     1st Qu.:25875.7   1st Qu.:15633
>>>  Median :37776     Median :36668.2   Median :22374
>>>  Mean   :41314     Mean   :36882.2   Mean   :29737
>>>  3rd Qu.:52236     3rd Qu.:49452.6   3rd Qu.:48116
>>>  Max.   :65396     Max.   :65417.5   Max.   :63826
>>>  NA's   :    3
>>>
>>> While there are no NA in the raster brick:
>>>> summary(SGRGBF40)
>>> Cells:  4646400
>>> NAs  :  0 0 0
>>>
>>>            1     2     3
>>> Min.     1777     3     3
>>> 1st Qu. 34910 22460 18640
>>> Median  42660 27540 24420
>>> Mean    41070 28570 25170
>>> 3rd Qu. 48470 33320 29060
>>> Max.    65520 65530 64220
>>>
>>>
>>> A possibility would be a division by 0 (weird, because all polygons
>>> are large), so I've tried:
>>>> vl <- extract(SGRGBF40, calibf2,fun=length,nl=3)
>>> but...
>>> Error in t(sapply(res[!i], function(x) apply(x, 2, FUN = fun, na.rm =
>>> na.rm))) :
>>>  error in evaluating the argument 'x' in selecting a method for
>>> function 't': Error in FUN(newX[, i], ...) :
>>>  2 arguments passed to 'length' which requires 1
>>> Calls: sapply -> lapply -> FUN -> apply
>>> Calls: extract -> extract -> .polygonValues -> t
>>>
>>> Any hint?
>>>
>>> Data:
>>> http://dl.dropbox.com/u/3180464/SGRGBWBPS125F40.tif
>>> http://dl.dropbox.com/u/3180464/SGRGBWBPS125F40.tfw
>>> http://dl.dropbox.com/u/3180464/calibf2.zip
>>>
>>>> sessionInfo()
>>> R version 2.13.0 (2011-04-13)
>>> Platform: x86_64-pc-linux-gnu (64-bit)
>>>
>>> locale:
>>>  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C
>>>  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8
>>>  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8
>>>  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8
>>>  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8
>>> [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
>>>
>>> attached base packages:
>>> [1] grid      stats     graphics  grDevices utils     datasets  methods
>>> [8] base
>>>
>>> other attached packages:
>>> [1] gridExtra_0.7 ggplot2_0.8.9 rgdal_0.6-33  raster_1.7-29 sp_0.9-76
>>> [6] reshape_0.8.4 plyr_1.4      proto_0.3-8   rkward_0.5.6
>>>
>>> loaded via a namespace (and not attached):
>>> [1] lattice_0.19-26 tools_2.13.0
>>>>
>>>
>>> Agus
>>>
>>
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster or rgdal: problem with geometry of tif + tifw file

[Agustin Lobo]

Hi!

I read a tif fle (with a companion tifw file) for which gdalinfo states:

alobo@alobo-laptop:/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2$
gdalinfo SDIM0114.tif
Warning 1: TIFFReadDirectory:Unknown field with tag 37717 (0x9355) encountered
Driver: GTiff/GeoTIFF
Files: SDIM0114.tif
Size is 2640, 1760
Coordinate System is `'
GeoTransform =
  445548.7745816645, -0.109570548549058, 0.216912817625298
  4628418.279379116, 0.216912817625298, 0.109570548549058
Metadata:
  TIFFTAG_SOFTWARE=SIGMA Photo Pro 3.5.2.
  TIFFTAG_DATETIME=2010:02:27 15:40:55
  TIFFTAG_XRESOLUTION=180
  TIFFTAG_YRESOLUTION=180
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  445548.775, 4628418.279)
Lower Left  (  445930.541, 4628611.124)
Upper Right (  445259.508, 4628990.929)
Lower Right (  445641.275, 4629183.773)
Center  (  445595.025, 4628801.026)


Then I read it in R:
> test= brick("/media/TRANSCEND2/VegcamPro/MATA20090729/TIFFOri/SDIM0114.tif")
> show(test)
class   : RasterBrick
dimensions  : 1760, 2640, 3  (nrow, ncol, nlayers)
resolution  : 0.1095705, 0.1095705  (x, y)
extent  : 445548.8, 445838, 4628418, 4628611  (xmin, xmax, ymin, ymax)
projection  : NA
values  : /media/TRANSCEND2/VegcamPro/MATA20090729/TIFFOri/SDIM0114.tif
min values  : 0 0 0
max values  : 65535 65535 65535

and write it back to disk:

> writeRaster(test,file="test.tif", format="GTiff", 
> overwrite=TRUE,datatype='FLT8S')

But gdalinfo states a different bounding box for the new file:
alobo@alobo-laptop:/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2$
gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 2640, 1760
Coordinate System is `'
Origin = (445548.774581664009020,4628611.233115110546350)
Pixel Size = (0.109570548549241,-0.109570548548469)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  445548.775, 4628611.233)
Lower Left  (  445548.775, 4628418.389)
Upper Right (  445838.041, 4628611.233)
Lower Right (  445838.041, 4628418.389)
Center  (  445693.408, 4628514.811)
.../...

Somehow the geometric properties are modified in R

I've made a single geotif file with gdal_transform (no tifw created)
and read it in R in the same way and get the same result, so the
problem does not come
from using a tifw file as I initially thought.

Perhaps a bug in brick() ?

Thanks

Agus

Files:
http://dl.dropbox.com/u/3180464/SDIM0114.tif
http://dl.dropbox.com/u/3180464/SDIM0114.tifw

> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: i486-pc-linux-gnu (32-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C
LC_TIME=en_US.UTF-8
 [4] LC_COLLATE=en_US.UTF-8LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8  LC_NAME=en_US.UTF-8
LC_ADDRESS=en_US.UTF-8
[10] LC_TELEPHONE=en_US.UTF-8  LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] grid  stats graphics  grDevices utils datasets
methods   base

other attached packages:
[1] gridExtra_0.7 ggplot2_0.8.8 rgdal_0.6-31  raster_1.8-15 sp_0.9-72
   reshape_0.8.3
[7] plyr_1.2.1proto_0.3-8   rkward_0.5.6

loaded via a namespace (and not attached):
[1] lattice_0.19-26 tools_2.13.0

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster or rgdal: problem with geometry of tif + tifw file

[Agustin Lobo]

Robert,

I think a warning is better than an error. There are many operations
that can be performed
ignoring the rotation, it should be up to the user.

In my case the problem is actually created at writing (writeRaster()),
> writeRaster(round(testcirpred),file="/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2/delme.tif",
>  format="GTiff", overwrite=TRUE,datatype='INT2U')

because the new tif file gets wrong  geometric information that has
priority over the one provided by the tfw file.

I circumvent the problem by using
listgeo -no_norm SDIM0114.tif > original.geo
in a directory where the tfw file is not present. Then I apply this
geometric information to the output from R:
geotifcp -g original.geo delme.tif SDIM0114L1.tif
cp SDIM0114.tfw SDIM0114L1.tfw

more or less as suggested in
http://remotesensing.org/geotiff/faq.html

It would be good being able to just omit the geometric info for the
output tif file in writeRaster(), so the original tfw file would be
directly valid for the output file
to be displayed in QGIS, thus avoiding the need of the fix with
listgeo and geotifcp

Agus

2011/5/8 Robert J. Hijmans :
> Hi Agus,
>
> The problem is that you raster is rotated (GeoTransform[3] and [5] are not 
> zero)
>
> GeoTransform =
>  445548.7745816645, -0.109570548549058, 0.216912817625298
>  4628418.279379116, 0.216912817625298, 0.109570548549058
>
> raster ignored the rotation. I have fixed that (in version 1-8-19) by
> throwing an error for such a file.
>
> readGDAL will read the file, but as points (SpatialPointsDataFrame),
> which you could interpolate to a non-rotated raster. I do not plan to
> support such files in raster, but I might include a function for
> memory-safe transformation to a non-rotated file.
>
> Thanks for catching & reporting this,
>
> Robert
>
>
> On Sun, May 8, 2011 at 9:50 AM, Agustin Lobo  wrote:
>> Hi!
>>
>> I read a tif fle (with a companion tifw file) for which gdalinfo states:
>>
>> alobo@alobo-laptop:/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2$
>> gdalinfo SDIM0114.tif
>> Warning 1: TIFFReadDirectory:Unknown field with tag 37717 (0x9355) 
>> encountered
>> Driver: GTiff/GeoTIFF
>> Files: SDIM0114.tif
>> Size is 2640, 1760
>> Coordinate System is `'
>> GeoTransform =
>>  445548.7745816645, -0.109570548549058, 0.216912817625298
>>  4628418.279379116, 0.216912817625298, 0.109570548549058
>> Metadata:
>>  TIFFTAG_SOFTWARE=SIGMA Photo Pro 3.5.2.
>>  TIFFTAG_DATETIME=2010:02:27 15:40:55
>>  TIFFTAG_XRESOLUTION=180
>>  TIFFTAG_YRESOLUTION=180
>>  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
>> Image Structure Metadata:
>>  INTERLEAVE=PIXEL
>> Corner Coordinates:
>> Upper Left  (  445548.775, 4628418.279)
>> Lower Left  (  445930.541, 4628611.124)
>> Upper Right (  445259.508, 4628990.929)
>> Lower Right (  445641.275, 4629183.773)
>> Center      (  445595.025, 4628801.026)
>>
>>
>> Then I read it in R:
>>> test= brick("/media/TRANSCEND2/VegcamPro/MATA20090729/TIFFOri/SDIM0114.tif")
>>> show(test)
>> class       : RasterBrick
>> dimensions  : 1760, 2640, 3  (nrow, ncol, nlayers)
>> resolution  : 0.1095705, 0.1095705  (x, y)
>> extent      : 445548.8, 445838, 4628418, 4628611  (xmin, xmax, ymin, ymax)
>> projection  : NA
>> values      : /media/TRANSCEND2/VegcamPro/MATA20090729/TIFFOri/SDIM0114.tif
>> min values  : 0 0 0
>> max values  : 65535 65535 65535
>>
>> and write it back to disk:
>>
>>> writeRaster(test,file="test.tif", format="GTiff", 
>>> overwrite=TRUE,datatype='FLT8S')
>>
>> But gdalinfo states a different bounding box for the new file:
>> alobo@alobo-laptop:/media/TRANSCEND2/VegcamPro/MATA20090729/Proc2$
>> gdalinfo test.tif
>> Driver: GTiff/GeoTIFF
>> Files: test.tif
>> Size is 2640, 1760
>> Coordinate System is `'
>> Origin = (445548.774581664009020,4628611.233115110546350)
>> Pixel Size = (0.109570548549241,-0.109570548548469)
>> Image Structure Metadata:
>>  COMPRESSION=LZW
>>  INTERLEAVE=PIXEL
>> Corner Coordinates:
>> Upper Left  (  445548.775, 4628611.233)
>> Lower Left  (  445548.775, 4628418.389)
>> Upper Right (  445838.041, 4628611.233)
>> Lower Right (  445838.041, 4628418.389)
>> Center      (  445693.408, 4628514.811)
>> .../...
>>
>> Somehow the geometric properties are modified in R
>>
>> I've made a single geotif file with gdal_transform (no tifw created)
>> and read it in R in the same way and get the same result, so the
>> problem d

[R-sig-Geo] raster: creating a layer of NA

[Agustin Lobo]

HI!

I must create a layer with all missing values. I try:
> show(D10N2009)
class   : RasterBrick
dimensions  : 393, 606, 35  (nrow, ncol, nlayers)
resolution  : 0.008928571, 0.008928571  (x, y)
extent  : -2.004464, 3.40625, 40.49554, 44.00446  (xmin, xmax, ymin, ymax)
projection  : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
values  : /media/Iomega_HDD/FLUXPYR/VGTFLUXPYR1/D10/FLXP_NDD10_2009.vrt
min values  : 0 0 0 0 0 0 0 0 0 0 ...
max values  : 255 255 255 255 255 255 255 255 255 255 ...

> delme = D10N2009[[1]]*0
> show(delme)
class   : RasterLayer
dimensions  : 393, 606, 238158  (nrow, ncol, ncell)
resolution  : 0.008928571, 0.008928571  (x, y)
extent  : -2.004464, 3.40625, 40.49554, 44.00446  (xmin, xmax, ymin, ymax)
projection  : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
values  : in memory
min value   : 0
max value   : 0

But
> delme = D10N2009[[1]]*NA
Error in D10N2009[[1]] * NA : non-numeric argument to binary operator

I also try
> NAvalue(delme) <- 0
but no trace of the NAs:
> show(delme)
class   : RasterLayer
dimensions  : 393, 606, 238158  (nrow, ncol, ncell)
resolution  : 0.008928571, 0.008928571  (x, y)
extent  : -2.004464, 3.40625, 40.49554, 44.00446  (xmin, xmax, ymin, ymax)
projection  : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
values  : in memory
min value   : 0
max value   : 0

By now I write with NAflag=0 and read back, but there must be a
cleaner way of generating a layer with all NA values
(it's a missing layer for a particular time)

Thanks!

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Simulated rasters?

[Agustin Lobo]

Hi!
Could raster inherit from grf objects (package feoR)?
This would be useful

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] simulated rasters?

[Agustin Lobo]

Hi!

Could raster inherit from grf objects (package geoR)?
It would be useful to generate simulated raster fields.Or perhpas is there
a way of doing this with another package?
> sim <- grf(441, grid="reg", cov.pars=c(1, .25), nsim=4)
> simr = raster(sim)
Error in function (classes, fdef, mtable)  :
  unable to find an inherited method for function "raster", for signature "grf"

By now I use
> a = (sim$data)
> dim(a) <- c(21,21)
> simr = raster(t(a)[21:1,])

Thanks,

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster::plot() common color scale?

[Agustin Lobo]

Hi!
I'm doing

par(mfrow=c(2,2))
plot(simr)
plot(simrlisa)
a2 = simr>l$lisamax
plot(l$lisamax)
plot(simr*a2)

where simr, l%lisamax and a2 are raster objects. Is there any way to
get a common color scheme
(i.e., setting common min and max values)?
The color table is the same, but the scaling is independent for each
raster object and cannot be compared.

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] rastrer::raster no extent argument?

[Agustin Lobo]

Hi!

In order to define an arbitrary test raster object, I thought I could do:
> foo4 <- raster(nrows=20, ncols=20,crs=NA,extent=c(0,20,0,20))
> values(foo4) <- runif(400,-100,150)

but extent is not a valid argument, so I use another step
> foo4 <- raster(nrows=20, ncols=20,crs=NA)
> values(foo4) <- runif(400,-100,150)
> extent(foo4) = extent(c(0,20,0,20))

Would not make sense adding this feature to raster() itself?

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster::plot() common color scale?

[Agustin Lobo]

Paul,

While this is certainly the most elegant solution, would it work with
really large raster objects? Note that you
make a data.frame, and the most important characteristic of package
raster is its ability to do not incorporate
the actual raster values to the raster objects to save memory.

Would this problem be avoided by the new wrapper provided by Robert?

Thanks

Agus

2011/6/10 Paul Hiemstra :
>  Hi,
>
> I would recommend using one of the more advanced plotting facilities in
> R, ggplot. An example:
>
> library(ggplot2)
> library(raster)
> theme_set(theme_bw())
>
> r = raster(system.file("external/test.grd", package="raster"))
>
> # Convert to data.frame
> r_df = as.data.frame(as(r, 'SpatialPixelsDataFrame'))
> # create new column 10 times larger
> # here you could extract data from other grids that
> # you want to show at the same time
> r_df$valuesx10 = r_df$values*2
>
> # Reshape the data for ggplot
> plotData = melt(r_df, id.vars = c('x','y'))
>
> ggplot(aes(x = x, y = y), data = plotData) +
>    geom_tile(aes(fill = value)) + facet_wrap(~ variable) +
>    scale_fill_gradient(low = 'white', high = 'blue') +
>    coord_equal()
>
> This is a basic ggplot example of visualizing rasters. I realize that it
> is a short example with a lot of things specific to ggplot, but I hope
> you can figure out why I use the code that I use. In my view, the
> investment in learning to use ggplot is worth it. A good place to start
> is the ggplot website [1].
>
> cheers.
> Paul
>
> [1] http://had.co.nz/ggplot2/
>
> On 06/09/2011 01:59 PM, Agustin Lobo wrote:
>> Hi!
>> I'm doing
>>
>> par(mfrow=c(2,2))
>> plot(simr)
>> plot(simrlisa)
>> a2 = simr>l$lisamax
>> plot(l$lisamax)
>> plot(simr*a2)
>>
>> where simr, l%lisamax and a2 are raster objects. Is there any way to
>> get a common color scheme
>> (i.e., setting common min and max values)?
>> The color table is the same, but the scaling is independent for each
>> raster object and cannot be compared.
>>
>> Thanks
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
> --
> Paul Hiemstra, Ph.D.
> Global Climate Division
> Royal Netherlands Meteorological Institute (KNMI)
> Wilhelminalaan 10 | 3732 GK | De Bilt | Kamer B 3.39
> P.O. Box 201 | 3730 AE | De Bilt
> tel: +31 30 2206 494
>
> http://intamap.geo.uu.nl/~paul
> http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770
>
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster::plot() common color scale?

[Agustin Lobo]

Robert,

I've tried with raster_1.8-35
but:

other attached packages:
[1] raster_1.8-35 ggplot2_0.8.9 proto_0.3-8   reshape_0.8.4 plyr_1.4
[6] plotrix_3.2-2 sp_0.9-76 rkward_0.5.6

loaded via a namespace (and not attached):
[1] digest_0.4.2lattice_0.19-26 tools_2.13.0


> foo1 <- foo2 <- foo3<- foo4 <- raster(nrows=20, ncols=20)
> values(foo1) <- runif(400,-100.100)
> values(foo2) <- runif(400,-150,100)
> values(foo3) <- runif(400,0,50)
> values(foo4) <- runif(400,-100,150)
> extent(foo1) = extent(foo2) = extent(foo3) = extent(foo4) = 
> extent(c(0,20,0,20))
> s <- stack(foo1,foo2,foo3,foo4)

> spplot(as(x, 'SpatialGridDataFrame'))
Error in rk.record.plot$.save.tlo.in.hP() :
  could not find function "trellis.last.object"
Calls: print ... print -> print.trellis -> printFunction -> 
> spplot(s)
Error in function (classes, fdef, mtable)  :
  unable to find an inherited method for function "spplot", for
signature "RasterStack"
Calls: spplot -> 
> theme_set(theme_bw())
+ gplot(s) + geom_tile(aes(fill = value)) + facet_wrap(~ variable) +
+scale_fill_gradient(low = 'white', high = 'blue') +
+ coord_equal()
Error: could not find function "gplot"
> theme_set(theme_bw())
+ ggplot(s) + geom_tile(aes(fill = value)) + facet_wrap(~ variable) +
+scale_fill_gradient(low = 'white', high = 'blue') +
+ coord_equal()
Error: ggplot2 doesn't know how to deal with data of class RasterStack
>


Agus

2011/6/11 Robert Hijmans :
>> Is there any way to
>> get a common color scheme
>> (i.e., setting common min and max values)?
>
> Here is yet another way:
>
> library(raster)
> foo <- raster(nrows=20, ncols=20)
> values(foo) <- runif(400)
> bar <- raster(nrows=20, ncols=20)
> values(bar) <- runif(400)
> s <- stack(foo, bar, foo*2)
>
> # first sample large rasters:
> x <- sampleRegular(s, 25000, asRaster=T)
>
> # then use spplot
> spplot(as(x, 'SpatialGridDataFrame'))
>
>
> # I have added two functions to raster version 1.8-33
> # a generic function for spplot and Raster objects
> # that does the same as the example above in one short line:
>
> spplot(s)
>
> # and a wrapper around ggplot (called gplot),
> # based on Paul's example, that allows you to do things like:
>
> theme_set(theme_bw())
>
> gplot(s) + geom_tile(aes(fill = value)) + facet_wrap(~ variable) +
>            scale_fill_gradient(low = 'white', high = 'blue') +
> coord_equal()
>
>
> Robert
>
>
> --
> View this message in context: 
> http://r-sig-geo.2731867.n2.nabble.com/raster-plot-common-color-scale-tp6457926p6464547.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster: cannot extract object from list

[Agustin Lobo]

A function of mine outputs several raster objects in a list "l":
> class(l$lisaobsv)
[1] "RasterLayer"
attr(,"package")
[1] "raster"

But by some reason I get this weird error:
> a = l$lisaobsv
Error in slot(x, slotname) :
  no slot of name "z" for this object of class "RasterLayer"
Calls:  -> lapply -> FUN -> slot -> .Call
Warning message:
In .rk.get.structure.global("a") : failure to get object a

Is this a bug or is there an special way to extract raster objects from a list?

Thanks

Agus

> sessionInfo()
R version 2.13.0 (2011-04-13)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8   LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8   LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8  LC_NAME=en_US.UTF-8
 [9] LC_ADDRESS=en_US.UTF-8LC_TELEPHONE=en_US.UTF-8
[11] LC_MEASUREMENT=en_US.UTF-8LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] grid  stats graphics  grDevices utils datasets  methods
[8] base

other attached packages:
 [1] rgdal_0.6-33  geoR_1.6-27   plotrix_3.2-2 ggplot2_0.8.9 proto_0.3-8
 [6] reshape_0.8.4 plyr_1.4  raster_1.8-40 sp_0.9-76 rkward_0.5.6

loaded via a namespace (and not attached):
[1] lattice_0.19-26 tools_2.13.0

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] raster: cannot extract object from list

[Agustin Lobo]

Solved: it was caused by having upgraded to raster 1.8.40 but using
raster objects created under 1.8.38

Agus

2011/7/6 Agustin Lobo :
> A function of mine outputs several raster objects in a list "l":
>> class(l$lisaobsv)
> [1] "RasterLayer"
> attr(,"package")
> [1] "raster"
>
> But by some reason I get this weird error:
>> a = l$lisaobsv
> Error in slot(x, slotname) :
>  no slot of name "z" for this object of class "RasterLayer"
> Calls:  -> lapply -> FUN -> slot -> .Call
> Warning message:
> In .rk.get.structure.global("a") : failure to get object a
>
> Is this a bug or is there an special way to extract raster objects from a 
> list?
>
> Thanks
>
> Agus
>
>> sessionInfo()
> R version 2.13.0 (2011-04-13)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C
>  [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8
>  [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8
> [11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8
>
> attached base packages:
> [1] grid      stats     graphics  grDevices utils     datasets  methods
> [8] base
>
> other attached packages:
>  [1] rgdal_0.6-33  geoR_1.6-27   plotrix_3.2-2 ggplot2_0.8.9 proto_0.3-8
>  [6] reshape_0.8.4 plyr_1.4      raster_1.8-40 sp_0.9-76     rkward_0.5.6
>
> loaded via a namespace (and not attached):
> [1] lattice_0.19-26 tools_2.13.0
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Summary: how to plot several raster layers on the same page and with the same color scale

[Agustin Lobo]

I summarize here the different solutions I got to my question regarding
how to plot several raster layers on the same page and with the same
color scale.
Sorry it took a long time, I was interrupted by something else for weeks.
Hope this is useful.
Many thanks to all.


require(raster)
require(ggplot2)
require(plotrix)
require(rasterVis)
require(lattice)

foo1 <- foo2 <- foo3<- foo4 <- raster(nrows=20, ncols=20)
values(foo1) <- runif(400,-100.100)
values(foo2) <- runif(400,-150,100)
values(foo3) <- runif(400,0,50)
values(foo4) <- runif(400,-100,150)
extent(foo1) = extent(foo2) = extent(foo3) = extent(foo4) = extent(c(0,20,0,20))


#Based on comment by Jeremy Raw and andy.b...@wwu.edu
par(mfrow=c(2,2))
plot(foo1,col=rev(heat.colors(32)),zlim=c(-150,150))
plot(foo2,col=rev(heat.colors(32)),zlim=c(-150,150))
plot(foo3,col=rev(heat.colors(32)),zlim=c(-150,150))
plot(foo4,col=rev(heat.colors(32)),zlim=c(-150,150))

#Note by ALOBO: simplest solution

#Based on comment by Adam Sparks
require(plotrix)
minValue = -150
maxValue = 150
colstep=10
brk = c(minValue, maxValue, by = colstep)
brk = seq(minValue, maxValue, by = colstep)
par(mfrow=c(2,2))
plot(foo1, breaks = brk, col = rev(heat.colors(30)), axes = FALSE,
xlab = '', ylab = '', legend = FALSE)
color.legend(xl=20.5,yb=1,xr=21.5,yt=18, legend = seq(minValue,
maxValue, by = 3*colstep),
  rect.col = rev(heat.colors(30)),pos=4,gradient="y",cex=0.55,offset=2,
col="black")
plot(foo2, breaks = brk, col = rev(heat.colors(30)),legend = FALSE,
xlab = '', ylab = '',)
color.legend(xl=20.5,yb=1,xr=21.5,yt=18, legend = seq(minValue,
maxValue, by = 3*colstep),
  rect.col = rev(heat.colors(30)),pos=4,gradient="y",cex=0.55,offset=2,
col="black")
plot(foo3, breaks = brk, col = rev(heat.colors(30)),legend = FALSE,
xlab = '', ylab = '',)
color.legend(xl=20.5,yb=1,xr=21.5,yt=18, legend = seq(minValue,
maxValue, by = 3*colstep),
  rect.col = rev(heat.colors(30)),pos=4,gradient="y",cex=0.55,offset=2,
col="black")
plot(foo3, breaks = brk, col = rev(heat.colors(30)),legend = FALSE,
xlab = '', ylab = '',)
color.legend(xl=20.5,yb=1,xr=21.5,yt=18, legend = seq(minValue,
maxValue, by = 3*colstep),
  rect.col = rev(heat.colors(30)),pos=4,gradient="y",cex=0.55,offset=2,
col="black")


#Note by ALOBO: in practice, you must use hist() to determine
meaningful min and max values in both cases


#Based on comment by  paul.hiems...@knmi.nl
library(ggplot2)
library(raster)
theme_set(theme_bw())

r = foo1

# Convert to data.frame
r_df = as.data.frame(as(r, 'SpatialPixelsDataFrame'))
print(str(r_df))
# create new column 10 times larger
# here you could extract data from other grids that
# you want to show at the same time
r_df$foo2 = values(foo2)
r_df$foo3 = values(foo3)
r_df$foo4 = values(foo4)

print(str(r_df))

# Reshape the data for ggplot
plotData = melt(r_df, id.vars = c('x','y'))

ggplot(aes(x = x, y = y), data = plotData) +
geom_tile(aes(fill = value)) + facet_wrap(~ variable) +
scale_fill_gradient(low = 'white', high = 'blue') +
coord_equal()

#Note by ALOBO: while this is the most aesthetic solution, would it
work for large raster objects?


#Based on comment by Robert Hijmans:
#Requires raster version >= 1.8-33 and rasterVis

s <- stack(foo1,foo2, foo3,foo4)
spplot(s)
theme_set(theme_bw())
gplot(s) + geom_tile(aes(fill = value)) + facet_wrap(~ variable) +
   scale_fill_gradient(low = 'white', high = 'blue') + coord_equal()

#Note by ALOBO: assumed it automatically makes a sampleRegular() for
large raster objects?
#In such a case, it would combine aesthetic and performance

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] sp::overlay() not conserving row.names ?

[Agustin Lobo]

I'm using sp::overlay between a SpPolDF and a SpPointsDF:

> delme <- overlay(rndp,q1km, fn = mean)
> class(q1km)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
> class(rndp)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"
> class(delme)
[1] "data.frame"

but I'm confused by the fact that the row.names(delme) are not the same
as the row.names(q1km):
> delme[1:10,]
  ID   id  val
X1  356. 357. 643.3162
X2  454. 455. 653.1833
X3  305. 306. 101.8930
X4  486.7143 487.7143 448.4612
X5  456.2500 457.2500 544.1646
X6  255.4000 256.4000 139.7507
X7  403.6667 404.6667 659.8114
X8  462.5000 463.5000 537.7209
X9  629.2000 630.2000 399.1040
X10 416.5000 417.5000 278.2879
> q1km@data[1:10,]
   idname ID   XMIN   XMAX   YMIN   YMAX
X0 X0  0 -8.994 -8.594 43.595 43.995
X1 X1  1 -8.594 -8.194 43.595 43.995
X2 X2  2 -8.194 -7.794 43.595 43.995
X3 X3  3 -7.794 -7.394 43.595 43.995
X4 X4  4 -7.394 -6.994 43.595 43.995
X5 X5  5 -6.994 -6.594 43.595 43.995
X6 X6  6 -6.594 -6.194 43.595 43.995
X7 X7  7 -6.194 -5.794 43.595 43.995
X8 X8  8 -5.794 -5.394 43.595 43.995
X9 X9  9 -5.394 -4.994 43.595 43.995

Is this the way it has to be? If it is, then it's kind of inconvenient
as I would like to join the new table "delme" to q1km@data.
It seems I just have to subtract 1 to the numeric part of
row.names(delme), but would like to make sure.

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] sp::overlay() not conserving row.names ?

[Agustin Lobo]

Yes!

I'm doing:
> delme2 <- over(q1km, rndp,fn = mean)

And now not only the row.names are correct and the order of the 2
arguments makes more sense, but the dimensions of the result are also
the same
that that of the first argument
(which implies rows of NA for those polygons with no points: the user
can get rid of
these rows with na.omit() later, for me having the NA is a lot better):

> delme2[1:3,]
ID  val
0 356. 760.5662
1 454. 541.5580
2 305. 467.0864
> q1km@data[1:3,]
  ID   XMIN   XMAX   YMIN   YMAX
0  0 -8.994 -8.594 43.595 43.995
1  1 -8.594 -8.194 43.595 43.995
2  2 -8.194 -7.794 43.595 43.995

> dim(delme2)
[1] 299   3
> dim(q1km@data)
[1] 299   6

I can just join this table to the one of q1km:
> delme2 = cbind(idname=row.names(delme2),delme2)
> delme2[1:3,]
  idname   ID  val
0  0 356. 760.5662
1  1 454. 541.5580
2  2 305. 467.0864

> q1kmv2 = q1km
> q1kmv2 <- mijoin(q1kmv2,delme2[,-2], by.x=1,by.y=1)
> q1kmv2@data[1:3,]
  ID   XMIN   XMAX   YMIN   YMAX idname  val
0  0 -8.994 -8.594 43.595 43.995  0 760.5662
1  1 -8.594 -8.194 43.595 43.995  1 541.5580
2  2 -8.194 -7.794 43.595 43.995  2 467.0864
3  3 -7.794 -7.394 43.595 43.995  3 487.6461

Thanks!
Agus

2011/7/6 Edzer Pebesma :
> Agus, sp::overlay will be deprecated at some stage in favour of
> sp::over, which is a more consistent and complete approach to the same
> problem. For instance, overlay(x,y) would do the same as overlay(y,x),
> which is not good, if you think about it a bit longer (my own mistake,
> long time ago).
>
> Could you please check if sp::over has the same problems? See also
> http://cran.r-project.org/web/packages/sp/vignettes/over.pdf
>
> On 07/06/2011 05:28 PM, Agustin Lobo wrote:
>> I'm using sp::overlay between a SpPolDF and a SpPointsDF:
>>
>>> delme <- overlay(rndp,q1km, fn = mean)
>>> class(q1km)
>> [1] "SpatialPolygonsDataFrame"
>> attr(,"package")
>> [1] "sp"
>>> class(rndp)
>> [1] "SpatialPointsDataFrame"
>> attr(,"package")
>> [1] "sp"
>>> class(delme)
>> [1] "data.frame"
>>
>> but I'm confused by the fact that the row.names(delme) are not the same
>> as the row.names(q1km):
>>> delme[1:10,]
>>           ID       id      val
>> X1  356. 357. 643.3162
>> X2  454. 455. 653.1833
>> X3  305. 306. 101.8930
>> X4  486.7143 487.7143 448.4612
>> X5  456.2500 457.2500 544.1646
>> X6  255.4000 256.4000 139.7507
>> X7  403.6667 404.6667 659.8114
>> X8  462.5000 463.5000 537.7209
>> X9  629.2000 630.2000 399.1040
>> X10 416.5000 417.5000 278.2879
>>> q1km@data[1:10,]
>>    idname ID   XMIN   XMAX   YMIN   YMAX
>> X0     X0  0 -8.994 -8.594 43.595 43.995
>> X1     X1  1 -8.594 -8.194 43.595 43.995
>> X2     X2  2 -8.194 -7.794 43.595 43.995
>> X3     X3  3 -7.794 -7.394 43.595 43.995
>> X4     X4  4 -7.394 -6.994 43.595 43.995
>> X5     X5  5 -6.994 -6.594 43.595 43.995
>> X6     X6  6 -6.594 -6.194 43.595 43.995
>> X7     X7  7 -6.194 -5.794 43.595 43.995
>> X8     X8  8 -5.794 -5.394 43.595 43.995
>> X9     X9  9 -5.394 -4.994 43.595 43.995
>>
>> Is this the way it has to be? If it is, then it's kind of inconvenient
>> as I would like to join the new table "delme" to q1km@data.
>> It seems I just have to subtract 1 to the numeric part of
>> row.names(delme), but would like to make sure.
>>
>> Agus
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> --
> Edzer Pebesma
> Institute for Geoinformatics (ifgi), University of Münster
> Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
> 8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
> http://www.52north.org/geostatistics      e.pebe...@wwu.de
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] raster: croping and aggregating

[Agustin Lobo]

I have 2 raster objects

> show(CODEQP)
class   : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent  : 438048.4, 439048.4, 4633415, 4634415  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +ellps=intl +units=m +no_defs
values  : in memory
min value   : 1
max value   : 1

> dem = raster("mde10mny.tif")
> show(dem)
class   : RasterLayer
dimensions  : 2333, 2788, 6504404  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent  : 435926.9, 463806.9, 4613063, 4636393  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values  : /media/Iomega_HDD/DIPU2011/RVISITANTS/mde10mny.tif
min : ?
max : ?

and want to make a thrid object "dem100" perfectly aligned to CODEQP.

I've tried with crop() but, as noted in the help page, the input
extent is actually modified to get aligned to the input
raster pixels and I actually would like the opposite:

> dem100 = crop(dem,extent(CODEQP))
> show(dem100)
class   : RasterLayer
dimensions  : 100, 100, 1  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent  : 438046.9, 439046.9, 4633413, 4634413  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values  : in memory
min value   : 549.4215
max value   : 627.9377

> show(CODEQP)
class   : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent  : 438048.4, 439048.4, 4633415, 4634415  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +ellps=intl +units=m +no_defs
values  : in memory
min value   : 1
max value   : 1

Note they are not identical.

Is there any way (besides using grass from within R) to get dem100 to
have the same geometry (extent and pixel size) than CODEQP?

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] create raster using extent and resolution?

[Agustin Lobo]

Hi!

Is it possible to create a raster object using extent and resolution,
thus having
raster() to calculate the required nb. of rows and columns? Something like,

a = raster (ext=miextent, res=100)

(Obviously the nb. of pixels would have to be rounded)

Thanks!

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] create raster using extent and resolution?

[Agustin Lobo]

Barry,

I tried that first, but by some reason (my raster version?), the
result here is different than yours:
> CODEQP <- raster(extent(PNMlimit))
> show(CODEQP)
class   : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 2586.382, 1886.275  (x, y)
extent  : 438048.4, 463912.3, 4615552, 4634415  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values  : none

> res(CODEQP)<-c(100,100)
> show(CODEQP)
class   : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent  : 438048.4, 439048.4, 4633415, 4634415  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values  : none

Note res() keeps the same nb. of pixels and changes (a lot) the extent.
In your example:
> a = raster(e)
> show(a)
class   : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 320, 460  (x, y)
extent  : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values  : none

> res(a)=c(100,100)
> show(a)
class   : RasterLayer
dimensions  : 10, 10, 100  (nrow, ncol, ncell)
resolution  : 100, 100  (x, y)
extent  : 178400, 179400, 333000, 334000  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values  : none

Agus

> sessionInfo()
R version 2.13.1 (2011-07-08)
Platform: i486-pc-linux-gnu (32-bit)

locale:
 [1] LC_CTYPE=en_US.UTF-8  LC_NUMERIC=C
LC_TIME=en_US.UTF-8
 [4] LC_COLLATE=en_US.UTF-8LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8  LC_NAME=en_US.UTF-8
LC_ADDRESS=en_US.UTF-8
[10] LC_TELEPHONE=en_US.UTF-8  LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] stats graphics  grDevices utils datasets  methods   base

other attached packages:
[1] rgdal_0.6-31  raster_1.8-42 sp_0.9-83 rkward_0.5.6

loaded via a namespace (and not attached):
[1] grid_2.13.1 lattice_0.19-30 tools_2.13.1


Agus


2011/7/15 Barry Rowlingson :
> On Fri, Jul 15, 2011 at 12:00 PM, Agustin Lobo  wrote:
>> Hi!
>>
>> Is it possible to create a raster object using extent and resolution,
>> thus having
>> raster() to calculate the required nb. of rows and columns? Something like,
>>
>> a = raster (ext=miextent, res=100)
>>
>> (Obviously the nb. of pixels would have to be rounded)
>
> In two steps: create the raster and then set the resolution using "res<-"
>
>> e
> class       : Extent
> xmin        : 178400
> xmax        : 181600
> ymin        : 329400
> ymax        : 334000
>> a = raster(e); res(a)=c(100,100)
>> a
> class       : RasterLayer
> dimensions  : 46, 32, 1472  (nrow, ncol, ncell)
> resolution  : 100, 100  (x, y)
> extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> values      : none
>
> - it chooses the dimensions, rounding if necessary.
>
> Barry
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] initGRASS() under linux

[Agustin Lobo]

Is it necessary to start R from within the GRASS console to use
spgrass6 under linux?
My understanding was that initGRASS() would let you set the grass environment
not having GRASS actually running.

Thanks

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] empty brick?

[Agustin Lobo]

I have a process that creates a set of (large) raster layers.
Does it make sense defining an empty brick object
and then putting each layer as band of the brick?
Or is it better saving each layer and then make the brick?
Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] overlay with 2 bricks?

[Agustin Lobo]

I'm trying to combine 2 operations into one to speed up a process.
It's simple calculating a median raster across layers of a brick object, but
the values that must be considered NA are defined by a second brick.
Normally, we first set the NA in the brick and then calculate the median:

fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
decadaN_v2 <- overlay(decadaN, decadaS, fun=fun1, overwrite=TRUE)
a <- overlay(decadaN_v2, fun=median, na.rm=T, overwrite=TRUE)

But this is too slow (we actually try to do this for 36 bricks within
a loop). So I tried:

mimedian <- function(x,y) {
 x[y!=248 & y!=232] <- NA
 median(x,na.rm=T)
}

And then try to use overlay only once:

a <- overlay(decadaN, decadaS, fun=mimedian, na.rm=T, overwrite=TRUE)

but raster complains that the function is not vectorized.
Is there any way around? I suspect the problem is that overlay works
across layers if one brick is provided,
but not if 2 bricks are provided. In other words, we want x in
mimedian to be a vector of 1 pixel across layers but
actually is a brick object.

Agus

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] overlay with 2 bricks?

[Agustin Lobo]

I think the problem is that I do not quite understand
the relationship between the objects passed within overlay (s1 and s2 in
your example) and what
the function within overlay is actually getting.

In your example,
x3 <- overlay(s1, s2, fun=fun3)

s1 and s3 are stacks of 10x10x3 but what is fun3 getting?
A matrix? Note we want to calculate the median across layers,
thus vectors of 3 elements for each pixel. As you use
apply(x,1,median,na.rm=TRUE)

it seems that x within fun3 is a matrix of 10x10 rows and 3 columns,
is this correct?

Also, why are you doing
fun3(s1[1:3], s2[1:3])

instead of
fun3(s1, s2)
?

Agus

2011/7/27 Robert J. Hijmans :
> Agus,
>
> This problem arises because you are combing what are (for calc anyway)
> different types of functions. The second part is called with apply, the
> first part not. I think you can fix that using my 'fun3' (see below).
>
>
> r1 <- r2 <- raster(nc=10, nr=10)
> r1[] <- round(runif(ncell(r1))) * 248
> r2[] <- round(runif(ncell(r2))) * 232
>
> s1 <- stack(r1, r1, r2)
> s2 <- stack(r2, r2, r1)
>
> fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
>
> fun2 <- function(x,y) {
> x[y!=248 & y!=232] <- NA
> median(x,na.rm=T)
> }
>
> fun3 <- function(x,y) {
> x[y!=248 & y!=232] <- NA
> apply(x,1,median,na.rm=TRUE)
> }
>
> fun1(s1[1:3], s2[1:3])
> fun2(s1[1:3], s2[1:3]) # not good. only returns a single value
> fun3(s1[1:3], s2[1:3]) # looks good
>
> x1 <- overlay(s1, s2, fun=fun1)
> x2 <- overlay(s1, s2, fun=fun2) # not good
> x3 <- overlay(s1, s2, fun=fun3) # OK
>
>
>
> OK means that no error is thrown. I have not checked if the results are
> meaningful, but I expect they are.
>
> Robert
>
>
>
> On Tue, Jul 26, 2011 at 8:27 AM, Agustin Lobo  wrote:
>
>> I'm trying to combine 2 operations into one to speed up a process.
>> It's simple calculating a median raster across layers of a brick object,
>> but
>> the values that must be considered NA are defined by a second brick.
>> Normally, we first set the NA in the brick and then calculate the median:
>>
>> fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
>> decadaN_v2 <- overlay(decadaN, decadaS, fun=fun1, overwrite=TRUE)
>> a <- overlay(decadaN_v2, fun=median, na.rm=T, overwrite=TRUE)
>>
>> But this is too slow (we actually try to do this for 36 bricks within
>> a loop). So I tried:
>>
>> mimedian <- function(x,y) {
>>  x[y!=248 & y!=232] <- NA
>>  median(x,na.rm=T)
>> }
>>
>> And then try to use overlay only once:
>>
>> a <- overlay(decadaN, decadaS, fun=mimedian, na.rm=T, overwrite=TRUE)
>>
>> but raster complains that the function is not vectorized.
>> Is there any way around? I suspect the problem is that overlay works
>> across layers if one brick is provided,
>> but not if 2 bricks are provided. In other words, we want x in
>> mimedian to be a vector of 1 pixel across layers but
>> actually is a brick object.
>>
>> Agus
>>
>
>        [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] overlay with 2 bricks?

[Agustin Lobo]

ok, then, as I understand it,
given a raster stack  or brick s1 with dimensions n*c*b
(thus, 3D, being n,c,b, respectively, nb of rows, columns and layers),
overlay(s1,fun(x),...) would reshape s1 and pass it to fun as  a 2D
matrix x of dimensions
(n*c) rows and b layers,
and
given r1 as a n*c raster layer
(being n,c, respectively, nb of rows and columns),
overlay(r1,fun(x),...) would reshape r1 and pass it to fun as a vector
x of dimensions
(n*c) rows and 1 layers,

Correct?

Thus, if we want to apply an eigenvector matrix u of dimensions b*b to
s1 of dimensions (n*c*b),
we could just use

funortho <- function(x,u){ x%*%u}

a <- overlay(s1,u,fun=funortho(x,u))

as x would have dimensions (n*c) rows and b columns

Correct?

Thanks!

Agus

2011/7/27 Robert J. Hijmans :
> On Wed, Jul 27, 2011 at 2:23 PM, Agustin Lobo  wrote:
>
>> I think the problem is that I do not quite understand
>> the relationship between the objects passed within overlay (s1 and s2 in
>> your example) and what
>> the function within overlay is actually getting.
>>
>
> It gets a matrix where the columns are layers and the rows are cells.
>
>
>>
>> In your example,
>> x3 <- overlay(s1, s2, fun=fun3)
>>
>> s1 and s3 are stacks of 10x10x3 but what is fun3 getting?
>> A matrix?
>
>
> yes, two matrices of 100*3  (or a chunk thereof if the raster is very large)
>
>  Note we want to calculate the median across layers,
>
>> thus vectors of 3 elements for each pixel. As you use
>> apply(x,1,median,na.rm=TRUE)
>>
>> it seems that x within fun3 is a matrix of 10x10 rows and 3 columns,
>> is this correct?
>>
>> Also, why are you doing
>> fun3(s1[1:3], s2[1:3])
>>
>>
> Because I am trying it for three cells, using the type of matrix that calc
> gets:  s1[1:3]
>
>
>> instead of
>> fun3(s1, s2)
>> ?
>>
>>
> Because fun3 may not know what to do with a RasterStack.
>
>
>
>> Agus
>>
>>
> Hth, Robert
>
>
>
>> 2011/7/27 Robert J. Hijmans :
>> > Agus,
>> >
>> > This problem arises because you are combing what are (for calc anyway)
>> > different types of functions. The second part is called with apply, the
>> > first part not. I think you can fix that using my 'fun3' (see below).
>> >
>> >
>> > r1 <- r2 <- raster(nc=10, nr=10)
>> > r1[] <- round(runif(ncell(r1))) * 248
>> > r2[] <- round(runif(ncell(r2))) * 232
>> >
>> > s1 <- stack(r1, r1, r2)
>> > s2 <- stack(r2, r2, r1)
>> >
>> > fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
>> >
>> > fun2 <- function(x,y) {
>> > x[y!=248 & y!=232] <- NA
>> > median(x,na.rm=T)
>> > }
>> >
>> > fun3 <- function(x,y) {
>> > x[y!=248 & y!=232] <- NA
>> > apply(x,1,median,na.rm=TRUE)
>> > }
>> >
>> > fun1(s1[1:3], s2[1:3])
>> > fun2(s1[1:3], s2[1:3]) # not good. only returns a single value
>> > fun3(s1[1:3], s2[1:3]) # looks good
>> >
>> > x1 <- overlay(s1, s2, fun=fun1)
>> > x2 <- overlay(s1, s2, fun=fun2) # not good
>> > x3 <- overlay(s1, s2, fun=fun3) # OK
>> >
>> >
>> >
>> > OK means that no error is thrown. I have not checked if the results are
>> > meaningful, but I expect they are.
>> >
>> > Robert
>> >
>> >
>> >
>> > On Tue, Jul 26, 2011 at 8:27 AM, Agustin Lobo 
>> wrote:
>> >
>> >> I'm trying to combine 2 operations into one to speed up a process.
>> >> It's simple calculating a median raster across layers of a brick object,
>> >> but
>> >> the values that must be considered NA are defined by a second brick.
>> >> Normally, we first set the NA in the brick and then calculate the
>> median:
>> >>
>> >> fun1 <- function(x,y) { x[y!=248 & y!=232] <- NA; return(x)}
>> >> decadaN_v2 <- overlay(decadaN, decadaS, fun=fun1, overwrite=TRUE)
>> >> a <- overlay(decadaN_v2, fun=median, na.rm=T, overwrite=TRUE)
>> >>
>> >> But this is too slow (we actually try to do this for 36 bricks within
>> >> a loop). So I tried:
>> >>
>> >> mimedian <- function(x,y) {
>> >>  x[y!=248 & y!=232] <- NA
>> >>  median(x,na.rm=T)
>> >> }
>> >>
>> >> And then try to use overlay only once:
>> >>
>> >> a <- overlay(decadaN, decadaS, fun=mimedian, na.rm=T, overwrite=TRUE)
>> >>
>> >> but raster complains that the function is not vectorized.
>> >> Is there any way around? I suspect the problem is that overlay works
>> >> across layers if one brick is provided,
>> >> but not if 2 bricks are provided. In other words, we want x in
>> >> mimedian to be a vector of 1 pixel across layers but
>> >> actually is a brick object.
>> >>
>> >> Agus
>> >>
>> >
>> >        [[alternative HTML version deleted]]
>> >
>> > ___
>> > R-sig-Geo mailing list
>> > R-sig-Geo@r-project.org
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>
>        [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo