ENH: Consider adding lat, lon and alt option in georeferencing

[syedhamidali]

• xradar version: 0.8.1
• Python version: 3.12.7
• Operating System: MacOS Sequoia 15.1 (24B83)

Description

The current georeferencing in xradar provides Cartesian coordinates (x, y, z), but there’s no direct option to include geographic coordinates (lat, lon, and alt). Adding this feature would enhance usability for users needing georeferenced datasets in geographic space.

What I Did

Attempted to use georeferencing with geo=True

geo = radar.xradar.georeference(geo=True)

[kmuehlbauer]

@syedhamidali Thanks for this proposal.

In #25 (comment) I've proposed to just implement radar coordinates as AEQD coordinates, and leave projecting/transforming coordinates to specialized libraries (eg. Py-ART/wradlib or even more projection dedicated ones like pyproj).

During enhancing plotting capabilities we've finally implemented a way to use the CRS coordinate (eg. spatial_ref, crs_wkt) via pyproj to have cartopy do the transformation.

There is also work still going on to build this CRS as Index accessor-based into xarray (pydata/xarray#2288), but it's hard to predict, when this will be in a usable state (need to re-check!)

From my perspective we should not implement the conversion/transformation into xradar as proposed in #244, but guide users how to do this with these external libraries (which they might use anyway in further processing). When using our new map_over_sweeps accessor this is only one or two additional lines of code.

If we (ping @openradar/developers) still want to have this functionality available from .xradar.georeference() I'd suggest a slightly different API:

ds = ds.xradar.georeference(crs=trg_crs, names={"x": "lon", "y": "lat", "z": "alt"})

Internally we would setup pyproj.Transformer and do the transform from xradars AEQD to any possible projection/coord system.

Here is the functionality as MCVE:

import numpy
import xarray as xr
import xradar
from pyproj import CRS, Transformer, transform
from open_radar_data import DATASETS

filename = DATASETS.fetch("cfrad.20080604_002217_000_SPOL_v36_SUR.nc")
radar = xd.io.open_cfradial1_datatree(filename, first_dim="auto")

swp = radar["sweep_0"].to_dataset().xradar.georeference()

src_crs = swp.xradar.get_crs()
trg_crs = CRS.from_user_input(4326)
names = {"x": "lon", "y": "lat", "z": "alt"}

transformer = Transformer.from_crs(src_crs, trg_crs)

(trg_x, trg_y, trg_z) = transformer.transform(swp.x, swp.y, swp.z)

swp = swp.assign_coords({
    names.get("x", "x"): (swp.x.dims, trg_x),
    names.get("y", "y"): (swp.y.dims, trg_y),
    names.get("z", "z"): (swp.z.dims, trg_z),
})
     
display(swp) 

• In the process we would need to set/override the attributes as needed for that particular item.
• We would also have to set/add/override the projection coordinate.

[kmuehlbauer]

Also pinging @openradar/xradar for feedback.

[syedhamidali]

src_crs = swp.xradar.get_crs()
trg_crs = CRS.from_user_input(4326)
names = {"x": "lon", "y": "lat", "z": "alt"}

transformer = Transformer.from_crs(src_crs, trg_crs)

(trg_x, trg_y, trg_z) = transformer.transform(swp.x, swp.y, swp.z)  

@kmuehlbauer
I like this approach. It’s simple and avoids adding extra code to the package. We can just include an example for it.

[kmuehlbauer]

src_crs = swp.xradar.get_crs()
trg_crs = CRS.from_user_input(4326)
names = {"x": "lon", "y": "lat", "z": "alt"}

transformer = Transformer.from_crs(src_crs, trg_crs)

(trg_x, trg_y, trg_z) = transformer.transform(swp.x, swp.y, swp.z)  

@kmuehlbauer I like this approach. It’s simple and avoids adding extra code to the package. We can just include an example for it.

The major pain point is the last two bullet points of my above comment which also have to be addressed. I have hope, that those issues would be fixed in a possible dedicated package as laid out in that other xarray issue pydata/xarray#2288.

[benbovy]

I have hope, that those issues would be fixed in a possible dedicated package as laid out in that other xarray issue pydata/xarray#2288.

There you go: https://github.com/benbovy/xproj :-) Still at the stage of a rough proof-of-concept, though. Feedback is welcome!

[kmuehlbauer]

Thanks @benbovy for advertising xproj. That's already on my watch-list. We'll definitely give it a try!

[dcherian]

I agree, would be interesting to explore the new CoordinateTransformIndex & Xproj. There seem to be many parallels with the RasterIndex work in https://github.com/dcherian/rasterix/

[mgrover1]

Fully agree here @dcherian 👍