0

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')

1 Answer 1

0

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

enter image description here

Initial geodataframe

enter image description here

Final geodataframe with spi

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.