[R-sig-Geo] readGDAL loses datum
Oliver Soong
osoong+r at gmail.com
Fri Nov 2 21:14:58 CET 2012
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, Nov 2, 2012 at 2:24 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
More information about the R-sig-Geo
mailing list