Create new column containing values of a xarray variable selected by point coordinates

Question

I have a xarray and a geodataframe. I want to extract spi_1 values of xarray from point geometries of geodataframe by creating a new column spi_1 in geodataframe.

print(xr)

<xarray.DataArray 'spi_1' (latitude: 130, longitude: 90)>
array([[        nan,         nan,         nan, ...,         nan,
                nan,         nan],
       [        nan,         nan,         nan, ...,         nan,
                nan,         nan],
       [        nan,         nan,         nan, ...,         nan,
                nan,         nan],
       ...,
       [-1.33405736, -1.34589297, -1.37533219, ..., -1.17554468,
        -1.18997407, -1.2228405 ],
       [-1.31988994, -1.35825529, -1.37751363, ..., -1.15944975,
        -1.17545966, -1.24460598],
       [-1.31240037, -1.34033954, -1.33462096, ..., -1.11622841,
        -1.12031903, -1.16772469]])
Coordinates:
  * latitude   (latitude) float32 4.725 4.775 4.825 4.875 ... 11.07 11.12 11.17
  * longitude  (longitude) float32 -3.275 -3.225 -3.175 ... 1.075 1.125 1.175
    time       datetime64[ns] 2022-02-01
    spi_1      (latitude, longitude) float64 nan nan nan ... -1.116 -1.12 -1.168

print(gdf)

    latitude  longitude  area_in_meters  confidence                  geometry
0   5.465193  -0.838154         18.8943      0.7155  POINT (-0.83815 5.46519)
1   6.153525  -2.062720         21.3032      0.8067  POINT (-2.06272 6.15353)
2   5.293025  -0.820021         17.8368      0.7078  POINT (-0.82002 5.29303)
3   5.532924  -1.489137         22.5540      0.6939  POINT (-1.48914 5.53292)
4   5.969935  -1.867445         14.5647      0.7480  POINT (-1.86744 5.96993)
..       ...        ...             ...         ...                       ...
84  5.432539  -1.310526         18.9966      0.6254  POINT (-1.31053 5.43254)
85  5.543290  -0.933647         11.2356      0.7201  POINT (-0.93365 5.54329)
86  6.028192  -1.803719         28.8116      0.8215  POINT (-1.80372 6.02819)
87  5.476935  -0.895145         22.8107      0.8535  POINT (-0.89514 5.47693)
88  5.592796  -0.804110         33.4757      0.7537  POINT (-0.80411 5.59280)

I tried the code below but it only selected points with coordinates near gdf points inside the xarray.

xr.sel(longitude=gdf['longitude'], latitude=gdf['latitude'], method='nearest')

Answer 1

This is not the solution but it was the best I could do.

import rioxarray
from shapely.geometry import mapping
from shapely.geometry import Point

# add projection system to nc
xr= xr.rio.write_crs('EPSG:4326', inplace=True)

# create new xarray containing spi_1 values only for selected by building coordinates
xr_spi = xr.rio.clip(gdf.geometry.apply(mapping), gdf.crs, drop=False)

# convert xarray to dataframe and drop missing values
df_spi = xr_spi.to_dataframe().dropna().reset_index()

# drop duplicated spi_1 columns
df_spi = df_spi.T.drop_duplicates().T

# convert dataframe to geodataframe
geometry = [Point(xy) for xy in zip(df_spi['longitude'], df_spi['latitude'])]
crs = 'EPSG:4326'
gdf_spi = gpd.GeoDataFrame(df_spi, crs=crs, geometry=geometry)

Initial xarray and geodataframe

Initial geodataframe

Final geodataframe with spi

The next step into merging the gdf with gdf_spi might be adding a buffer to gdf_spi and merging it with gdf