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

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?

• I don't know if PROJ.4 supports a latitude of origin for the Mercator projection. Try adding +lat_0=-31.6997 and +lon_0=41.313 to the definition. As it is, the y origin is at the equation, thus the very different latitude result. The longitude value is also off a little so changing to the "lonor" value may fix that. – mkennedy Jul 28 '16 at 20:40
• Is it correct that all points are in the ocean south of Madagascar? – AndreJ Jul 30 '16 at 8:45
• yes, that is correct. – juh Jul 30 '16 at 11:01

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

projection_merc_tmp = " ".join(("+proj=merc +lon_0=",str(lon0),
"+k=1 +x_0=0 +y_0=0",
"+ellps=sphere",
"+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",
"+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 \$ \$

• You can use accent grave instead of \$ to mark variables to `x_0`. – AndreJ Jul 30 '16 at 17:28

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
``````
• Thank you very much. I've tried your solution, tmerc, using a sphere instead of the geoid, but it is not working for points far from the coordinates I gave at the beginning. The problem is: I know that the projection used is a Mercator Projection (it's written in the netcdf file i've taken data from) and I can retrieve the coordinates with the mathematical formulas. But, I'd like to make it work with pyproj so I can use the WGS84 geoid without loosing too much speed at computation. – juh Jul 30 '16 at 11:11
• See my revised answer. – AndreJ Jul 30 '16 at 13:48