raster (was readGDAL) loses datum
Just to be clear, rgdal::getProjectionRef has a bunch of conditionals, and the part I showed is the part that's relevant for this particular case. I guess I could just set the environment variable myself: Sys.setenv(OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = "NO") Oliver
On Wed, Nov 7, 2012 at 11:40 AM, Oliver Soong <osoong+r at gmail.com> wrote:
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