[Geotiff] geos projection sweep parameter

Discussion:

David Hoese

2017-07-06 22:12:44 UTC

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.

Dave

Even Rouault

2017-07-07 00:18:52 UTC

Permalink

Post by David Hoese
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
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

David Hoese

2017-07-10 15:10:48 UTC

Permalink

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

Post by Even Rouault

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)
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
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
Version: 1
Key_Revision: 1.0
6012.01960799896 0 0 -5434865.725632
0 6012.01960799896 0 -721442.352958935
0 0 0 0
0 0 0 1
End_Of_Tags.
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

Even Rouault

2017-07-10 15:35:18 UTC

Permalink

Post by David Hoese
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".

Yes. But ESRI PE string is vendor specific, so a strict GeoTIFF reader will generally not understand it.

Post by David Hoese
You made it sound like ArcGIS uses its own form of WKT?

WKT v1 was not strongly standardized, so various vendors have come with their own variations,
to decide how to name projection methods, parameters, etc...
Some preview :
https://web.archive.org/web/20130728081442/http://home.gdal.org:80/projects/opengis/wktproblems.html

Post by David Hoese
Does it also use the "ESRI PE String" key?

I guess so ! I don't recall the details, but I believe that support for GeoTIFF ESRI PE string in GDAL
probably originates from merging this part of the ESRI fork of GDAL to GDAL mainline.

Post by David Hoese
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.

For GDAL-based software, potentially no change on the reading side would be required.
If the ESRI PE string contains the EXTENSION PROJ4 node (which is rather hacky), GDAL will use it.
I've just tested it by manually crafting such a GeoTIFF. Works at least with GDAL 1.11 or later.

Post by David Hoese
Perhaps I will start with contacting some of the working group and see
what path forward they might recommend.

I actually forwarded your email to one of the members of the OGC GeoTIFF SWG

Even

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