Additional info on "Reversed raster after readGDAL()" (Fwd: Overlay spatial lines on raster)

Robert,

Sorry I could not react earlier to this mystery, I've found a clue now:
the raster is reversed in R in respect to QGIS (see attachment).

I've uploaded test2.tif (identical to test.tif but with the
appropriate projection info)
to
http://sites.google.com/site/eospansite/dummy/test2.tif

The same problem arises with raster() and with readGDAL(), so we
should CC to sig-geo at some point?

This is what I do:

r2 <- raster("test2.tif")
plot(r2)
GDALinfo("test2.tif")

rows ? ? ? ?256
columns ? ? 256
bands ? ? ? 1
origin.x ? ? ? ?424389
origin.y ? ? ? ?4636016
res.x ? ? ? 345.0078
res.y ? ? ? 194.0195
oblique.x ? 0
oblique.y ? 0
driver ? ? ?GTiff
projection ?+proj=utm +zone=31 +ellps=intl +units=m +no_defs
file ? ? ? ?test2.tif
apparent band summary:
? GDType ? ? ? ?Bmin ? ? ? Bmax
1 Float64 -4294967295 4294967295
Metadata:
AREA_OR_POINT=Area

r3 <- readGDAL("test2.tif")

test2.tif has GDAL driver GTiff
and has 256 rows and 256 columns

image(r3)

I'm not sure I'm not doing something really stupid here, I apologize
in that case!

Thanks

Agus

2010/2/1, Robert J. Hijmans <r.hijmans at gmail.com>:

Agus, I am not sure what is going on here. But this example seems to
work for me:

library(raster)
setwd('c:/downloads')
load("xter.rda")
r <- raster("test.tif")
projection(r) <- projection(xter)

# to speed things up for testing
r = crop(r, extent(423517, ?454000, 4675610 , 4686260 ))

lr = linesToRaster(xter, r, field = "ID_GRAFIC", ?progress='text')

z = zonal(r, lr, mean)
z$mean = as.numeric(z$mean) ?# this will not be necessary in the next
version

d <- xter at data
dd <- merge(d, z, by.x="ID_GRAFIC", by.y="zone", sort=F, all.x=T)
xter at data <- dd[order(dd[,1]),]

# just for testing
test = linesToRaster(xter, r, field = "mean", progress='text')

par(mfrow=c(1,3))
plot(lr)
plot(r)
lines(xter)
plot(test)

# looks good, also in spplot
x11()
spplot(xter, 'mean')

writeOGR(xter, dsn="xter",layer="xtertest",driver="ESRI Shapefile")





On Mon, Feb 1, 2010 at 1:21 AM, Agustin Lobo <Agustin.Lobo at ija.csic.es>
wrote:

Robert,

