I have a 2D array that's 144000 x 64000 that has some values that are irrelevant to the question. Let's say the array can be used like an X, Y coordinate system. I somehow need to get the Latitude and Longitude based from the Mollweide projection.

I took a look at the Wiki inverse transformation formula, but I believe I'm doing the math incorrectly. Let's say I have some X,Y coordinate (130000, 50000), how would I get the latitude and longitude from it?

Note: the actual array is composed of several arrays that are 4000 x 4000, with the number of arrays being 36 x 16.

My current math is based on this SO question:

R = (36 * 4000) / 2 / sqrt(2) = 50911.68825 radians
theta = arcsin(50000 / R / sqrt(2) = 43.9829 radians
latitude = arcsin[(2 * theta + sin(2 * theta)) / pi] = DOMAIN ERROR 
    => since (2 * theta + sin(2 * theta)) / pi = 28.31853329
longitude = 130000 * pi / (2 * R * sqrt(2) * cos(theta)) = 3.941594775

So neither latitude nor longitude make ANY sense. Am I supposed to convert theta to something below 2 * Pi? What math am I doing incorrectly?

  • Are you sure you want to do this the hard way? An open source lib like GDAL (implementing OSR) would be much faster and a whole lot more trustworthy. Feb 25 '20 at 4:51
  • @MichaelStimson Please point me in the right direction :) I've been using rasterio for some population density .tif files and it returns these numpy arrays. Is there anything specific that GDAL can offer that helps? Rasterio also uses GDAL and this is my first time touching anything related to GIS.
    – acwpython
    Feb 25 '20 at 4:54
  • Isn't rasterio built on GDAL? Either way a point geometry can be generated, a spatial reference defined then projected, have a read of gdal.org/tutorials/osr_api_tut.html and see if that helps. This one too, for python pcjericks.github.io/py-gdalogr-cookbook/projection.html Feb 25 '20 at 4:56
  • 1
    Also, the projection origin is at the center of your image, not a corner, so you need to take that into account as well in the calculation by subtracting half of your image width and height to your pixel x,y coordinates before passing your coordinates in the mollweide formulas. Ex: to evaluate pixel (130000,50000) use (130000-72000,50000-32000) -> (58000,18000)
    – FSimardGIS
    Feb 25 '20 at 5:31
  • 1
    My apologies, sqrt(2) was correct
    – FSimardGIS
    Feb 25 '20 at 5:39

For the calculation of R, you will have to use half of your image width (144,000 / 2) = 72,000, which gives an R value of 25455.8441227.

Also, since the projection origin is at the center of the image, not a corner, you have to take this into account as well before passing your x,y values in the Mollweide formulas. Ex: to evaluate pixel (130000,50000) use (130000-72000,50000-36000) -> (58000,18000)

Also, image pixel values in the y axis are often computed from the top down, so that's something to be wary of. You might have to reverse the sign for the latitude value, if your pixel y coordinate origin is at the top of your image.

  • Note: From discussion, the 144000 x 64000 array is incorrect. The GHS-POP dataset that this question was based on had a visual that did not show Antarctica; there should've been a 2:1 ratio since this was a Mollweide projection. So the actual array should be 144000 x 72000.
    – acwpython
    Feb 25 '20 at 6:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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