Convert LR/LOLA data to latitude and longitude

Question

I have height map data from here LR LOLA DEM from the Moon which seems to be "simple spherical" and I want to convert them to latitude and longitude with this: (It's a bit pseudo code)

from affine import Affine
import rasterio
from pyproj import Proj, transform, Transformer, CRS

infile = rasterio.open(GEOTIFF)
p = (0, 0)  # pixel position

srcProj = Proj(infile.crs)
destProj = ???  # what goes here?

transformer = Transformer.from_crs(srcProj, destProj)

p = p * infile.transform

lo, la = transformer.transform(p[0], p[1])

What destination projection should I use? I tried different ones, but got various exceptions (for example: src and dest are not the same celestial body, etc.)

Addendum: Python describes the infield.crs with the following, maybe that helps:

PROJCS["SimpleCylindrical Moon",
GEOGCS["GCS_Moon",DATUM["D_Moon",SPHEROID["Moon",1737400,0]],
PRIMEM["Reference_Meridian",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Equirectangular"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",0],
PARAMETER["standard_parallel_1",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]]]

Answer 1

Using pyproj.CRS you can get the geodetic CRS of the projected CRS.

from pyproj import CRS

crs = CRS(infile.crs)
print(repr(crs.geodetic_crs))

<Geographic 2D CRS: GEOGCRS["GCS_Moon",DATUM["D_Moon",ELLIPSOID["Moon" ...>
Name: GCS_Moon
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- undefined
Datum: D_Moon
- Ellipsoid: Moon
- Prime Meridian: Reference_Meridian

If you want to reproject the raster to a new grid based on latitude & longitude, then I would follow either one of these instructions:

• https://rasterio.readthedocs.io/en/stable/topics/reproject.html
• https://corteva.github.io/rioxarray/stable/examples/reproject.html

Here is an example using rioxarray:

import rioxarray
from pyproj import Transformer


rds = rioxarray.open_rasterio(GEOTIFF)
rds_latlon = rds.rio.reproject(crs.geodetic_crs.to_wkt("WKT1_GDAL"))
rds_latlon.rio.to_raster(LATLON_GEOTIFF)

If you want to reproject the point's coordinates and don't care about the grid being a regularly spaced grid, then this should work:

rioxarray example converting all points:

import rioxarray
from pyproj import Transformer


rds = rioxarray.open_rasterio(GEOTIFF)
transformer = Transformer.from_crs(crs, crs.geodetic_crs)
lon, lat = transformer.transform(rds.x.values, rds.y.values)

rasterio example for a single point:

from pyproj import Transformer

transformer = Transformer.from_crs(crs, crs.geodetic_crs)
lon, lat = transformer.transform(*infile.xy(0, 0))