Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Le mercredi 15 avril 2015 22:48:45, Even Rouault a écrit : The way I see it, there are two different ways to model Google Mercator: 1. As a WGS84 datum/ellipsoid with a custom Mercator-based projection which uses only the semi-major axis, as seems to be currently done in EPSG: http://epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857 2. As the Google spheroid with Mercator 1SP and no transformation when going to WGS 84. I don't think #1 is currently an option given the list of coordinate transformation codes available in the GeoTIFF specification. Right, there's no such coordinate transformation for GeoTIFF Thinking a bit more about this issue, a dedicated coordinate transformation would be the cleanest way, both in GeoTIFF and proj.4 (or a flag in proj.4 to indicate that the spherical version of the projection method, if it exists and that's the case for mercator, must be used even if the ellipsoid definition is not spherical) The fact that we manage currently to do WebMercator with proj.4 is due to WebMercator being based on WGS84 datum, and WGS84 being the pivot datum used by proj.4 when doing datum transformations. The +nadgrids=@null basically means that the transform from the datum defined with ellipsoid parameters +a=6378137 +b=6378137 and grid @null to datum WGS84 is the identity (for that part the values of +a and +b are completely ignored. They are only used when going from (long,lat) to (x,y)). If we ever wanted to do spherical mercator on another datum than WGS84, there's currently no way to express that with proj.4 (except through a grid from the datum to WGS84 perhaps) All that said, that doesn't really help me solving my issue, but makes me believe that with what currently exists in the GeoTIFF spec, there's no clean to express WebMercator. - ProjectedCSTypeGeoKey = 3857 is probably the best one, although it is an extension of the original GeoTIFF spec, and that isn't understood properly by some readers - All formulations that try to expand the definition with ProjCoordTransGeoKey = CT_Mercator, its projection parameters and GCS parameters aren't really appropriate, since there's no way of capturing that its a spherical mercator that must be used. For #2, I would think/hope that if you specified the identity 7-parameter transformation in the GeoTIFF header, ArcGIS would properly reproject the image: Hum, this will not work I'm afraid. See my response to Dmitry's to see why +towgs84=0,0,0,0,0,0,0 is not appropriate. (Plus the fact that the +towgs84 in geotiff is probably not very widely supported.) http://trac.osgeo.org/geotiff/ticket/26 André ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Hi Pepijn, How does this relate to http://trac.osgeo.org/geotiff/wiki/RefiningGeoTIFF? http://trac.osgeo.org/geotiff/wiki/RefiningGeoTIFF? Very much (although I'm foolishly attempting to reach interoperability with existing implementations) Frank wrote there that accepted industry practice has been to accept newer EPSG PCS and GCS codes even though they are not explicitly listed in the GeoTIFF specification It’s not clear to me from the GeoTIFF spec what the correct way to convert EPSG PCS codes to GeoTIFF PCS codes is, but it does seem to indicate that this is done somehow already. The relevant section is : http://www.remotesensing.org/geotiff/spec/geotiff6.html#6.3.3 6.3.3.1 Projected CS Type Codes Ranges: [1, 1000] = Obsolete EPSG/POSC Projection System Codes [2, 32760] = EPSG Projection System codes 32767= user-defined [32768, 65535] = Private User Implementations But there's a big hole for ]1000,2[, and 3857 is in it. And even if it was in the [2000,3276] range, it didn't exist at the time of the latest GeoTIFF spec. Frank's accepted industry practice is to just make Projected CS Type Code = EPSG PCS code. OGC started a SWG to do this GeoTIFF spec refinement work. Might be worthwhile to ask that group for input, although looking at their mailing list, github repo and wiki pages there doesn’t seem to be much going on at the moment (or at least there’s no record of it). Yes, I've contacted one of the members. I can't see a geotiff list at https://lists.opengeospatial.org/mailman/listinfo . Is it a public one ? Thanks, Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Hi Even, How does this relate to http://trac.osgeo.org/geotiff/wiki/RefiningGeoTIFF? http://trac.osgeo.org/geotiff/wiki/RefiningGeoTIFF? Frank wrote there that accepted industry practice has been to accept newer EPSG PCS and GCS codes even though they are not explicitly listed in the GeoTIFF specification It’s not clear to me from the GeoTIFF spec what the correct way to convert EPSG PCS codes to GeoTIFF PCS codes is, but it does seem to indicate that this is done somehow already. OGC started a SWG to do this GeoTIFF spec refinement work. Might be worthwhile to ask that group for input, although looking at their mailing list, github repo and wiki pages there doesn’t seem to be much going on at the moment (or at least there’s no record of it). Best regards, Pepijn___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
- All formulations that try to expand the definition with ProjCoordTransGeoKey = CT_Mercator, its projection parameters and GCS parameters aren't really appropriate, since there's no way of capturing that its a spherical mercator that must be used. Well, I agree that most readers probably won't be able to handle it, that being said, if the de facto convention is to fall-back to EPSG for codes that are not understood, setting that key might be your best option. Setting the ProjCoordTransGeoKey to 1024 as well might be a good idea (with WGS 84 as the ellipsoid). The EPSG operation methods, which currently fall in the [1000-1] range, are in the range of the supported transformation codes: 6.3.3.3 Coordinate Transformation Codes Ranges: 0 = undefined [1, 16383] = GeoTIFF Coordinate Transformation codes [16384, 32766] = Reserved by GeoTIFF 32767 = user-defined [32768, 65535] = Private User Implementations That is, fallback to an EPSG code for the projection if the projection does not map to one that is defined in the specification. Of course, current readers probably will not handle it correctly. Hum, the issue is that the values currently defined by GeoTIFF in the [1, 16383] range has nothing to do with the EPSG conversion methods. Your proposal might be reasonable, but it is AFAIK not at all a current industry practice. This would also imply definiing if it is allowed to use the EPSG conversion method when there are corresponding existing GeoTIFF enumerated values Then again, I cannot see how any readers other than the ones who lookup the ProjectedCSTypeGeoKey in the EPSG registry could realistically work correctly at the moment. True. That's more or less what GDAL does currently. In fact, it is more complicated than that. GDAL relies on the libgeotiff function GTIFGetDefn() which, from the raw GeoTIFF keys and the EPSG dictionnaries as .csv files, tries to build a fully expanded SRS definition using GeoTIFF enumerations. Which explains why GDAL can write some PCS codes but not read a correct SRS back, because there might be no GeoTIFF projection method corresponding. GDAL could (should) likely just use its importFromEPSG() method in those cases. Not sure why we don't do that already. -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
On 2015-04-16 12:07, Even Rouault wrote: - All formulations that try to expand the definition with ProjCoordTransGeoKey = CT_Mercator, its projection parameters and GCS parameters aren't really appropriate, since there's no way of capturing that its a spherical mercator that must be used. Well, I agree that most readers probably won't be able to handle it, that being said, if the de facto convention is to fall-back to EPSG for codes that are not understood, setting that key might be your best option. Setting the ProjCoordTransGeoKey to 1024 as well might be a good idea (with WGS 84 as the ellipsoid). The EPSG operation methods, which currently fall in the [1000-1] range, are in the range of the supported transformation codes: 6.3.3.3 Coordinate Transformation Codes Ranges: 0 = undefined [1, 16383] = GeoTIFF Coordinate Transformation codes [16384, 32766] = Reserved by GeoTIFF 32767 = user-defined [32768, 65535] = Private User Implementations That is, fallback to an EPSG code for the projection if the projection does not map to one that is defined in the specification. Of course, current readers probably will not handle it correctly. Hum, the issue is that the values currently defined by GeoTIFF in the [1, 16383] range has nothing to do with the EPSG conversion methods. Your proposal might be reasonable, but it is AFAIK not at all a current industry practice. This would also imply definiing if it is allowed to use the EPSG conversion method when there are corresponding existing GeoTIFF enumerated values Yeah, I was probably overreaching with that. Probably better to just make sure the next version of the specification has a transformation code available for Google Mercator as a coordinate transformation. I just don't like having a file with information that will not transform correctly in the event the client does not lookup that projected CRS code in EPSG. Having the spheroid as the datum with the Mercator 1SP projection leads to way too many special cases. I'd rather have a projection that some readers might not understand. Seems safer and more explicit. Then again, I cannot see how any readers other than the ones who lookup the ProjectedCSTypeGeoKey in the EPSG registry could realistically work correctly at the moment. True. That's more or less what GDAL does currently. In fact, it is more complicated than that. GDAL relies on the libgeotiff function GTIFGetDefn() which, from the raw GeoTIFF keys and the EPSG dictionnaries as .csv files, tries to build a fully expanded SRS definition using GeoTIFF enumerations. Which explains why GDAL can write some PCS codes but not read a correct SRS back, because there might be no GeoTIFF projection method corresponding. GDAL could (should) likely just use its importFromEPSG() method in those cases. Not sure why we don't do that already. Trying to import from EPSG first makes a lot of sense to me. Cheers, André ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Ranges: [1, 1000] = Obsolete EPSG/POSC Projection System Codes [2, 32760] = EPSG Projection System codes 32767= user-defined [32768, 65535] = Private User Implementations But there's a big hole for ]1000,2[, and 3857 is in it. And even if it was in the [2000,3276] range, it didn't exist at the time of the latest GeoTIFF spec. Frank's accepted industry practice is to just make Projected CS Type Code = EPSG PCS code. Ok. I was looking for some rules that translate the [2, 32760] range to EPSG codes. OGC started a SWG to do this GeoTIFF spec refinement work. Might be worthwhile to ask that group for input, although looking at their mailing list, github repo and wiki pages there doesn’t seem to be much going on at the moment (or at least there’s no record of it). Yes, I've contacted one of the members. I can't see a geotiff list at https://lists.opengeospatial.org/mailman/listinfo . Is it a public one ? The archives are private so you probably need to be an OGC member to access them. I’m not working at Luciad anymore so I can’t be of much help myself; I haven’t registered as an individual member. My old Luciad account is still active which gives me access to the archives, but I’m not sure how long that will remain the case. Best regards, Pepijn ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
On 2015-04-16 07:08, Even Rouault wrote: Le mercredi 15 avril 2015 22:48:45, Even Rouault a écrit : The way I see it, there are two different ways to model Google Mercator: 1. As a WGS84 datum/ellipsoid with a custom Mercator-based projection which uses only the semi-major axis, as seems to be currently done in EPSG: http://epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857 2. As the Google spheroid with Mercator 1SP and no transformation when going to WGS 84. I don't think #1 is currently an option given the list of coordinate transformation codes available in the GeoTIFF specification. Right, there's no such coordinate transformation for GeoTIFF Thinking a bit more about this issue, a dedicated coordinate transformation would be the cleanest way, both in GeoTIFF and proj.4 (or a flag in proj.4 to indicate that the spherical version of the projection method, if it exists and that's the case for mercator, must be used even if the ellipsoid definition is not spherical) The fact that we manage currently to do WebMercator with proj.4 is due to WebMercator being based on WGS84 datum, and WGS84 being the pivot datum used by proj.4 when doing datum transformations. The +nadgrids=@null basically means that the transform from the datum defined with ellipsoid parameters +a=6378137 +b=6378137 and grid @null to datum WGS84 is the identity (for that part the values of +a and +b are completely ignored. They are only used when going from (long,lat) to (x,y)). If we ever wanted to do spherical mercator on another datum than WGS84, there's currently no way to express that with proj.4 (except through a grid from the datum to WGS84 perhaps) All that said, that doesn't really help me solving my issue, but makes me believe that with what currently exists in the GeoTIFF spec, there's no clean to express WebMercator. - ProjectedCSTypeGeoKey = 3857 is probably the best one, although it is an extension of the original GeoTIFF spec, and that isn't understood properly by some readers - All formulations that try to expand the definition with ProjCoordTransGeoKey = CT_Mercator, its projection parameters and GCS parameters aren't really appropriate, since there's no way of capturing that its a spherical mercator that must be used. Well, I agree that most readers probably won't be able to handle it, that being said, if the de facto convention is to fall-back to EPSG for codes that are not understood, setting that key might be your best option. Setting the ProjCoordTransGeoKey to 1024 as well might be a good idea (with WGS 84 as the ellipsoid). The EPSG operation methods, which currently fall in the [1000-1] range, are in the range of the supported transformation codes: 6.3.3.3 Coordinate Transformation Codes Ranges: 0 = undefined [1, 16383] = GeoTIFF Coordinate Transformation codes [16384, 32766] = Reserved by GeoTIFF 32767 = user-defined [32768, 65535] = Private User Implementations That is, fallback to an EPSG code for the projection if the projection does not map to one that is defined in the specification. Of course, current readers probably will not handle it correctly. Then again, I cannot see how any readers other than the ones who lookup the ProjectedCSTypeGeoKey in the EPSG registry could realistically work correctly at the moment. Based on the outcome of http://lists.osgeo.org/pipermail/gdal-dev/2015-April/041532.html, WKT readers should be able to handle a WKT string with a projection that is named Popular Visualisation Pseudo Mercator if they do in fact support it, so the WKT string for 3857 could also be updated to be modeled the same way. For #2, I would think/hope that if you specified the identity 7-parameter transformation in the GeoTIFF header, ArcGIS would properly reproject the image: Hum, this will not work I'm afraid. See my response to Dmitry's to see why +towgs84=0,0,0,0,0,0,0 is not appropriate. (Plus the fact that the +towgs84 in geotiff is probably not very widely supported.) Right, from some reason, I my mind towgs84=0,0,0,0,0,0,0 was the identity transformation yesterday as I wrote my reply, but obviously it still would need to go to/from geocentric which would cause problems. André ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Dis you try with 102113 or 102100 instead of 3857 ? I had once to use those code in order to be discovered as 3857 in arcgis. when i created geotiff with gdal warp/transform. Le 15 avr. 2015 18:16, Even Rouault even.roua...@spatialys.com a écrit : Hi, I've collected in http://trac.osgeo.org/gdal/ticket/5924 a summary of issues and findings when trying to write GeoTIFF files in EPSG:3857 in a standard way (using ProjectedCSTypeGeoKey = 3857), while making them correctly displayed by ArcGIS (and ideally by other independant implementations). The current situation is that there's no such way (AFAIK). I'd appreciate any feedback/suggestion related to that issue. Best regards, Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Hi, I've collected in http://trac.osgeo.org/gdal/ticket/5924 a summary of issues and findings when trying to write GeoTIFF files in EPSG:3857 in a standard way (using ProjectedCSTypeGeoKey = 3857), while making them correctly displayed by ArcGIS (and ideally by other independant implementations). The current situation is that there's no such way (AFAIK). I'd appreciate any feedback/suggestion related to that issue. Best regards, Even -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Le mercredi 15 avril 2015 22:14:37, Dmitry Baryshnikov a écrit : As for me the main problem is, that gdal write wrong WKT for 3857 (http://trac.osgeo.org/gdal/ticket/3962). The Geographic SRS use WGS84 ellipsoid (SPHEROID[WGS_1984,6378137,298.257223563]]) instead sphere (EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs],) and all the data shifts. I do think that the above SPHEROID WKT and Proj.4 are both OK, even if they don't look consistent. The datum is technically WGS84 datum (ellipsoidal), even if the projection only uses the semi major axis (spherical). This is all the hell of this WebMercator projection ! We use a hack of proj.4 to use classical mercator transformation with spherical parameters, while still being considered as a WGS84 datum. When reprojecting from a non-WGS 84 datum (let's say EPSG:4322) to EPSG:3857, the datum transformation should be from this non-WGS 84 to ellipsoid WGS 84, and then the WebMercator projection uses the spherical parameters to compute the easting/northing. Look: 1) Normal EPSG:3857 transformation (with +nadgrids=@null) : $ echo 2 49 | gdaltransform -s_srs EPSG:4322 -t_srs EPSG:3857 222656.112419298 6274866.20215882 2.9538588412106 -- correct. Same result as operating in 2 steps : $ echo 2 49 | gdaltransform -s_srs EPSG:4322 -t_srs EPSG:4326 | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857 222656.112419298 6274866.20215883 2.9538588412106 2) Incorrect definition using +towgs84=0,0,0,0,0,0,0 : $ gdaltransform -s_srs EPSG:4322 -t_srs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +towgs84=0,0,0,0,0,0,0 2 49 222656.112419297 6242580.32725437 -12133.4419494262 2 49 is transformed from EPSG:4322 datum to a spherical datum, instead to the ellipsoid WGS84 datum. 3) Incorrect definition without +towgs84 or +nadgrids=@null : no datum transformation at all is done since its only an ellipsoid definition, and proj.4 will not do any ellipsoid transformation in that case $ gdaltransform -s_srs EPSG:4322 -t_srs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m 2 49 222638.981586547 6274861.39400658 0 Best regards, Dmitry 15.04.2015 19:15, Even Rouault пишет: Hi, I've collected in http://trac.osgeo.org/gdal/ticket/5924 a summary of issues and findings when trying to write GeoTIFF files in EPSG:3857 in a standard way (using ProjectedCSTypeGeoKey = 3857), while making them correctly displayed by ArcGIS (and ideally by other independant implementations). The current situation is that there's no such way (AFAIK). I'd appreciate any feedback/suggestion related to that issue. Best regards, Even ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
The way I see it, there are two different ways to model Google Mercator: 1. As a WGS84 datum/ellipsoid with a custom Mercator-based projection which uses only the semi-major axis, as seems to be currently done in EPSG: http://epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857 2. As the Google spheroid with Mercator 1SP and no transformation when going to WGS 84. I don't think #1 is currently an option given the list of coordinate transformation codes available in the GeoTIFF specification. Right, there's no such coordinate transformation for GeoTIFF For #2, I would think/hope that if you specified the identity 7-parameter transformation in the GeoTIFF header, ArcGIS would properly reproject the image: Hum, this will not work I'm afraid. See my response to Dmitry's to see why +towgs84=0,0,0,0,0,0,0 is not appropriate. (Plus the fact that the +towgs84 in geotiff is probably not very widely supported.) http://trac.osgeo.org/geotiff/ticket/26 André ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
As for me the main problem is, that gdal write wrong WKT for 3857 (http://trac.osgeo.org/gdal/ticket/3962). The Geographic SRS use WGS84 ellipsoid (SPHEROID[WGS_1984,6378137,298.257223563]]) instead sphere (EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs],) and all the data shifts. Best regards, Dmitry 15.04.2015 19:15, Even Rouault пишет: Hi, I've collected in http://trac.osgeo.org/gdal/ticket/5924 a summary of issues and findings when trying to write GeoTIFF files in EPSG:3857 in a standard way (using ProjectedCSTypeGeoKey = 3857), while making them correctly displayed by ArcGIS (and ideally by other independant implementations). The current situation is that there's no such way (AFAIK). I'd appreciate any feedback/suggestion related to that issue. Best regards, Even ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
On 2015-04-15 14:59, Even Rouault wrote:
Le mercredi 15 avril 2015 18:21:51, xavier lhomme a écrit :
Dis you try with 102113 or 102100 instead of 3857 ? I had once to use
those code in order to be discovered as 3857 in arcgis. when i created
geotiff with gdal warp/transform.
Xavier,
thanks for the hint. I've just tried it and thinks it qualifies as yet another
not completely satisfactory solution.
Using -a_srs EPSG:102113 effectively results in a correctly placed geotiff in
ArcGIS (with a warning about it being in GCS_WGS_1984_Major_Auxiliary_Sphere),
but there are several points that are not that great :
- this is not a standard EPSG code
- due to being 32767, it is not encoded in ProjectedCSTypeGeoKey, which was
one of my requirement
- when read by GDAL, the expanded GeoTIFF definition isn't properly identified
as WebMercator, consequently the proj.4 string lacks the +nadgrids=@null
stuff, which makes it a non-complete datum definition and cause issues when
reprojecting such GeoTIFF with GDAL (lack of datum transformation from/to the
WGS84 datum of WebMercator). That could admitedly being tweaked in GDAL.
For the record, the GeoTIFF keys written by GDAL are:
{{{
GTModelTypeGeoKey (Short,1): ModelTypeProjected
GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
GTCitationGeoKey (Ascii,22): WGS_1984_Web_Mercator
GeographicTypeGeoKey (Short,1): User-Defined
GeogCitationGeoKey (Ascii,151): GCS Name =
GCS_WGS_1984_Major_Auxiliary_Sphere|Datum = WGS_1984_Major_Auxiliary_Sphere|
Ellipsoid = WGS_1984_Major_Auxiliary_Sphere|Primem = Greenwich|
GeogGeodeticDatumGeoKey (Short,1): User-Defined
GeogAngularUnitsGeoKey (Short,1): Angular_Degree
GeogEllipsoidGeoKey (Short,1): User-Defined
GeogSemiMajorAxisGeoKey (Double,1): 6378137
GeogSemiMinorAxisGeoKey (Double,1): 6378137
GeogPrimeMeridianLongGeoKey (Double,1): 0
ProjectedCSTypeGeoKey (Short,1): User-Defined
ProjectionGeoKey (Short,1): User-Defined
ProjCoordTransGeoKey (Short,1): CT_Mercator
ProjLinearUnitsGeoKey (Short,1): Linear_Meter
ProjNatOriginLongGeoKey (Double,1): 0
ProjNatOriginLatGeoKey (Double,1): 0
ProjFalseEastingGeoKey (Double,1): 0
ProjFalseNorthingGeoKey (Double,1): 0
ProjScaleAtNatOriginGeoKey (Double,1): 1
}}}
And GDAL reads that as :
{{{
PROJCS[WGS_1984_Web_Mercator,
GEOGCS[GCS_WGS_1984_Major_Auxiliary_Sphere,
DATUM[WGS_1984_Major_Auxiliary_Sphere,
SPHEROID[WGS_1984_Major_Auxiliary_Sphere,6378137,0]],
PRIMEM[Greenwich,0],
UNIT[degree,0.0174532925199433]],
PROJECTION[Mercator_1SP],
PARAMETER[central_meridian,0],
PARAMETER[scale_factor,1],
PARAMETER[false_easting,0],
PARAMETER[false_northing,0],
UNIT[metre,1,
AUTHORITY[EPSG,9001]]]
}}}
(not sure why it lacks the latitude_of_origin, but that's another minor
issue)
ArcCatalog identifies the SRS as:
{{{
WGS_1984_Web_Mercator
WKID: 3785 Authority: EPSG
Projection: Mercator
false_easting: 0.0
false_northing: 0.0
central_meridian: 0.0
standard_parallel_1: 0.0
Linear Unit: Meter (1.0)
Geographic Coordinate System: GCS_WGS_1984_Major_Auxiliary_Sphere
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984_Major_Auxiliary_Sphere
Spheroid: WGS_1984_Major_Auxiliary_Sphere
Semimajor Axis: 6378137.0
Semiminor Axis: 6378137.0
Inverse Flattening: 0.0
}}}
Even
The way I see it, there are two different ways to model Google Mercator:
1. As a WGS84 datum/ellipsoid with a custom Mercator-based projection
which uses only the semi-major axis, as seems to be currently done in
EPSG: http://epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857
2. As the Google spheroid with Mercator 1SP and no transformation when
going to WGS 84.
I don't think #1 is currently an option given the list of coordinate
transformation codes available in the GeoTIFF specification.
For #2, I would think/hope that if you specified the identity
7-parameter transformation in the GeoTIFF header, ArcGIS would properly
reproject the image:
http://trac.osgeo.org/geotiff/ticket/26
André
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
15.04.2015 23:46, Even Rouault пишет: Le mercredi 15 avril 2015 22:14:37, Dmitry Baryshnikov a écrit : As for me the main problem is, that gdal write wrong WKT for 3857 (http://trac.osgeo.org/gdal/ticket/3962). The Geographic SRS use WGS84 ellipsoid (SPHEROID[WGS_1984,6378137,298.257223563]]) instead sphere (EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs],) and all the data shifts. I do think that the above SPHEROID WKT and Proj.4 are both OK, even if they don't look consistent. The ArcGIS (or QGIS in qpj) save 3857 in such prj file: PROJCS[WGS 84 / Pseudo-Mercator, 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[Mercator_1SP], PARAMETER[central_meridian,0], PARAMETER[scale_factor,1], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AXIS[X,EAST], AXIS[Y,NORTH], EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs], AUTHORITY[EPSG,3857]] And GDAL opens such shape file and set the spheroid. But if I save the shape file in 3857 via GDAL I get ellipsoid: PROJCS[WGS_84_Pseudo_Mercator, GEOGCS[GCS_WGS_1984, DATUM[D_WGS_1984, SPHEROID[WGS_1984,6378137,298.257223563]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]], PROJECTION[Mercator], PARAMETER[central_meridian,0], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[Meter,1], PARAMETER[standard_parallel_1,0.0]] And both GDAL(QGIS) and ArcGIS open this shape file and set the ellipsoid! The datum is technically WGS84 datum (ellipsoidal), even if the projection only uses the semi major axis (spherical). This is all the hell of this WebMercator projection ! We use a hack of proj.4 to use classical mercator transformation with spherical parameters, while still being considered as a WGS84 datum. When reprojecting from a non-WGS 84 datum (let's say EPSG:4322) to EPSG:3857, the datum transformation should be from this non-WGS 84 to ellipsoid WGS 84, and then the WebMercator projection uses the spherical parameters to compute the easting/northing. Look: 1) Normal EPSG:3857 transformation (with +nadgrids=@null) : $ echo 2 49 | gdaltransform -s_srs EPSG:4322 -t_srs EPSG:3857 222656.112419298 6274866.20215882 2.9538588412106 -- correct. Same result as operating in 2 steps : $ echo 2 49 | gdaltransform -s_srs EPSG:4322 -t_srs EPSG:4326 | gdaltransform -s_srs EPSG:4326 -t_srs EPSG:3857 222656.112419298 6274866.20215883 2.9538588412106 2) Incorrect definition using +towgs84=0,0,0,0,0,0,0 : $ gdaltransform -s_srs EPSG:4322 -t_srs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +towgs84=0,0,0,0,0,0,0 2 49 222656.112419297 6242580.32725437 -12133.4419494262 2 49 is transformed from EPSG:4322 datum to a spherical datum, instead to the ellipsoid WGS84 datum. 3) Incorrect definition without +towgs84 or +nadgrids=@null : no datum transformation at all is done since its only an ellipsoid definition, and proj.4 will not do any ellipsoid transformation in that case $ gdaltransform -s_srs EPSG:4322 -t_srs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m 2 49 222638.981586547 6274861.39400658 0 Best regards, Dmitry 15.04.2015 19:15, Even Rouault пишет: Hi, I've collected in http://trac.osgeo.org/gdal/ticket/5924 a summary of issues and findings when trying to write GeoTIFF files in EPSG:3857 in a standard way (using ProjectedCSTypeGeoKey = 3857), while making them correctly displayed by ArcGIS (and ideally by other independant implementations). The current situation is that there's no such way (AFAIK). I'd appreciate any feedback/suggestion related to that issue. Best regards, Even ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev Best regards, Dmitry ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
The ArcGIS (or QGIS in qpj) save 3857 in such prj file: I hoped that this thread could remain focused just about GeoTIFF encoding. Shapefile encoding issues are a different matter (see http://trac.osgeo.org/gdal/ticket/3962) PROJCS[WGS 84 / Pseudo-Mercator, 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[Mercator_1SP], PARAMETER[central_meridian,0], PARAMETER[scale_factor,1], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AXIS[X,EAST], AXIS[Y,NORTH], EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs], AUTHORITY[EPSG,3857]] Are you 100% positive that ArcGIS generates such a .prj ? This doesn't look like classical ESRI WKT (which generally lacks most AUTHORITY nodes and has no EXTENSION PROJ4 stuff), but like GDAL WKT. And GDAL opens such shape file and set the spheroid. But if I save the shape file in 3857 via GDAL I get ellipsoid: PROJCS[WGS_84_Pseudo_Mercator, GEOGCS[GCS_WGS_1984, DATUM[D_WGS_1984, SPHEROID[WGS_1984,6378137,298.257223563]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]], PROJECTION[Mercator], PARAMETER[central_meridian,0], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[Meter,1], PARAMETER[standard_parallel_1,0.0]] Yes, that's the result after morphToESRI(). Which is probably not the right ESRI WKT for WebMercator, anyway... Should rather be something like the following according to http://trac.osgeo.org/gdal/ticket/3962 : PROJCS[WGS_1984_Web_Mercator_Auxiliary_Sphere, GEOGCS[GCS_WGS_1984, DATUM[D_WGS_1984, SPHEROID[WGS_1984,6378137.0,298.257223563]], PRIMEM[Greenwich,0.0], UNIT[Degree,0.0174532925199433]], PROJECTION[Mercator_Auxiliary_Sphere], PARAMETER[False_Easting,0.0], PARAMETER[False_Northing,0.0], PARAMETER[Central_Meridian,0.0], PARAMETER[Standard_Parallel_1,0.0], PARAMETER[Auxiliary_Sphere_Type,0.0], UNIT[Meter,1.0], AUTHORITY[ESRI,102100]] And both GDAL(QGIS) and ArcGIS open this shape file and set the ellipsoid! I don't understand what you mean here. -- Spatialys - Geospatial professional services http://www.spatialys.com ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Encoding EPSG:3857 (WebMercator) in GeoTIFF, and ArcGIS interoperability
Sorry for some spam in GeoTIFF topic. I'll create some shape files and screenshots and attach them to the ticket #3962 to explain what I mean. Best regards, Dmitry 16.04.2015 00:10, Even Rouault пишет: The ArcGIS (or QGIS in qpj) save 3857 in such prj file: I hoped that this thread could remain focused just about GeoTIFF encoding. Shapefile encoding issues are a different matter (see http://trac.osgeo.org/gdal/ticket/3962) PROJCS[WGS 84 / Pseudo-Mercator, 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[Mercator_1SP], PARAMETER[central_meridian,0], PARAMETER[scale_factor,1], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]], AXIS[X,EAST], AXIS[Y,NORTH], EXTENSION[PROJ4,+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs], AUTHORITY[EPSG,3857]] Are you 100% positive that ArcGIS generates such a .prj ? This doesn't look like classical ESRI WKT (which generally lacks most AUTHORITY nodes and has no EXTENSION PROJ4 stuff), but like GDAL WKT. And GDAL opens such shape file and set the spheroid. But if I save the shape file in 3857 via GDAL I get ellipsoid: PROJCS[WGS_84_Pseudo_Mercator, GEOGCS[GCS_WGS_1984, DATUM[D_WGS_1984, SPHEROID[WGS_1984,6378137,298.257223563]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]], PROJECTION[Mercator], PARAMETER[central_meridian,0], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[Meter,1], PARAMETER[standard_parallel_1,0.0]] Yes, that's the result after morphToESRI(). Which is probably not the right ESRI WKT for WebMercator, anyway... Should rather be something like the following according to http://trac.osgeo.org/gdal/ticket/3962 : PROJCS[WGS_1984_Web_Mercator_Auxiliary_Sphere, GEOGCS[GCS_WGS_1984, DATUM[D_WGS_1984, SPHEROID[WGS_1984,6378137.0,298.257223563]], PRIMEM[Greenwich,0.0], UNIT[Degree,0.0174532925199433]], PROJECTION[Mercator_Auxiliary_Sphere], PARAMETER[False_Easting,0.0], PARAMETER[False_Northing,0.0], PARAMETER[Central_Meridian,0.0], PARAMETER[Standard_Parallel_1,0.0], PARAMETER[Auxiliary_Sphere_Type,0.0], UNIT[Meter,1.0], AUTHORITY[ESRI,102100]] And both GDAL(QGIS) and ArcGIS open this shape file and set the ellipsoid! I don't understand what you mean here. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
