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

[Oliver Soong]

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
>