Thursday, 5 November 2015

gdal - Transform an image with projection EPSG:4326 to EPSG:3857


I have an image with known boundary and projection EPSG:4326 and want to overlay it to openstreetmap using Leaflet.js.


I have tried using function imageoverlay, but the image is not aligned with the basemap.


L.imageOverlay(imgUrl, [[-15, 75], [45, 145]], {opacity: 0.6, autoZIndex: true});

So, I use gdalwarp to transform the image projection to EPSG:3857.


gdal_translate -of Gtiff -a_ullr -15 145 45 75 -a_srs EPSG:4326 test.png out_4326.tiff

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -ts 850 728 out_4326.tiff out_3857.tiff
gdal_translate -of png out_3857.tiff final.png

However, the final image looks like distorted. Here is the original image. enter image description here


Here is the final image enter image description here



Answer



gdal_translate expects the extent in the format -a_ullr ulx uly lrx lry. Form the picture I guess you swapped x (Easting) and y (Northing), and upper left and lower right.


I get the right picture with:


gdal_translate -of Gtiff -a_ullr 75 45 145 -15 -a_srs EPSG:4326 CJGMY.png out_4326.tiff
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -ts 850 728 out_4326.tiff out_3857.tiff

gdal_translate -of png out_3857.tiff final.png

fitting to Openstreetmap background:


enter image description here


For the correct leaflet syntax, see http://leafletjs.com/reference.html#imageoverlay


var imageUrl = 'http://www.lib.utexas.edu/maps/historical/newark_nj_1922.jpg',
imageBounds = [[40.712216, -74.22655], [40.773941, -74.12544]];

For an image of Newark NJ, this seems to be lower left, upper right in Northing-Easting order.


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