Sunday 21 June 2015

coordinate system - Hammer projection in GDAL


I am working on Ubuntu 18.04 with PROJ 4.9.3 and GDAL 2.2.3 (installed from the ubuntu-gis PPA).


In a previous question on Hammer's projection in QGIS it was concluded that with older PROJ versions the inverse of this projection was not implemented and therefore not supported by GDAL. However, with the software versions I have installed, GDAL still does not support Hammer's projection.


An example with PROJ applying Hammer's projection and its inverse:



$ cat wgs.md
15 -15
$ cs2cs +init=epsg:4326 +to +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs wgs.md > hammer.md
$ cat hammer.md
1625591.45 -1668538.36 0.00
$ cs2cs +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs +to +init=epsg:4326 hammer.md
15dE 15dS 0.000

Trying the same with GDAL:


$ gdaltransform -s_srs "+init=epsg:4326" -t_srs "+proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs" wgs.md

ERROR 1: Translating source or target SRS failed:
+proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs

Why is GDAL failing to use this projection?



Answer



Caveat


The following descriptions are the result of a quick glance through the GDAL/OGR source codes. Statements may contain errors - if so, please inform. The answer may also be updated after further reviews on the GDAL/OGR source codes.


The Short Answer


GDAL/OGR depends on PROJ.4 for all its map projection needs. However, while PROJ.4 "speaks" in PROJ.4 string format, e.g., "+proj=hammer", GDAL/OGR "speaks" in PROJCS WKT instead. The PROJCS WKT format was originally introduced by ESRI, and GDAL/OGR adopted it with some modifications/extensions.


All the GDAL/OGR apps (i.e., gdaltransform, gdalwarp, gdal_translate, and etc) insist on being able to convert User's PROJ.4 string to PROJCS WKT internally. This means all the projection's parameter values must be explicitly determined, i.e., to take whatever values given, and to assign default values to whatever values not given. In the case of "+proj=hammer", GDAL/OGR does not yet have such codes and hence it wasn't able to successfully perform the conversion. (Of course, there are also situations where PROJ.4 strings just can't be converted to PROJCS WKT.)



The authors of GDAL/OGR knew the above problem, and they built in an escape via "+wktext" to handle those situations. When GDAL/OGR apps detect this escape, they convert it to an escaped PROJCS WKT. You can try this by typing


gdalsrsinfo "+proj=hammer +wktext"


The Misc Answers




  • The small differences in GDAL's PROJCS WKT and ESRI's PROJCS WKT are enough to cause problems for either ArcGIS or GDAL/OGR-dependent apps (e.g., GeoServer). When you encounter "weird" problems, this is a good place to troubleshoot.





  • GDAL/OGR and PROJ.4 have different "default" parameter values for quite a number of projections if the User did not specify them. This is another source of confusion because User normally refers to PROJ.4's documentation, while GDAL/OGR does not state this anywhere except in the source codes.




  • Strictly-speaking, it is OGR that speaks PROJCS WKT and uses PROJ.4.




Follow Up Comment from Questioner



Quick reaction: does not work either with the +wktest parameter, but the error message is different. Now it complains the file format is not recogised. I am feeding it the same file I fed to cs2cs>




The gdaltransform help page stated that srcfile is:


File with source projection definition or GCP's. If not given, source
projection is read from the command-line -s_srs or -gcp parameters

Unfortunately, this wasn't helpful. Running gdalinfo --debug on hints the reason for wgs.md's failure.


enter image description here


It appears that GDAL does not have support for plain text file. OGR does not have support for plain text file too although it does have support for CSV. (On the other hand, cs2cs and proj are part of the PROJ.4 project.)


The screen shot below shows the result I obtained using the "WGS84" coordinates as used in your example, ie 15 deg East, 15 deg South:-


enter image description here


Usually, I like to keep in mind that GDAL tools are mainly for rasters while OGR tools are for vectors.



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...