[R-sig-Geo] readGDAL loses datum

Roger Bivand Roger.Bivand at nhh.no
Fri Nov 2 22:51:04 CET 2012 

On Fri, 2 Nov 2012, Oliver Soong wrote:

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

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

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

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