Message-ID: <CAPQ=hVJRfRLn_G1s3efGTCUktckr19viWdQAFCFbrdWWXQ6j2Q@mail.gmail.com>
Date: 2012-11-07T19:40:53Z
From: Oliver Soong
Subject: raster (was readGDAL) loses datum
In-Reply-To: <CANtt_hy5cqpjt7TwbYvyZNugZVF=s=dbrrDK6gVVcnA-T7mSjw@mail.gmail.com>
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
>
>
>