raster (was readGDAL) loses datum
Not surprisingly, you guys are ahead of me. The dev version does indeed work as expected. I do like how the necessary changes to raster were submitted months before the changes to rgdal. Thanks for all your patience and effort. Oliver
On Thu, Nov 8, 2012 at 1:46 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
On Wed, 7 Nov 2012, Oliver Soong wrote:
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")
I assume that you are running released raster on CRAN - have you tried installing the development version from R-forge, in which many of these steps have been taken? Could you check that first? The conditionals are there to try to avoid trashing OVERRIDE_PROJ_DATUM_WITH_TOWGS84 possibly set by other software, among other things. The development raster code goes through rgdal::GDALinfo() which respects the rgdal option. Roger
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
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no