[R-sig-Geo] readGDAL loses datum

[Roger Bivand]

On Fri, 2 Nov 2012, Oliver Soong wrote:

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

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:

> img2 at proj4string
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:

> system(paste("gdalinfo", img1.file))
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:

> system(paste("gdalinfo", img2.file))
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)

> system(paste("gdalinfo", 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
>
> _______________________________________________
> 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