[Geotiff] geos projection sweep parameter

Mon Jul 10 10:10:48 EST 2017

Even,

So to summarize, a typical geotiff reading client should look for the 
standard geotiff projection information, if not found look for "ESRI PE 
String". Right? You made it sound like ArcGIS uses its own form of WKT? 
Does it also use the "ESRI PE String" key?

As you stated the two "possible" solutions are to get the "geos" 
projection (and others) added to the geotiff specification or to get 
"sweep" added to the ESRI WKT specification. Both of which would require 
client libraries to be updated to properly parse and use this 
information which could take years to reach a majority of the GIS 
community; maybe less given the right documentation and google search 
results.

Perhaps I will start with contacting some of the working group and see 
what path forward they might recommend. Thank you for your help.

Dave

On 7/6/17 7:18 PM, Even Rouault wrote:
> On jeudi 6 juillet 2017 17:12:44 CEST David Hoese wrote:
> 
>  > Hi,
> 
>  >
> 
>  > I work with data from the GOES-16 satellite ABI instrument. The official
> 
>  > ABI Level 1B product consists of data mapped to the GEOS projection with
> 
>  > a sweep angle axis parameter set to 'x' instead of the default 'y'. More
> 
>  > information on this projection can be found here:
> 
>  > http://proj4.org/projections/geos.html
> 
>  >
> 
>  > As of GDAL v2.1+ the "sweep" parameter is handled internally and I can't
> 
>  > remember what version of PROJ.4 added support for this. As far as I can
> 
>  > tell there is no standard way to store the "sweep" parameter in a
> 
>  > GeoTIFF file. I'm not sure if this is due to a limitation of the WKT
> 
>  > specification or if it is something that only requires changes in the
> 
>  > GeoTIFF library. Does anyone know what needs to happen to get "sweep"
> 
>  > added to a GeoTIFF? If this requires changes to libgeotiff I am willing
> 
>  > to add them or provide patches to add the functionality. Thanks for any
> 
>  > information you can provide.
> 
>  >
> 
> David,
> 
> The issue is more fundamental than the sweep parameter being stripped off.
> 
> The main issue is that the GEOS projection does not exist at all in 
> GeoTIFF representation
> 
> (GeoTIFF encoding is not nominally based on WKT, but consists of keys 
> with codified values)
> 
> The list of GeoTIFFprojection methods is :
> 
> http://web.archive.org/web/20160407200550/http://www.remotesensing.org:80/geotiff/spec/geotiff6.html#6.3.3.3
> 
> What makes things seem to work (except the sweep parameter I mean) is 
> that GDAL
> 
> embeds a form of WKT in the PCSCitationGeoKey as a fallback when there 
> is no proper
> 
> GeoTIFF way of representing the SRS
> 
> Given a source netCDF file with :
> 
> PROJCS["unnamed",
> 
> GEOGCS["WGS 84",
> 
> DATUM["WGS_1984",
> 
> SPHEROID["WGS 84",6378137,298.257223563,
> 
> AUTHORITY["EPSG","7030"]],
> 
> AUTHORITY["EPSG","6326"]],
> 
> PRIMEM["Greenwich",0,
> 
> AUTHORITY["EPSG","8901"]],
> 
> UNIT["degree",0.0174532925199433,
> 
> AUTHORITY["EPSG","9122"]],
> 
> AUTHORITY["EPSG","4326"]],
> 
> PROJECTION["Geostationary_Satellite"],
> 
> PARAMETER["central_meridian",-89.5],
> 
> PARAMETER["satellite_height",35785831],
> 
> PARAMETER["false_easting",0],
> 
> PARAMETER["false_northing",0],
> 
> EXTENSION["PROJ4","+proj=geos +lon_0=-89.5 +h=35785831 +x_0=0 +y_0=0 
> +datum=WGS84 +units=m +no_defs +sweep=x"]]
> 
> Origin = (-5434865.725631997920573,-721442.352958934963681)
> 
> Pixel Size = (6012.019607998958236,6012.019607998958236)
> 
> gdal_translat'ing to TIF gives a GeoTIFF whose listgeo dump is
> 
> Geotiff_Information:
> 
> Version: 1
> 
> Key_Revision: 1.0
> 
> Tagged_Information:
> 
> ModelTransformationTag (4,4):
> 
> 6012.01960799896 0 0 -5434865.725632
> 
> 0 6012.01960799896 0 -721442.352958935
> 
> 0 0 0 0
> 
> 0 0 0 1
> 
> End_Of_Tags.
> 
> Keyed_Information:
> 
> GTModelTypeGeoKey (Short,1): User-Defined
> 
> GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
> 
> GTCitationGeoKey (Ascii,24): "Geostationary_Satellite"
> 
> GeographicTypeGeoKey (Short,1): GCS_WGS_84
> 
> GeogCitationGeoKey (Ascii,13): "GCS_WGS_1984"
> 
> GeogAngularUnitsGeoKey (Short,1): Angular_Degree
> 
> GeogSemiMajorAxisGeoKey (Double,1): 6378137
> 
> GeogInvFlatteningGeoKey (Double,1): 298.257223563
> 
> PCSCitationGeoKey (Ascii,383): "ESRI PE String = 
> PROJCS["Geostationary_Satellite",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Geostationary_Satellite"],PARAMETER["central_meridian",-89.5],PARAMETER["satellite_height",35785831],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1]]"
> 
> ProjLinearUnitsGeoKey (Short,1): Linear_Meter
> 
> End_Of_Keys.
> 
> End_Of_Geotiff.
> 
> The WKT form used is the ESRI variant of WKT. This is the way GDAL tries 
> to fit a SRS in GeoTIFF when
> 
> there is no proper GeoTIFF encoding for it (for supported GeoTIFF SRS, 
> this ESRI WKT
> 
> is not set). This is what ArcGIS is supposed to use too (not
> 
> completely sure if the massaging of OGC WKT to ESRI WKT completely 
> matches ArcGIS though,
> 
> especially for geos which is quite new).
> 
> The issue is that ESRI WKT doesn't support the EXTENSION PROJ4 node, and 
> hence the sweep
> 
> parameter is lost.
> 
> If there is an official way in ESRI WKT to represent the sweep 
> parameter, then we could use it in GDAL.
> 
> So actions could be to try investigating with ESRI how they would 
> represent that. Or as a fallback
> 
> preserve the EXTENSION node in that case, hoping that this doesn't cause 
> too much issues to
> 
> other software (but given that the current generated ESRI WKT already 
> lacks an important information,
> 
> interoperability is already non existent...)
> 
> A proper solution would be of course to map GEOS to GeoTIFF, but 
> unfortunately there has not
> 
> been activity in years to make official evolutions to the specification 
> (OGC created a charter for a GeoTIFF
> 
> SWG , http://www.opengeospatial.org/projects/groups/geotiffswg , but we 
> haven't seen any
> 
> public activity coming out from that)
> 
> For the sake of completness, you can workaround the issue by creating a 
> GDAL PAM .aux.xml
> 
> sidecar file
> 
> gdal_translate in.nc out.tif
> 
> gdal_translate in.nc out.png -of PNG (just to get a out.png.aux.xml with 
> the GDAL WKT)
> 
> mv out.png.aux.xml out.tif.aux.xml
> 
> gdalinfo out.tif
> 
> Even
> 
> -- 
> 
> Spatialys - Geospatial professional services
> 
> http://www.spatialys.com
> 