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

[Even Rouault]

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.

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

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com