Monday, 13 April 2015

gdalwarp - Making geotiff and reprojecting to EPSG:4326 using GDAL?


My goal is to download a gif file from this webpage:


http://aviationweather.gov/data/iffdp/2104.gif


Then crop it to the known latlon bounds of 75N,45S,35W,70E, convert to geotif, and then project to EPSG:4326. The map says "Mer", but how does one figure out what EPSG to use in gdal when converting to geotiff? I've tried many: 3395,3857,3785, etc, but none seem to work properly. Here are the commands I've used:


gdal_translate -of GTiff -a_srs EPSG:3857 -a_ullr -35 75 70 -45 2104_cropped.png 2104_3857.tiff

gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 2104_3857.tiff 2104_4326.tiff

I get no warning with gdalwarp, but the image never warps, it stays the same. When i overlay them on the maps that expect EPSG:4326, the maps to don't line up.




Answer



Firstly, crop the image source (coords are expressed in pixels here) with:


gdal_translate -srcwin 115 18 1360 2156 2104.gif 2104_cropped.tif

Then, transform the known WGS84 coordinates of the upper left and lower right corners to the "WGS 84 / World Mercator" projection (EPSG:3395):


cs2cs +init=epsg:4326 +to +init=epsg:3395
-35 75
-3896182.18 12890914.14 0.00
70 -45
7792364.36 -5591295.92 0.00


Assign the "WGS 84 / World Mercator" projection and the coordinates of the upper left and lower right corners in the same SRS:


gdal_translate -a_srs EPSG:3395 -a_ullr -3896182.18 12890914.14 7792364.36 -5591295.92 2104_cropped.tif 2104_3395.tif

Finally, transform the image from the "WGS 84 / World Mercator" projection to EPSG:4326:


gdalwarp -t_srs EPSG:4326 2104_3395.tif 2104_4326.tif

This is the graphical result in gvSIG:


enter image description here


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