gdal - Re-project raster image from Mercator to Equirectangular

I have a raster image of the world map in Mercator projection, square 1:1 format. I want transform the image into an Equirectangular projection accurately.

I assumed from my research that GDAL would be the way to go, specifically gdalwarp. I used gdaltranslate to add gcp points to the image on the four corners, and then used gdalwarp to make the transformation.

It produced an image with the correct size – 2:1 ratio – but the projection has not been accurately transformed. It looks exactly the same as if I had taken it into Photoshop and simply stretched it to that shape. The proportional relationship between Mercator and Equirectangular is obviously more complicated than a simple, global stretch like this.

So my question is: is there a way to accurately transform between projections for raster images?

Edit / Additional Information

Here is the starting image: http://imgur.com/48wtsgn

I then applied to GCP points using the following gdal_translate command:

gdal_translate -gcp 0.0 0.0 -180.0 90.0 -gcp 1000.0 0.0 180.0 90.0 -gcp 1000.0 1000.0 180.0 -90.0 -gcp 0.0 1000.0 -180.0 -90.0 Mercator.tif MercatorGCP.tif

Then I used the following gdalwarp command to make the intended reprojection from Mercator (ESPG:3857) to Equirectangular (ESPG:32662):

gdalwarp -s_srs EPSG:3857 -t_srs EPSG:32662 MercatorGCP.tif Equirectangular.tif

Here is the resulting image: http://imgur.com/5LGk16o

Compare it to a map with proper Equirectangular projection, like this: http://imgur.com/3ncIW2R

You can see that the proportions have not been properly adjusted, the map not actually reprojected at all, simply squashed uniformly.

Am I doing something wrong, or is this simply not possible with gdal? If not, is there another way to do this?

Thanks!

Answer

For the Mercator projection, the extent can not reach North and South pole for mathematical reasons.

The standard Google and Openstreetmap mercator projection is limited to 85.011° North and South to get a square map.

See http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#X_and_Y for explanation.

Using EPSG:3857, the extent of a map is -20037508,-20037508,20037508,20037508 in CRS units. Using degree coordinates at this point make little sense.

You have to change your GCP to reflect the proper extent coordinates.

---

EDIT 1

Adding the GCPs does not change the corner coordinates. Instead apply the correct bounds to the file with -a_ullr. This way it works:

gdal_translate -a_srs EPSG:3857 -a_ullr -20037508 20037508 20037508 -20037508 Mercator.tif MercatorULLR-CRS.tif
gdalwarp -s_srs EPSG:3857 -t_srs EPSG:4326 MercatorULLR-CRS.tif Equirectangular.tif

And the result looks like this:

Note that I have chosen EPSG:4326 as output. Really a geographic CRS, but displayes the same way a equirectangular projection would. The blue marble background exceeds the part between 85.011° and 90°.

---

EDIT 2

GDAL also works with GCP, but that requires a further step:

1. Apply the GCP


2. Transform the picture using the GCP to World Mercator


3. Reproject from World Mercator to your desired CRS

To apply the GCP correctly, you have to know the local coordinate systems used for native tif and your CRS. The tif has the origin in the upper left corner, X to the right, and Y to the bottom. The World Mercator has the origin at 0°E/0°N, X to the East, and Y to the North.

So the command lines are:

gdal_translate -gcp 0.0 0.0 -20037508 20037508 -gcp 1000.0 0.0 20037508 20037508 -gcp 1000.0 1000.0 20037508 -20037508 -gcp 0.0 1000.0 -20037508 -20037508 Mercator.tif MercatorGCP.tif
gdalwarp -r bilinear -t_srs EPSG:3857 MercatorGCP.tif MercatorCRS.tif
gdalwarp -t_srs EPSG:4326 MercatorCRS.tif Equirectangular.tif