readGDAL loses datum
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