Skip to content

readGDAL loses datum

5 messages · Oliver Soong, Roger Bivand

Oliver Soong
# 
R 2.15.1 32-bit, rgdal 0.7.20, Windows 7.

grid <- GridTopology(c(-2100000, 1200000), c(100, 100), c(1, 1))
p4s <- CRS("+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 +ellps=GRS80
+towgs84=0,0,0")
img1 <- SpatialGridDataFrame(grid, data.frame(band1 = 1), p4s)
img1.file <- file.path(tempdir(), "img1.tif")
writeGDAL(img1, img1.file)
img2 <- readGDAL(img1.file)
img2.file <- file.path(tempdir(), "img2.tif")
writeGDAL(img2, img2.file)
img1 at proj4string
img2 at proj4string

For me, img1 at proj4string has +datum=NAD83 and img2 at proj4string does
not.  Not surprisingly, if I look at both files in Arc, img1 has a
defined datum and img2 does not.

Am I doing anything wrong?

Oliver
Roger Bivand
#
On Fri, 2 Nov 2012, Oliver Soong wrote:

            

      
      
      
    

        
No, but it isn't obvious:

p4s <- CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96
+x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs")

gives on re-reading:

            

      
      
      
    

        
CRS arguments:
  +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0
+datum=NAD27 +units=m +no_defs +ellps=clrk66
+nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat

So:

   oSRS.importFromWkt( &pszSRS_WKT );
   oSRS.exportToProj4( &pszSRS_WKT );

in RGDAL_GetProjectionRef() in src/gdal-bindings.cpp sees that input 
+datum=NAD83 is equivalent to +towgs84=0,0,0,0,0,0,0 on output. So the 
descriptions are not string-equivalent, but are equivalent through 
+towgs84=0,0,0,0,0,0,0. If you run with your p4s, on a system with 
gdalinfo:

            

      
      
      
    

        
Driver: GTiff/GeoTIFF
Files: /tmp/Rtmp9j2vOr/img1.tif
Size is 1, 1
Coordinate System is:
PROJCS["unnamed",
     GEOGCS["NAD83",
         DATUM["North_American_Datum_1983",
             SPHEROID["GRS 1980",6378137,298.2572221010002,
                 AUTHORITY["EPSG","7019"]],
             TOWGS84[0,0,0,0,0,0,0],
             AUTHORITY["EPSG","6269"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4269"]],
     PROJECTION["Albers_Conic_Equal_Area"],
     PARAMETER["standard_parallel_1",29.5],
     PARAMETER["standard_parallel_2",45.5],
     PARAMETER["latitude_of_center",23],
     PARAMETER["longitude_of_center",-96],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

...

so the simplification is happening on conversion to Proj4 on reading.

I agree that on re-export that the WKT and Proj4 versions diverge, so:

            

      
      
      
    

        
Driver: GTiff/GeoTIFF
Files: /tmp/RtmpMhPgqf/img2.tif
Size is 1, 1
Coordinate System is:
PROJCS["unnamed",
     GEOGCS["GRS 1980(IUGG, 1980)",
         DATUM["unknown",
             SPHEROID["GRS80",6378137,298.257222101],
             TOWGS84[0,0,0,0,0,0,0]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Albers_Conic_Equal_Area"],
     PARAMETER["standard_parallel_1",29.5],
     PARAMETER["standard_parallel_2",45.5],
     PARAMETER["latitude_of_center",23],
     PARAMETER["longitude_of_center",-96],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

with the correct parameters, but no datum name tag. You get around this 
manually by adding the +datum= back in:

proj4string(img2) <- CRS(paste(proj4string(img2), "+datum=NAD83"))
writeGDAL(img2, img2.file)

            

      
      
      
    

        
Driver: GTiff/GeoTIFF
Files: /tmp/RtmpMhPgqf/img2.tif
Size is 1, 1
Coordinate System is:
PROJCS["unnamed",
     GEOGCS["NAD83",
         DATUM["North_American_Datum_1983",
             SPHEROID["GRS 1980",6378137,298.2572221010002,
                 AUTHORITY["EPSG","7019"]],
             TOWGS84[0,0,0,0,0,0,0],
             AUTHORITY["EPSG","6269"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4269"]],
     PROJECTION["Albers_Conic_Equal_Area"],
     PARAMETER["standard_parallel_1",29.5],
     PARAMETER["standard_parallel_2",45.5],
     PARAMETER["latitude_of_center",23],
     PARAMETER["longitude_of_center",-96],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]

I would appeal to any programmer with a little time to see how the step 
between:

   oSRS.importFromWkt( &pszSRS_WKT );
   oSRS.exportToProj4( &pszSRS_WKT );

and the R output might be checked. The content of pszSRS_WKT is OK before 
entering importFromWkt(), but is simplified on exit from exportToProj4(). 
The comparable part of gdal/gdal-1.9.2/apps/gdalinfo.c is around lines 
263-274.

The writing operation appears to be OK from your example.

Roger

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Oliver Soong
    

    

      
      #
    

  


  

      
I agree it seems to be happening when converting WKT to Proj4.
However, is this more of a GDAL bug?

            

      
      
      
    

        
[...]
PROJ.4 string is:
'+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

It looks to me like the datum is getting dropped by:
 OSRExportToProj4( hSRS, &pszProj4 );
This being what rgdal seems to use.

I'm sadly less familiar with Proj4 than I ought to be to be talking
about this, but it strikes me that if this is indeed a bug and not a
"feature", then it would make more sense to fix
OSRExportToProj4/OGRSpatialReference::exportToProj4.

Oliver

        
On Fri, Nov 2, 2012 at 2:24 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 2 Nov 2012, Oliver Soong wrote:

        

            

      
      
      
    

        
Yes, but do we know that hSRS contains these tags, since they have been 
imported from WKT to the internal representation? The description of the 
import from WKT process suggests that it will terminate before processing 
the whole string if the description is already "complete":

http://www.gdal.org/ogr/classOGRSpatialReference.html#ab74cfc985bd05404a4c61d2d633a6343

I tried adding morphFromESRI() before exporting to Proj4, but the problem 
is not resolved. Perhaps the gdal-dev list is where to ask?

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Fri, 2 Nov 2012, Roger Bivand wrote:

        

            

      
      
      
    

        
I've committed to the rgdal R-Forge project a user argument to relevant 
functions in rgdal for setting the behaviour to datum-preserving, either 
on a case-by-case or global level. If the environment variable is present, 
it will have precedence and will not be overwritten. I'd welcome reports 
from those who can try out the source checked out from R-forge.

In your case, I now see:

            

      
      
      
    

        
img1.tif has GDAL driver GTiff
and has 1 rows and 1 columns
[1] " +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 
+ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

            

      
      
      
    

        
OVERRIDE_PROJ_DATUM_WITH_TOWGS84=FALSE))
img1.tif has GDAL driver GTiff
and has 1 rows and 1 columns
[1] " +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 +ellps=GRS80 +towgs84=0,0,0"

However, there does not seem to be a clear way to auto-detect and set the 
switch. The package also has a cached variable, so setting:

set_OVERRIDE_PROJ_DATUM_WITH_TOWGS84(FALSE)

in an R session will use that value when raster projections are read until 
the setting is changed. The argument will not be used for GDAL < 1.8.0, 
because its use first appeared there. It will also not be used if an 
environment variable "OVERRIDE_PROJ_DATUM_WITH_TOWGS84" is found, to avoid 
overwriting its value.

This is the result of discussions on the gdal-dev list, thread starting 
at:

http://lists.osgeo.org/pipermail/gdal-dev/2012-November/034550.html

and details in:

http://trac.osgeo.org/gdal/ticket/4880

Roger