[R-sig-Geo] doubt with readGDAL

[Barry Rowlingson]

2009/10/4 "José M. Blanco Moreno" <jmblanco at ub.edu>:
> 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
>
> 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).

 I just cut and pasted your example - firstly there's no tiff in that
zip, there's a .asc (Arc GRID) file instead. I read that in with
readGDAL, and got:

> bbox(cata)
            min        max
x -4.166183e-04   5.000417
y -1.577894e+01 -10.778109
> bbox(cata2)
            min        max
x -4.166183e-04   5.000417
y -1.577894e+01 -10.778109

 which looks okay. So maybe its to do with the tiff file which you
must have made somehow - projections maybe?

Barry