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

[Oliver Soong]

I had tried that, but it unfortunately did not fix the problem with raster.

library(rgdal)
library(raster)
set_OVERRIDE_PROJ_DATUM_WITH_TOWGS84(FALSE)
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)
projection(img1)
projection(img2)
projection(img3)

I think it might have something to do with this line in
raster:::.readRasterFromGDAL:

projection(r) <- .Call("RGDAL_GetProjectionRef", x, PACKAGE = "rgdal")

By contrast, rgdal::getProjectionRef does:
Sys.setenv(OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = "NO")
res <- .Call("RGDAL_GetProjectionRef", dataset, PACKAGE = "rgdal")
Sys.unsetenv("OVERRIDE_PROJ_DATUM_WITH_TOWGS84")

Could you use rgdal::getProjectionRef, rather than directly calling
RGDAL_GetProjectionRef?

Oliver


On Wed, Nov 7, 2012 at 11:09 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
> Oliver,
>
> I think starting your script with this should do it:
>
> library(rgdal)
> set_OVERRIDE_PROJ_DATUM_WITH_TOWGS84(FALSE)
>
> Robert
>
>
>