0

I'm working in a xarray dataset with rectilinear grids (latitude constant over x-axis / longitude constant over y-axis) which have 2D coordinates. I needed to subset this dataset to a particular lat/lon value, but unfortunately subsetting a dataset according to lat/lon variables in multidimensional coordinates are a little less intuitive and became quite exhausting for me, since we cannot use ds.sel(lat=y, lon=x, method='nearest'), for example. Ultimately, I resort to ds.where() option but ended up having to find the nearest lat/lon coordinates in the dataset manually. Here's a print of my dataset:

xarray.DataArray 'u' (MT: 366, Depth: 1, Y: 266, X: 99)

[9638244 values with dtype=float32]

Coordinates:
    MT (MT) datetime64[ns] 2012-01-01 ... 2012-12-31
    Depth (Depth) float32 0.0
    Latitude (Y, X) float32 array([[-21.92 -21.92 ... -1.28 -1.28]])
    Longitude (Y, X) float32 array([[319.6 319.7 319.8 ... 327.4 327.4]])

Attributes:
    standard_name: eastward_sea_water_velocity
    units: m/s
    long_name: u-veloc. [90.9H]
    _ChunkSizes: [1 3 413 563]

As a result, I expect to convert the 2D coordinates into 1D lat/lon dimensions as below:

xarray.DataArray 'u' (MT: 366, Depth: 1, Lat: 266, Lon: 99)

[9638244 values with dtype=float32]

Coordinates:
    MT (MT) datetime64[ns] 2012-01-01 ... 2012-12-31
    Depth (Depth) float32 0.0
    Latitude (Lat) float32 -21.92 -21.92 ... -1.28 -1.28
    Longitude (Lon) float32 319.6 319.7 319.8 ... 327.4 327.4

Attributes:
    standard_name: eastward_sea_water_velocity
    units: m/s
    long_name: u-veloc. [90.9H]
    _ChunkSizes: [1 3 413 563]

After some research, I found the xESMF package which does a great job for regridding 2D curvilinear grids into 2D rectilinear grids, and also this thread comparing both xESMF and Scipy performances in terms of interpolating 2D to 1D coordinates. However, I didn't find any proper solution for converting these 2D coordinates into nominal 1D lat/lon dimensions. Is it possible to be done?

Gabriel Lucas
  • 170
  • 1
  • 2
  • 14
  • 1
    if the grid is evenly spaced, you convert the coordinate to index using division by the bin size. – Mad Physicist Jan 03 '22 at 21:56
  • 1
    In order for us to help you more effectively, please provide a [minimal reproducible example](https://stackoverflow.com/help/minimal-reproducible-example), including a small example input data and the corresponding expected result. – Pierre D Jan 03 '22 at 23:00
  • 1
    At least, please `print(ds)` and give us a more detailed definition of your grid spec including the `x, y` and `lat, lon` coordinates. additionally, you'll need to decide what you actually want to happen when reshaping your data. do you want nearest-neighbor remapping to the new grid? bilinear interpolation? something else? and what do you want your output lat/lons to be exactly? thanks! – Michael Delgado Jan 04 '22 at 02:10
  • Longitude 319.6? On what planet? – Tim Roberts Jan 04 '22 at 06:02
  • The longitudes in this dataset are given in decimal degrees from 0 to 360. To find the values between -180 to 180 you just need to subtract it from 360 (for western hemisphere), i.e. -39.4 or 39.4°W – Gabriel Lucas Jan 04 '22 at 06:25

0 Answers0