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

Roger Bivand Roger.Bivand at nhh.no
Wed Feb 10 18:50:05 CET 2010 

On Wed, 10 Feb 2010, Agustin Lobo wrote:

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

And there is a similar mechanism in sp. However, please find out and 
report to the plugin author why the plugin is flipping GTiffs.

Roger

> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> 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.Bivand at nhh.no

More information about the R-sig-Geo mailing list