High resolution x,y dataset for Mollweide projection convert to latitude and longitude coordinates

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. 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
• 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) Feb 25 '20 at 5:31
• My apologies, sqrt(2) was correct Feb 25 '20 at 5:39