[R-sig-Geo] Additional info on "Reversed raster after readGDAL()" (Fwd: Overlay spatial lines on raster)

Agustin Lobo alobolistas at gmail.com
Wed Feb 10 18:45:11 CET 2010 

FYI
This problem was actually identified because I was getting wrong
results at calculating
averaged raster values for intersecting line segments.
I include the last message sent by Robert. His trick is a good
workaround by now.
Agus


---------- Forwarded message ----------
From: "Robert J. Hijmans" <r.hijmans at gmail.com>
Date: Wed, 10 Feb 2010 09:01:30 -0800
Subject: Re: [R-sig-Geo] Overlay spatial lines on raster
To: Agustin Lobo <Agustin.Lobo at ija.csic.es>

Ah!

The fundamentals seems to be beyond my reach (in (r)gdal) but the
simple raster patch if you *know* this is happening is:

r = flip(r)

Robert


On Wed, Feb 10, 2010 at 4:19 AM, Agustin Lobo <Agustin.Lobo at ija.csic.es> wrote:
> 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]]
>>>>>>>>>>>
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> R-sig-Geo mailing list
>>>>>>>>>>> R-sig-Geo at stat.math.ethz.ch
>>>>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>>>>
>>>>> --
>>>>> 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
>



--
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

More information about the R-sig-Geo mailing list