1

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?

9
  • 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. Commented Feb 25, 2020 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
    Commented Feb 25, 2020 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 Commented Feb 25, 2020 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
    Commented Feb 25, 2020 at 5:31
  • 1
    My apologies, sqrt(2) was correct
    – FSimardGIS
    Commented Feb 25, 2020 at 5:39

1 Answer 1

2

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.

1
  • 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
    Commented Feb 25, 2020 at 6:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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