[gdal-dev] EPSG code not recognized from GeoTiff written by GDAL

[Roger Bivand]

On Sat, 28 May 2016, Even Rouault wrote:

> On Saturday 28 May 2016 09:02:49 Roger Bivand wrote:
>> Blumentrath, Stefan <Stefan.Blumentrath <at> nina.no> writes:
>>> Hi,
>>>
>>> And thanks Even and Jukka for your answers!
>>>
>>> Based on the info you provided I investigated a bit more and it seems to
>>
>> be an issue "downstreams" indeed.
>>
>>> gdalsrsinfo epsg:25832 gave the same result for me than for Jukka (the
>>
>> result from listgeo further down,
>>
>>> just for the record).
>>>
>>> In GRASS GIS 7.0.4 and 7.1.svn, g.proj -w gave me the projection
>>
>> information as in the GeoTIFF, though
>>
>>> without any AUTHORITY parameters (as it is supposed to according to the
>>
>> manual). The AUTHORITY tags seem
>>
>>> to be added by r.out.gdal, except for the authority parameter for the
>>
>> entire projection... Even if g.proj
>>
>>> -g returns epsg=25832...
>>>
>>> gdal_translate -a_srs EPSG:25832 in.tif out.tif adds the AUTHORITY
>>
>> parameter to the GeoTIFF, while
>>
>>> gdal_edit -a_srs EPSG:25832 does not.
>>
>> One route to writing from R (from a SpatialGridDataFrame or another object
>> coerced to that class) is:
>>
>> SEXP
>> RGDAL_SetProject(SEXP sxpDataset, SEXP proj4string) {
>>
>>   OGRSpatialReference oSRS;
>>   char *pszSRS_WKT = NULL;
>>
>>   GDALDataset *pDataset = getGDALDatasetPtr(sxpDataset);
>>
>>   installErrorHandler();
>>   oSRS.importFromProj4(CHAR(STRING_ELT(proj4string, 0)));
>>   oSRS.exportToWkt( &pszSRS_WKT );
>>   uninstallErrorHandlerAndTriggerError();
>>
>>   installErrorHandler();
>>   OGRErr err = pDataset->SetProjection(pszSRS_WKT);
>>   CPLFree( pszSRS_WKT );
>>
>>   if (err == CE_Failure)
>> 	warning("Failed to set projection\n");
>>   uninstallErrorHandlerAndTriggerError();
>>
>>   return(sxpDataset);
>> }
>>
>> It is as you see going from a proj4 string (used in Spatial objects to hold
>> the CRS/SRS) to WKT. Could this be re-written to provide the Authority, or
>> is Authority unrepresented in proj4?
>
> If you've an expanded proj.4 string definition like "+proj=tmerc +lat_0=0
> +lon_0=93 +k=1 +x_0=500000 +y_0=0 +ellps=krass
> +towgs84=15.8,-154.4,-82.3,0,0,0,0 +units=m +no_defs", then the link with the
> authority is lost. It would be nice to have some logic in GDAL to try to
> retrieve it by comparing with the EPSG database, but this isn't implemented
> yet (and there could be some ambiguities as you can find some EPSG codes that
> resolve to the same proj.4 string)
>
> Well, strings like "+init=epsg:21456" are also valid proj.4 strings and
> importFromProj4() will then be able to set the authority node. But that
> requires storing this un-expanded proj.4 string instead of its resolved
> definition.

Even, thanks!

This may be the resolution, where we change the definition of the R CRS 
class to add a slot to hold the un-expanded +init=epsg:* if available. I 
agree that guessing EPSG codes from proj4 strings is often unreliable, 
I've tried it. @Edzer - could you consider adding this to:

https://github.com/edzer/sfr

I've put this thread as an issue.

Best wishes,

Roger

>
> Actually proj.4 is also happy if the string contains both +init and +proj :
> "+init=epsg:21456 +proj=tmerc +lat_0=0 +lon_0=93 +k=1 +x_0=500000 +y_0=0
> +ellps=krass +towgs84=15.8,-154.4,-82.3,0,0,0,0 +units=m +no_defs". But you
> have to be careful that the +init and the rest of the definition are
> consistant, since I'd say it is pretty much undefined which part will win over
> this other one.
>
> Even
>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
http://depsy.org/person/434412