[R-sig-Geo] readGDAL loses datum
Roger Bivand
Roger.Bivand at nhh.no
Sat Nov 3 20:36:03 CET 2012
On Sat, 3 Nov 2012, Oliver Soong wrote:
> I did notice another thing. I've attached the .prj file I've been
> using along with a raster exported in that projection using ArcGIS.
>
>> system("gdalinfo -proj4 Arc.tif")
> [...]
> 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
> +datum=NAD83 +units=m +no_defs '
> [...]
>> readGDAL("Arc.tif")@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=NAD83
> +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
>
> There are similar issues with extra Proj4 parameters showing up with
> +proj=longlat +datum=WGS84 +no_defs, although these seem to be added
> on read and stripped on write, somehow.
>
> wgs84_1 <- SpatialGridDataFrame(GridTopology(c(0, 0), c(1, 1), c(1,
> 1)), data.frame(band = 1), CRS("+proj=longlat +datum=WGS84 +no_defs"))
> writeGDAL(wgs84_1, "wgs84_1.tif")
> system("gdalinfo -proj4 wgs84_1.tif")
> wgs84_2 <- readGDAL("wgs84_1.tif")
> wgs84_2 at proj4string
> writeGDAL(wgs84_2, "wgs84_2.tif")
> system("gdalinfo -proj4 wgs84_2.tif")
>
> Any idea where the extra +ellps +towgs84 terms are coming from?
Yes - they are part of the safety improvements made in GDAL 1.8.0 to make
sure that the spatial reference definitions are unambiguous. Extra tags
are not a problem as long as they are consistent - your problem was that
Arc doesn't like not being given a datum. Now we know how to do that, even
though +towgs84=0,0,0 implies both of NAD83 and WGS84.
Roger
>
> Oliver
>
>
> On Sat, Nov 3, 2012 at 6:59 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>> On Fri, 2 Nov 2012, Roger Bivand wrote:
>>
>>> 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?
>>
>>
>> I've committed to the rgdal R-Forge project a user argument to relevant
>> functions in rgdal for setting the behaviour to datum-preserving, either on
>> a case-by-case or global level. If the environment variable is present, it
>> will have precedence and will not be overwritten. I'd welcome reports from
>> those who can try out the source checked out from R-forge.
>>
>> In your case, I now see:
>>
>>> proj4string(readGDAL("img1.tif"))
>>
>> img1.tif has GDAL driver GTiff
>> and has 1 rows and 1 columns
>> [1] " +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"
>>>
>>> proj4string(readGDAL("img1.tif",
>>
>> OVERRIDE_PROJ_DATUM_WITH_TOWGS84=FALSE))
>> img1.tif has GDAL driver GTiff
>> and has 1 rows and 1 columns
>> [1] " +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"
>>
>> However, there does not seem to be a clear way to auto-detect and set the
>> switch. The package also has a cached variable, so setting:
>>
>> set_OVERRIDE_PROJ_DATUM_WITH_TOWGS84(FALSE)
>>
>> in an R session will use that value when raster projections are read until
>> the setting is changed. The argument will not be used for GDAL < 1.8.0,
>> because its use first appeared there. It will also not be used if an
>> environment variable "OVERRIDE_PROJ_DATUM_WITH_TOWGS84" is found, to avoid
>> overwriting its value.
>>
>> This is the result of discussions on the gdal-dev list, thread starting at:
>>
>> http://lists.osgeo.org/pipermail/gdal-dev/2012-November/034550.html
>>
>> and details in:
>>
>> http://trac.osgeo.org/gdal/ticket/4880
>>
>> Roger
>>
>>
>>>
>>> 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
>>
>
--
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