2

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"]]]

1 Answer 1

3

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:

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))

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.