Thanks. The mechanics work, but I get inconsistent results.
This is what I do (there is also a minor problem with writeOGR() that
I'm circumventing by now):

lr = linesToRaster(xter, lr, field = "ID_GRAFIC" )
xterpdata2 = zonal(r, lr, mean)
xterpdata2 = xterpdata2[order(xterpdata2[,1]),]
xterpdata3 <- na.omit(xterpdata2)
dim(xterpdata3)
dim(xterpdata2)
xter2 <- xter
xter2data <- xter2 at data
xter2data <-
merge(xter2data,xterpdata3,by.x="ID_GRAFIC",by.y="zone",sort=F,all.x=T)
xter2data <- xter2data[order(xter2data[,1]),]
xter2 at data <- xter2data
writeOGR(xter2, dsn="xter2",layer="xter2",driver="ESRI Shapefile")
#and get an error:
#Error in writeOGR(xter2, dsn = "xter2", layer = "xter2", driver = "ESRI
Shapefile") :
# ?Can't convert columns of class: array; column names: mean
#class(xter2 at data$mean)
#[1] "array"
#which I curcunvent by:
xter2 at data$mean <- as.numeric(xter2 at data$mean)
writeOGR(xter2, dsn="xter2",layer="xter2",driver="ESRI Shapefile")
writeRaster(lr,filename="lr",format="GTiff")


Now, we check for example the result of zonal()

xterpdata2[xterpdata2[,1]==8500,]

? ? zone ? ? mean
8500 8500 4.253456

But in the QGIS we see that the values of test.tif (r in R)
are very different (~19 in test (see to the left), while field "mean" is
4.253, see attached picture).
Could be the problem is in zonal() ?

Agus


Robert J. Hijmans wrote:

Hola Agus,

Thanks for sending these data. That was very helpful. This is a whole
different (subsetting) problem for line segments consisting of two
points only (with simple solution, "drop=FALSE"). Nice to have found
it.
Your example now runs for me, although it takes three minutes on my PC
(I'll speed it up eventually, let's first make sure it works). You can
use the attached updated package (or wait till tomorrow for R-Forge).

Robert


On Sun, Jan 31, 2010 at 3:18 PM, Agustin Lobo <Agustin.Lobo at ija.csic.es>
wrote:

Thanks, actually the sourceforge version was ready when I tried again
(few
minutes ago)
Still problems. I've uploaded the files for you to try if you wish:
http://sites.google.com/site/eospansite/dummy/xter.rda
http://sites.google.com/site/eospansite/dummy/test.tif

It would be great if you could run:
load("xter.rda")
r <- raster("test.tif")
projection(r) <- projection(xter)
lr = linesToRaster(xter, r, field = "ID_GRAFIC" )

I get:
Error in subset.matrix(xyxy, !((xyxy[, 2] > maxy & xyxy[, 4] > maxy) |
?:
?subscript out of bounds
Calls: linesToRaster -> .getCols -> subset -> subset.matrix


Agus

Robert J. Hijmans wrote:

Hi Agus, Attached is the linux (is that your o.s.?) package for local
install. Robert

On Sun, Jan 31, 2010 at 2:00 AM, Agustin Lobo <alobolistas at gmail.com>
wrote:

Looks great, is there any way I can try
right now? Otherwise I'll wait
until 0.9.9-5 is ready for installation. I tried
now and only 0.9.9-4 is available.

Thanks!

Agus

Robert J. Hijmans wrote:

Hi Agus, thank you. You were right about the cause of the problem ---
I patched it in version 0.9.9-5, now on R-Forge (and should be
available for automatic install within 24 hours); the function could
use some optimization for speed but at least it works now for
SpatialLines* with a larger extent then the target RasterLayer.

# simple example
library(raster)
library(maptools)
data(wrld_simpl)
r = raster(xmn=-10, xmx=20, ymn=37, ymx=43, ncol=360, nrow=72)
r = linesToRaster(wrld_simpl, r, progress='text')
plot(r)
plot(wrld_simpl, add=T)


# example using zonal as I proposed in your original question, but
here with polygonsToRaster
# What is (roughly) the mean latitude of each country?

x = raster()
# res(x) = 0.1 ? ?# higher res, slower, but more small countries are
included
x = polygonsToRaster(wrld_simpl, x, progress='text')
y = raster(x)
y[] = rep(yFromRow(y, 1:nrow(y)), each=ncol(y))
z = zonal(y, x, mean)
res= data.frame(wrld_simpl[z[,'zone'], 'NAME'], lat=z[,2])
res[order(res[,2]), ]

Robert



On Sat, Jan 30, 2010 at 5:03 AM, Agustin Lobo <alobolistas at gmail.com>
wrote:

Thanks Robert,

The raster package is a magician's chest!
I get an error, though:

class(xter)

[1] "SpatialLinesDataFrame"
attr(,"package")
[1] "sp"

class(r)

[1] "RasterLayer"
attr(,"package")
[1] "raster"

lr = linesToRaster(xter, r, field = "ID_GRAFIC" )

[1] "A"
? rows cols
[1,] ? ?1 ? -1
[2,] ? -1 ? -1
? ? ? [,1] ? ?[,2]
[1,] 413809.3 4685685
[1] ?413786.5 4685646.9 ?413810.1 4685686.4
Error in col1:col2 : NA/NaN argument
Calls: linesToRaster -> .getCols

I think because the extent of the raster is smaller than the extent
of
the
lines

extent(xter)

class ? ? ? : Extent
xmin ? ? ? ?: 412792
xmax ? ? ? ?: 517648.1
ymin ? ? ? ?: 4624794
ymax ? ? ? ?: 4698657

extent(r)

class ? ? ? : Extent
xmin ? ? ? ?: 424389
xmax ? ? ? ?: 512711
ymin ? ? ? ?: 4636016
ymax ? ? ? ?: 4685685

which I can visualize with:

plot(extent(xter))
plot(extent(r),col="red",add=T)

The problem is that I cannot find the way to clip xter to the extent
of
r.
pruneMap() assumes a map object and
I guess that we cannot convert from Spatial Lines to map.

I've tried

delme <- xter
delme at bbox <- bbox(r)

but still get the same problem:

lr = linesToRaster(delme, r, field = "ID_GRAFIC" )

[1] "A"
? rows cols
[1,] ? ?1 ? -1
[2,] ? -1 ? -1
? [,1] [,2]
[1,] ? NA ? NA
[1] ?413786.5 4685646.9 ?413810.1 4685686.4
Error in col1:col2 : NA/NaN argument
Calls: linesToRaster -> .getCols

Agus

2010/1/29 Robert J. Hijmans <r.hijmans at gmail.com>

How about something like:

library(raster)
r = linesToRaster(xter, xterpdata, field = "ID" )
v = zonal(xterpdata, r, mean)

v now shold has the mean value of the raster for each line ID and
you
can merge those back


On Fri, Jan 29, 2010 at 9:54 AM, Agustin Lobo
<alobolistas at gmail.com>
wrote:

Partially answering myself:
This is what I've done until now:

xter <-




readOGR(dsn="/media/Transcend/PROVI/MARICEL2/x_terRiverBasin",layer="x_ter")

class(xter)

[1] "SpatialLinesDataFrame"
attr(,"package")
[1] "sp"

xterp <- coordinates(xter)
xterp <- sapply(xterp, function(x) do.call("rbind", x))
xterp <- do.call("rbind",xterp)

as told by Roger.

Now:
require(raster)
r <- raster("test.tif")
xterpdata <- xyValues(r,xterp)
xterp2 <- data.frame(utmx=xterp[,1], utmy=xterp[,2],
NO3_2007=xterpdata)
coordinates(xterp2) <- c("utmx", "utmy")

So I get to step 2, but don't see how to convert from the points
to the lines again (to get xterl). The problem is that when I did:
xterp <- sapply(xterp, function(x) do.call("rbind", x))
xterp <- do.call("rbind",xterp)

I lost the IDs of the lines. So I would need to keep the IDs of
the
lines
for each point in the result of
coordinates() to be able of going back to the lines once I've put
the
raster
values in the data slot of the xterp

I'll keep on, but any help would be much appreciated.

Agus

2010/1/29 Agustin Lobo <alobolistas at gmail.com>

Hi!

My goal is to colorcode a river network according to the values
of
a raster map (which comes from an interpolation). The spatial
lines
are
imported
from a shapefile. The raster could be imported or I could run the
interpolation in R.

As what I'm planing could be rather involved, I prefer asking
before:

1. Convert spatial lines to spatial points as discussed earlier
this
week.
2. Overlay points on the raster and get the values for the points
in
a
Spatial Points Data Frame
3. Convert the Spatial Points Data Frame into Spatial Lines Data
Frame
4. writeOGR() to shapefile.

The main question is: is the data frame of the points going to be
transferred to the lines
at step 3?

THnaks

Agus

? ? [[alternative HTML version deleted]]

--
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster

--
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster

--
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
e-mail Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster