I have coordinates on a mercator grid projection from a NetCDF (CF) file. The metadata attributes in the file include : standard_parallel, longitude_of_projection_origin, earth_shape (spherical), earth_radius.

I want to use the gdal python API to transform the coordinates into lat-long. I need to set up the mercator coordinate system which I presume is best accomplished using SetMercator2SP. It takes arguments StdP1, CenterLat, CenterLong, FalseEasting and FalseNorthing.

How do I map these to the data I have?

I assume StdP1 is the standard parallel.

Is the longitude of projection origin the CenterLong?

What about CenterLat ... is that typically assumed to be zero if not listed?

Do I ignore false easting and northing?

  • Check this out: remotesensing.org/geotiff/proj_list/mercator_1sp.html Notably "PROJ.4 does not support a latitude of natural origin other than the equator", otherwise I thinkyour presumptions are right. (Apart from ignoring false easting/northing - you cannot do that if not zero, but it's a trivial adjustment). Anyway, in my experience it's rare for NetCDF to explicitly store non-longlat explicitly, so keen to learn more if you can share more details from the header. – mdsumner Nov 13 '14 at 12:02

Having worked this out, I should point out a flaw in my question: there is no SetMercator2SP in the python bindings. There is only one SetMercator and you use it to emulate SetMercator1SP or SetMercator2SP by choosing arguments. From what I have seen online, it is possible there are some programming mistakes ... not suprising since this isn't the most common projection.

The function SetMercator has arguments CenterLat, CenterLong, Scale, FalseEasting, FalseNorthing. The argument "CenterLat" seems to be the standard parallel (latitude of unit scale) not an origin, and I believe it is mutually exclusive of Scale. The CenterLong is the longitude of origin. The argument Scale would be ignored if you are emulating the "2SP" arguments, since use of the standard parallel and scale are mutually exclusive.

Note the comment by mdsummer about the only latitude of origin being the equator. My standard parallel was 30 degrees. Since this is substantially above the equator, this means that the data lying right at the standard parallel had a y coordinate significantly greater than zero (zero would be the value of the projection at the equator). However, my data came in this way ... the y coords at the bottom of my grid were appropriately above zero as if the equator were the y origin. Therefore I did not need a FalseNorthing.

Given the ambiguity, in the end I did this directly in proj4:

import pyproj
p=pyproj.Proj(r'+proj=merc +lat_ts=30.0 +lon_0=-128.735 +a=6370997.0 +b=6370997.0') 
print p(x_merc,y_merc,inverse=True)

Note that this use of a and b with these values is identical to "+ellps=sphere" which would shorten the command. However, in my case I had a slightly different earth_radius and I had to sub that value in.

Not the answer you're looking for? Browse other questions tagged or ask your own question.