Friday 23 October 2015

coordinate system - Hammer projection in QGIS


I would like to use Hammer's projection in QGIS. I am providing the follow proj string:



+proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs

But QGIS complains that it is invalid:


enter image description here


According to the proj documentation all parameters are optional, I would expected it to function like this.


I have QGIS 2.18.20 and proj 4.9.3.



Answer



While +hammer is implemented in PROJ.4, it is not in GDAL.


In a OSGEO4W shell or Linux terminal, you can put some lon lat coordinates in a textfile wgs.txt and run


cs2cs +init=epsg:4326 +to +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs wgs.txt >hammer.txt

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

The first one runs fine (using pure PROJ), while the second fails.


QGIS depends on GDAL, hence it can not use your projection string.


In most cases it is a matter of a missing inverse projection. Running


 proj -I +proj=hammer
cs2cs +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs +to +init=epsg:4326 hammer.txt >wgs84a.txt

with proj 4.8.0 returns an error in both lines:


 :

inverse projection not available
program abnormally terminated



UPDATE


It seems that the inverse projection has been implemented in Proj4 v4.9.3, so following an advice from http://osgeo-org.1560.x6.nabble.com/QGIS-Developer-EqualEarth-projection-support-in-QGIS-td5375972.html a projection not known to GDAL can be added to QGIS with


 +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +wktext +no_defs

which gives the expected image of the Earth:


enter image description here



and the command line tests work too:


 proj -I +proj=hammer +wktext  cs2cs +init=epsg:4326 +to +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +wktext +no_defs wgs.txt >hammer.txt
cs2cs +proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +wktext +no_defs +to +init=epsg:4326 hammer.txt >wgs84a.txt
gdaltransform -s_srs "+init=epsg:4326" -t_srs "+proj=hammer +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +wktext +no_defs" hammer.txt

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...