3

I want to know how to select values from the xarray DataArray based on the location (geo_df.geometry) and time (geo_df.plant_date & geo_df.cut_date) of rows in the geopandas GeoDataFrame. I want to join them as 'features' in an output GeoDataFrame.

My datasets:

Packages I'm using:

import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from shapely import geometry
import xarray as xr

I have a geodataframe storing lat/lon POINTS which corresponds to households. The index column is the id of the households.

geo_df.head()

Out[]:
  crop_name     xxx     cut_date plant_date                       geometry
0   SORGHUM  0.061029 2011-11-10 2011-11-10 POINT (37.89087631 14.35381619)
1    MILLET -0.104342 2011-10-19 2011-10-19 POINT (37.89087631 14.35381619)
2   SORGHUM -0.031697 2013-11-26 2013-11-26 POINT (37.89087631 14.35381619)

I have an xarray object storing GRIDDED vegetation health data (NDVI).

ndvi_df = xr.open_dataset(geo_data_dir+ndvi_dir).ndvi

Out[]: <xarray.DataArray 'ndvi' (time: 212, lat: 200, lon: 220)>
[9328000 values with dtype=float32]
Coordinates:
  * lon      (lon) float32 35.024994 35.074997 35.125 35.174988 35.22499 ...
  * lat      (lat) float32 14.974998 14.924995 14.875 14.824997 14.775002 ...
  * time     (time) datetime64[ns] 2000-02-14 2000-03-16 2000-04-15 ...
Attributes:
    long_name:   Normalized Difference Vegetation Index
    units:       1
    _fillvalue:  -3000

I have a geodataframe storing a POLYGON which corresponds to a country.

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
ethiopia = world.loc[world["name"] == "Ethiopia"]

Visual Summary:

My datasets plotted on top of one another look as follows (plotted annually for demonstration purposes).

(ndvi_df.loc[f'{year}-01-16T00:00:00.000000000':f'{year}-12-16T00:00:00.000000000']
 .mean(dim='time')
 .plot(cmap='gist_earth_r', vmin=-0.1, vmax=1)
)

ax = plt.gca()

ethiopia.plot(alpha=0.2, color='black', ax=ax)

(geo_df
 .loc[ (lsms_geo_1["cut_date"] > f'{year}-01-01') & (lsms_geo_1["cut_date"] < f'{year+1}-01-01') ]
 .plot(markersize=6 ,ax=ax, color="#FEF731")
)
ax.set_title(f'{year} Mean NDVI and Households')
plt.show()

Household data plotted on top of NDVI gridded product, with Ethiopia shapefile shaded.

Ideal Output:

I want as an output, a geodataframe with extra columns telling me the NDVI values in the PRECEDING MONTHS for the pixel which the households are inside.

The index column is the id of the households.

like this:

  crop_name     xxx     cut_date plant_date                       geometry  ndvi_month_0  ndvi_month_1  ndvi_month_2
0   SORGHUM  0.061029 2011-11-10 2011-11-10 POINT (37.89087631 14.35381619)          0.3           0.3           0.3
1    MILLET -0.104342 2011-10-19 2011-10-19 POINT (37.89087631 14.35381619)          0.6           0.6           0.6
2   SORGHUM -0.031697 2013-11-26 2013-11-26 POINT (37.89087631 14.35381619)          0.1           0.1           0.1

I would also like to know how to subset my data in xarray object by using the geodataframe polygon ethiopia.

NOTE: Reposted from stack overflow here because it seems a GIS related question.

1

You can do this with a combination of the standard pandas apply method and nearest neighbour lookups from xarray:

def ndvi_for_point(row):
    """Return a pandas series of ndvi values which will be indexed by the
       time index."""
    cut_date = row['cut_date']
    start_date = cut_date.replace(year=cut_date.year - 1)
    # use `.loc` to index by labels rather than position
    limited_ndvi = ndvi_df.loc[start_date: cut_date]

    point = row['geometry']

    return limited_ndvi.sel(lat=point.y, lon=point.x, method='nearest').to_series()

ndvi_extract = geo_df.apply(ndvi_for_point, axis=1)

The above code will generate a dataframe indexed by the point index, and with columns on the time index from the NDVI dataset, containing the data values from the netCDF file. apply on axis=1 means row by row, rather than column by column, which will give you access to both the cut date and the geometry.

In terms of masking the raster, if you're doing so simply to screen out points within Ethiopia, then you can simply use a geopandas spatial join to join the two vector data layers together:

gpd.sjoin(geo_df, ethiopia)

If you do need to mask the netCDF itself I'd suggest asking a separate distinct question, though you could start by looking at burning shapes into a raster with rasterio.

  • Thank you so much @om_henners. Is there anyway to speed it up by only selecting the ndvi_for_point for the 12 months before the cut_date? Ideally I would do this inside the ndvi_for_point() function. Most likely by calling another function to keep the computations separate. – Tommy Lees Jul 9 '18 at 14:10
  • 1
    @TommyLees Yes, you can slice by labels in xarray as well. Updated the above to show an example by replacing the year of the cut date with the year before. I should note as well, if you have different date ranges that come out because of different cut dates you'll get NaN values in your output dataframe. – om_henners Jul 10 '18 at 1:33
  • okay so is there a way to change the column headings to ndvi_cut_date_0, ndvi_cut_date_1, ndvi_cut_date_2 ... which would be consistent across all rows? I was trying to fix this yesterday but I had an awful time with ValueError: Plan shapes are not aligned after trying to reset the index to t-0, t-1, t-2 – Tommy Lees Jul 10 '18 at 7:37
  • 1
    okay fixed it with this (final three lines on your function)! series = limited_ndvi.sel(lat=point.y, lon=point.x, method='nearest') columns = [f"ndvi_month_{i}" for i in np.arange(len(series))] return pd.Series(series.values , index=columns) – Tommy Lees Jul 10 '18 at 7:51
  • @TommyLees Nice one. Glad to hear it – om_henners Jul 10 '18 at 10:36

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.