[R-sig-Geo] raster (was readGDAL) loses datum

[Oliver Soong]

So I'm also using the raster package, and the datum problems with
rgdal are inherited.  For those who haven't been following this
thread, there is a workaround as of rgdal 0.7.22 and setting
OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = FALSE.  Is there any way to have
raster respect this setting?  For small rasters, I know I can use
readGDAL and coerce the SpatialGridDataFrame to a RasterLayer, but
this requires being able to load the entire raster into memory.

R 2.15.1 32-bit, rgdal 0.7.22, raster 2.0.12, Windows XP/7.

Test code, translated from the rgdal example:

ext <- extent(-2100050, -2099950, 1199950, 1200050)
p4s <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0
+y_0=0 +datum=NAD83 +units=m +no_defs"
img1 <- raster(nrows = 1, ncols = 1, crs = p4s, ext = ext)
values(img1) <- 0
img1.file <- file.path(tempdir(), "img1.tif")
writeRaster(img1, img1.file, "GTiff", overwrite = TRUE)
img2 <- raster(img1.file)
img3 <- readGDAL(img1.file, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = FALSE)
projection(img1)
projection(img2)
projection(img3)

Thanks,
Oliver


On Fri, Nov 2, 2012 at 1:03 AM, Oliver Soong <osoong+r at gmail.com> wrote:
> R 2.15.1 32-bit, rgdal 0.7.20, Windows 7.
>
> grid <- GridTopology(c(-2100000, 1200000), c(100, 100), c(1, 1))
> p4s <- CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96
> +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80
> +towgs84=0,0,0")
> img1 <- SpatialGridDataFrame(grid, data.frame(band1 = 1), p4s)
> img1.file <- file.path(tempdir(), "img1.tif")
> writeGDAL(img1, img1.file)
> img2 <- readGDAL(img1.file)
> img2.file <- file.path(tempdir(), "img2.tif")
> writeGDAL(img2, img2.file)
> img1 at proj4string
> img2 at proj4string
>
> For me, img1 at proj4string has +datum=NAD83 and img2 at proj4string does
> not.  Not surprisingly, if I look at both files in Arc, img1 has a
> defined datum and img2 does not.
>
> Am I doing anything wrong?
>
> Oliver