[R-sig-Geo] doubt with readGDAL

[Roger Bivand]

On Sun, 4 Oct 2009, "José M. Blanco Moreno" wrote:

> Dear R-users community,
> I have been trying to figure out on my own, but no way...
> I am trying to read and decimate a geotiff:
>
> library(maptools)
> library(spatial)
> library(sp)
> library(rgdal)
>
> download.file(url='http://srtm.geog.kcl.ac.uk/portal/srtm41/srtm_data_arcascii/srtm_37_04.zip',destfile='srtm_37_04.zip')
> shell('unzip -o srtm_37_04.zip')
> cata <- 
> readGDAL('srtm_37_04.tif',offset=c(2900,0),region.dim=c(2600,5000),output.dim=c(260,500))
> cata2 <- cata
> cata2 at data$band1 <- ifelse(is.na(cata2 at data$band1),-1,cata2 at data$band1)
> cata2 at data$band1 <- 
> c(50,100,250,500,1000,1500,2000,3000)[cut(cata at data$band1,c(0,50,100,250,500,1000,1500,2000,3000))]
> image(cata2,col=gray(seq(.9,0.6,length=9)))
> contour(cata2,add=T,levels=c(0,50,100,250,500,1000,1500,2000,3000),col=gray(.7),drawlabels=F)
> axis(1)
> axis(2)
>
> But the coordinates don't come out: they are translated in the y-axis *a 
> lot*; bbox(cata) returns:
>           min        max
> x -4.166183e-04   5.000417
> y -1.577894e+01 -10.778109
>
> While in the non-decimated version they are:
>           min      max
> x -0.0004166183  4.16625
> y 40.4170835754 42.58375

Thanks for a clear report - the y range and the cell sizes were wrong when 
all three of offset=, region.dim= and output.dim= were used together. With 
current rgdal 0.6-18, using the Arc ASCII file (not GTiff), I see for the 
undecimated subscene:

> gridparameters(catb)
   cellcentre.offset     cellsize cells.dim
x      4.839851e-08 0.0008333333      5000
y      4.041750e+01 0.0008333333      2600
> bbox(catb)
             min      max
x -0.0004166183  4.16625
y 40.4170835754 42.58375

but for the decimated subscene:

> gridparameters(cata)
   cellcentre.offset   cellsize cells.dim
x       0.004584215 0.01000167       500
y     -15.769325078 0.01923397       260
> bbox(cata)
             min        max
x -4.166183e-04   5.000417
y -1.577894e+01 -10.778109


In forthcoming 0.6-19 (committed to the sourceforge rgdal repository) and 
the undecimated subscene:

> gridparameters(catb)
   cellcentre.offset     cellsize cells.dim
x      4.839851e-08 0.0008333333      5000
y      4.041750e+01 0.0008333333      2600
> bbox(catb)
             min      max
x -0.0004166183  4.16625
y 40.4170835754 42.58375

with the decimated subscene now corresponding:

> gridparameters(cata)
   cellcentre.offset    cellsize cells.dim
x       0.003750048 0.008333333       500
y      40.421250242 0.008333333       260
> bbox(cata)
             min      max
x -0.0004166183  4.16625
y 40.4170835754 42.58375

This is with this subscene on this data set - I would be very grateful if 
anyone could check from a source install after checking out from 
sourceforge for other, including Southern hemisphere decimated subscenes 
(there is a sign wrinkle). I'll wait for feedback before submitting to 
CRAN.

Again, thanks for a helpful report!

Roger

>
> Could anyone indicate a workaround to preserve the coordinates of the grid?
> Thank you for any information on this issue (which I suppose is generated by 
> my own incompetence).
>
> José M.
>
> sessionInfo()
> R version 2.8.1 (2008-12-22)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=Spanish_Spain.1252;LC_CTYPE=Spanish_Spain.1252;LC_MONETARY=Spanish_Spain.1252;LC_NUMERIC=C;LC_TIME=Spanish_Spain.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base 
> other attached packages:
> [1] rgdal_0.6-8     spatial_7.2-47  maptools_0.7-25 sp_0.9-37 
> foreign_0.8-36
>
> loaded via a namespace (and not attached):
> [1] grid_2.8.1      lattice_0.17-25
>
>
>

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