Seeking assist: reproject coordinate array values

1

I aim to reproject values from two arrays, ‘longitude’ and ‘latitude’ and treating each corresponding value as a coordinate pair reproject to a local map projection. I’m just wondering about the approach I use, and whether there’s a more idiomatic Python way to do this.

Here is a notebook that achieves my aim:

The NetCDF file in the notebook with longitude,latitude arrays shape = (2030, 1354) is downloaded from NASA ocean colour.

There’s code to open the file, and isolate the ‘/navigation_data’ group which has longitude and latitude.

Then

  • reshape the values into a matrix of longitude,latitude pairs
  • define a local projection from a central coordinate
  • reproject with pyproj in a function project_array()
  • set up ‘x’ and ‘y’ DataArrays with the same properties as ‘longitude’ and ‘latitude’ and fill with the projected values

My questions are, what is the most xarray-idiomatic way to do this? Is there anything especially wrong with the style here of how the reprojection is done?

I know that copying the longitude,latitude to fill with x,y is not ideal but I’m less concerned about that. Is there a better way to do the raw coordinate reprojection, but if there’s input into how the new arrays are added that will also help a lot.

Any feedback welcome, thanks! Michael Sumner

2

Have you looked at pyresample

3

I don’t see how it’s relevant because there’s no resampling here :folded_hands:

Though, it’s certainly tangentially useful in what I want this for, and I might use it to tell that part of the story sometime later :smiley:

4

There’s no need to transform back and forth between Xarray coordinate variables and a coordinate matrix I think. You could probably simply do something like this:

def project_xarray(
    xx_src: xarray.DataArray,
    yy_src: xarray.DataArray,
    source,
    target,
):
    transformer = Transformer.from_crs(source, target, always_xy=True)
    res = transformer.transform(xx_src, yy_src)
    xx_dest = xarray.DataArray(res[0], coords=xx_src.coords, dims=xx_src.dims)
    yy_dest = xarray.DataArray(res[1], coords=yy_src.coords, dims=yy_src.dims)
    return xx_dest, yy_dest

xx, yy = project_xarray(
    ds.longitude, ds.latitude, source="EPSG:4326", target="EPSG:3412"
)
ds = ds.assign_coords(xx=xx, yy=yy)
1 Like
5

Great! Thank you,that works. Here’s what I did in the end:

from pyproj import Transformer
import numpy

import xarray
import os

def project_xarray(
    xx_src: xarray.DataArray,
    yy_src: xarray.DataArray,
    source,
    target,
):
    transformer = Transformer.from_crs(source, target, always_xy=True)
    res = transformer.transform(xx_src, yy_src)
    xx_dest = xarray.DataArray(res[0], coords=xx_src.coords, dims=xx_src.dims)
    yy_dest = xarray.DataArray(res[1], coords=yy_src.coords, dims=yy_src.dims)
    return xx_dest, yy_dest

oc = "http://oceandata.sci.gsfc.nasa.gov/getfile/AQUA_MODIS.20250704T003501.L2.OC.NRT.nc"
#not sure how to automate this, but 11Mb file at this URL needs earthdata Auth
#!wget --header="Authorization: Bearer ey..." http://oceandata.sci.gsfc.nasa.gov/getfile/AQUA_MODIS.20250704T003501.L2.OC.NRT.nc

f = os.path.basename(oc)
coords = xarray.open_datatree(f).navigation_data

xx, yy = project_xarray(
    coords.longitude, coords.latitude, source="EPSG:4326", target="EPSG:3412"
)

gd = xarray.open_dataset(f, group = "geophysical_data")
gd.assign_coords(xx = xx, yy = yy)