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

[Oliver Soong]

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