readGDAL loses datum

I agree it seems to be happening when converting WKT to Proj4.
However, is this more of a GDAL bug?
system(paste("gdalinfo -proj4", img1.file))
[...]
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, 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

readGDAL loses datum

Thread (5 messages)