Skip to content

raster (was readGDAL) loses datum

6 messages · Robert J. Hijmans, Oliver Soong, Roger Bivand

Oliver Soong
# 
So I'm also using the raster package, and the datum problems with
rgdal are inherited.  For those who haven't been following this
thread, there is a workaround as of rgdal 0.7.22 and setting
OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = FALSE.  Is there any way to have
raster respect this setting?  For small rasters, I know I can use
readGDAL and coerce the SpatialGridDataFrame to a RasterLayer, but
this requires being able to load the entire raster into memory.

R 2.15.1 32-bit, rgdal 0.7.22, raster 2.0.12, Windows XP/7.

Test code, translated from the rgdal example:

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, OVERRIDE_PROJ_DATUM_WITH_TOWGS84 = FALSE)
projection(img1)
projection(img2)
projection(img3)

Thanks,
Oliver
On Fri, Nov 2, 2012 at 1:03 AM, Oliver Soong <osoong+r at gmail.com> wrote:
Oliver Soong
# 
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 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:
Roger Bivand
#
On Wed, 7 Nov 2012, Oliver Soong wrote:

            

      
      
      
    

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