Using pyproj and Mercator: x_0=0 and y_0=0?

Question

I'm having a hard time figuring out how to project correctly from a Mercator projection to LAT/Lon.

Basically, I have a netCDF file with the following variables:

• LON0 = 42.0 => reference longitude for the conformal projection (°E)
• LAT0 = -30.0 => reference latitude for the conformal projection (°N)
• LONOR = 41.313 => longitude of point x=0,y=0 in the conformal projection (°E)
• LATOR = -31.6997 => latitude of point x=0,y=0 in the conformal projection (°N)

I know that the point (476000,428000) is associated with geographical coordinates : ( 46.2559 , -27.8431 )

However, I can't figure out how to retrieve it by using pyproj transform object. Here is the code I'm using:

from pyproj import Proj,transform

projection_input = " ".join( ("+proj=merc +lon_0=42.0",
                              "+k=1 +x_0=0 +y_0=0",
                              "+ellps=WGS84 +datum=WGS84 +units=m",
                              "+no_defs")
                           )
inProj = Proj(projection_input)

projection_output = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
outProj = Proj(projection_output)

print transform(inProj, outProj, 476000, 428000)

Which gives me: ([46.27598102655299], [3.8677438926911805])

It has something to do with x_0,y_0 but I don't know what to put in it so I can get the right transformation (projection). How can this be solved by using the variables (only) listed at the top?

Answer 1

Well, I finally got it to work.

When I tried to perform a Mercator projection on the sphere, I forgot to take into account the latitude reference (which is -30.0 here). Therefore, the earth's radius is not 6371229m (the one used where I work) but rather, with (see picture)

** /!!!!\ ** HOWEVER, lat0 must not be taken into account when computing Mercator projection coordinates and transformation. The latitude reference will always be the equator. In this case, lat0 is only used to compute the radius of the spheroid used instead of using WGS84 geoid. ** /!!!!\ **

Once we know x_0 and y_0 and the right geoid or spheroid with correct +a and +b parameters, everythin works!

Here is the python code for retrieving x_0 and y_0:

from pyproj import Proj as ppProj
from pyproj import transform as pptransform
import numpy as np

lon0 = 42.0
lat0 = -30.0
lonor = 41.31303595758687
lator = -31.6996886147616
rad_earth = 6371229.0 * np.cos(np.radians(lat0))

projection_merc_tmp = " ".join(("+proj=merc +lon_0=",str(lon0),
                      "+k=1 +x_0=0 +y_0=0",
                      "+ellps=sphere",
                      "+a=",str(rad_earth),
                      "+b=",str(rad_earth),
                      "+units=m +no_defs")
proj_tmp = ppProj(projection_merc_tmp)
laloProj = pProj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
x0,y0 = pptransform(laloProj, proj_tmp, lonor, lator)
x0 = abs(x0)
y0 = abs(y0)

Then we can compute the correct projected coordinates for a given lat,lon:

projection_merc_str = " ".join(("+proj=merc +lon_0=",str(lon0),
                      "+k=1",
                      "+x_0=",str(x0),
                      "+y_0=",str(y0),
                      "+ellps=sphere",
                      "+a=",str(rad_earth),
                      "+b=",str(rad_earth),
                      "+units=m +no_defs")
                      )
proj_end = ppProj(projection_merc_str)

And

x,y = pptransform(laloProj, proj_end, 46.17279, -27.91653)

gives the correct result: x,y = 468000.,420000. (Long story short, I had to add 4km to each x,y to be able to compare those values with lat,lon values in the netcdf file (because of the 8km resolution and grid used by the model which produced the data).

Thanks a lot for all the answers, problem solved (hopefully).

PS: thanks AndreJ, I didn't know that one! Thanks a lot.

PS2: Well it seems that I cannot display math characters with $ $

Answer 2

While Proj.4 does not support Mercator projections with a latitude of origin different from the equator, there is a Mercator (2SP) projection that allows for a latitude of true scale different from the equator.

I have set up a file with the degree coordinates:

42 -30
41.31303595758687 -31.6996886147616
46.2559 -27.8431 

and one with the projected coordinates:

0 0
476000 428000

to run cs2cs on them.

Start off with a projection on the given center:

echo WGS84 - merc0 >out.txt
cs2cs +init=epsg:4326 +to +proj=merc +lat_ts=-30 +lon_0=42 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs -f "%%.4f" degree.txt >>out.txt
echo merc0 - WGS84 >>out.txt
cs2cs +proj=merc +lat_ts=-30 +lon_0=42 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +to +init=epsg:4326 -f "%%.4f"  merc.txt >>out.txt

which outputs:

WGS84 - merc0 
0.0000  -3018190.8809 0.0000
-66286.0745 -3208283.8198 0.0000
410635.9601 -2781624.9076 0.0000
merc0 - WGS84 
42.0000 0.0000 0.0000
46.9333 4.4612 0.0000

To get the right coordinates for false Easting and Northing, invert the sign of the coordinates in the second line, and add them as x_0 and y_0:

echo WGS84 - merc2 >>out.txt
cs2cs +init=epsg:4326 +to +proj=merc +lat_ts=-30 +lon_0=42 +x_0=66286.0745 +y_0=3208283.8198 +ellps=WGS84 +datum=WGS84 +units=m +no_defs -f "%%.4f" degree.txt >>out.txt
echo merc2 - WGS84 >>out.txt
cs2cs +proj=merc +lat_ts=-30 +lon_0=42 +x_0=66286.0745 +y_0=3208283.8198 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +to +init=epsg:4326 -f "%%.4f"  merc.txt >>out.txt

and you get:

WGS84 - merc2 
66286.0745  190092.9389 0.0000
-0.0000 -0.0000 0.0000
476922.0346 426658.9122 0.0000
merc2 - WGS84 
41.3130 -31.6997 0.0000
46.2463 -27.8307 0.0000

This is less than 2 km away from your desired point.

Using a sphere with R=6371229 meters takes it even closer:

cs2cs +proj=longlat +R=6371229 +to +proj=merc +lat_ts=-30 +lon_0=42 +x_0=66155.3878 +y_0=3221546.0900 +R=6371229 +units=m +no_defs -f "%%.4f" degree.txt >>out.txt

476003.2077 427999.9954 0.0000