Python Xarray – Fixing Wrong Coordinates When Saving DataArray

noaapythonrasterrioxarrayxarray

I have three 3×3 arrays that are downsampled data from NOAA's Global Hydro Estimator (smaller to make it quicker for testing). One array is lats, one is lons, and one is rainfall. Here's a sample

lons: 
[[-180    0  180]
 [-180    0  180]
 [-180    0  180]]

lats:
[[  60   60   60]
 [   0    0    0]
 [ -60  -60 -600]]

rainfall:
[[10 50  0]
 [ 0 20 20]
 [10 90  0]]

I would like to convert this into a raster with EPSG:4326 CRS. In order to do that I've chosen xarray and created a DataArray and saved to a raster like so

ds = xr.DataArray(rainfall,
                  coords={'lat':(['y','x'],lat_grid),
                          'lon':(['y','x'],lon_grid)},
                  dims=['y','x'])

ds.rio.write_crs("epsg:4326", inplace=True)
ds.rio.to_raster('test.tif')

However, when I look at the raster in qgis the lat/lon are on a scale from 0,3 instead of their respective values in the array above. I've tried messing around with the coords and dims and swapping them with no luck, does anyone have any idea of how to accomplish this?

Best Answer

Your latitude and longitude are 2D. They need to be converted to 1D. Here are some suggestions: https://github.com/corteva/rioxarray/issues/209

Related Solutions

Python – xarray Slicing Across the Antimeridian

Xarray is (intentionally) ignorant of coordinate systems, so it has no special handling for cyclic coordinates such as longitude. Because your longitude array has only increasing values, xarray interprets selections like slice(40, -80) in the same way that x[i:j] works if x is a NumPy array and i > j >= 0, and thus returns an empty selection.

The easiest way to work around this is to use boolean indexing instead:

ds.sel(lon=(ds.lon < -80) | (ds.lon > 40))

The inner parentheses are important, because | has higher operator precedence than < in Python.

Alternatively (maybe if you're studying the Pacific ocean), you might find it convenient to adjust the coordinate system of your data so it is centered over the anti-meridian instead. This is straightforward to accomplish with roll:

ds_rolled = ds.assign_coords(lon=(ds.lon % 360)).roll(lon=(ds.dims['lon'] // 2))

[GIS] Indexing coordinates in order to get a value by coordinates from xarray

You can use xarray.where. See: http://xarray.pydata.org/en/stable/indexing.html

xds = xarray.open_dataset("TT_199501_daymean.nc")
xds.where((xds.lon==5.8252) & (xds.lat==46.9359), drop=True)

UPDATE:

You can find the nearest neighbor using a KDTree. See: https://stackoverflow.com/questions/10818546/finding-index-of-nearest-point-in-numpy-arrays-of-x-and-y-coordinates

Then, you can select the point using .isel(...).

UPDATE: https://github.com/xarray-contrib/xoak

Best Answer

Related Solutions

Python – xarray Slicing Across the Antimeridian

[GIS] Indexing coordinates in order to get a value by coordinates from xarray

Related Question