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

[Roger Bivand]

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