This notebook teaches you how to read satellite imagery (Sentinel-2) from Google Earth Engine together with other data, e.g. WorldClim, and aggregate them for crop yield modeling in four districts of Nepal with a RandomForest Regression and XGBoost.

Nepal data courtesy of ICIMOD.

Setup the Notebook

!pip install -q geopandas shapely scikit-learn treeinterpreter rasterio rasterstats folium

%matplotlib inline
import os 
from os import makedirs, path as op
import matplotlib
import numpy as np
import rasterio
import geopandas as gpd
import pandas as pd
import shapely
import folium
import json
import time
import matplotlib.pyplot as plt

from sklearn import preprocessing
import xgboost as xgb
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import GridSearchCV 
from sklearn.model_selection import train_test_split

# If not on Colab you'll need install the earth-engine Python API
#!pip install earthengine-api #earth-engine Python API

import ee

from google.colab import drive
drive.mount('/content/drive/', force_remount=True)
# your root directory for outputs is set to your google drive
# you should create a sub-folder called data under your 'Colab Notebooks'.
# An example is in https://developmentseed.org/sat-ml-training/GettingStarted#Explore-your-drive
my_root_dir = "/content/drive/My Drive/Colab Notebooks/data"
crop_data = '/content/drive/Shared drives/servir-sat-ml/data/crop_yield/'
# crop_data = '/content/drive/My Drive/SERVIR-HKH-crop-yield/crop_yield'

Mounted at /content/drive/

Wheat yield data in 2018, Nepal

The following dataset contain 145 field measurement of wheat yield data.

# reading collected wheat productivity data from excel.
wheat_2018 = pd.read_excel(op.join(crop_data, 'Wheat yield_Field_Demonstrartion_sites_2018.xlsx'))
wheat_2018 = wheat_2018.drop(columns=['Farmer'])
wheat_2018.head(1)

# only keep the columns that interest us;
wheat_2018 = wheat_2018[['Lat', 'Long', 'District','wheat variety', 'Treat', '1000 Grain wt after oven dry + envelop', 'PH', 'OM']]
# rename the column names to be more readable
wheat_2018 = wheat_2018.rename(columns={'wheat variety': 'wheat_variety', '1000 Grain wt after oven dry + envelop':'wheat_yield'})
# turn the pandas dataframe into geopandas dataframe
wheat_2018_pts = gpd.GeoDataFrame(wheat_2018, geometry=gpd.points_from_xy(wheat_2018.Long, wheat_2018.Lat))
wheat_2018_pts.head(3)

wheat_2018_pts.to_file("wheat_2018_pts.shp")

Basic stats of the wheat yield field measurements

## the field measurements located in these few districts
pd.Series(wheat_2018.District).unique()

array(['Kanchanpur', 'Kailali', 'Bardiya', 'Banke'], dtype=object)

pd.Series(wheat_2018.wheat_variety).unique()

array(['Banganga', 'Tilganga', 'NHS 1755', 'NL 971', 'Vijay', 'Aditya',
       'Local', 'Tilottama', 'Unknown', 'Indian Hybrid 2285', 'Gautam'],
      dtype=object)

# see wheat yield differences by the districts by wheat variety
only_yield = wheat_2018.drop(columns=['Long', 'Lat', 'PH', 'OM', 'geometry'])

only_yield.groupby(by=['District', 'wheat_variety']).agg({'wheat_yield': ['mean', 'min']})

# see wheat yield differences by the districs by wheat treatment
only_yield.groupby(by=['District', 'Treat']).agg({'wheat_yield': ['mean', 'min', 'max']})

Quiz

What pattern do you see from above groupby functions over wheat yield?

Visualizing the field measurement on an interactive map

terai_bbox = gpd.GeoDataFrame(wheat_2018_pts).total_bounds
x_map=(terai_bbox[0] + terai_bbox[2])/2
y_map=(terai_bbox[1] + terai_bbox[3])/2
print(f'bouding box for Terai area is {terai_bbox}, with centroid point {x_map},{y_map}')

bouding box for Terai area is [80.165 27.999 81.648 28.986], with centroid point 80.9065,28.4925

mymap = folium.Map(location=[y_map, x_map], zoom_start=11,tiles=None)
folium.TileLayer('Stamen Terrain',name="Light Map",control=False).add_to(mymap)

tooltip = 'Click me!'
for lat, lon,  wheat_variety, wheat_yield in zip(wheat_2018.Lat, wheat_2018.Long, wheat_2018.wheat_variety, wheat_2018.wheat_yield):
  # print(lat, lon,  wheat_variety)
  folium.Marker([lat,lon], popup=f'Wheat variety: {wheat_variety}, \n Yield {wheat_yield}', tooltip=tooltip).add_to(mymap)

mymap

Sentinel-2 and Worldclim data extraction from Google Earth Engine

# Athenticate to your GEE account. 
!earthengine authenticate

# Earth Engine Python API
ee.Initialize()

Sentinel-2 imagery bands and indices composition

The purpose of this excersice is we want to have bands and indices information that may be helpful to modeling wheat yiled from temporal and spatial perspective

# From GEE
# If you are interested in using Landsat8 surface reflectance, see data over here https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_SR#bands
# The following steps showing how to extrac Sentinel-2 from GEE
aoi = ee.Geometry.Rectangle(terai_bbox.tolist())
band_sel = ('B2', 'B3', 'B4', 'B8', 'B11', 'B12')

sentinel_scenes = ee.ImageCollection("COPERNICUS/S2")\
    .filterBounds(aoi)\
    .filterDate('2018-01-01', '2018-12-31')\
    .select(band_sel)\
    .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE',10))

scenes = sentinel_scenes.getInfo()
[print(scene['id']) for scene in scenes["features"]]

sentinel_mosaic = sentinel_scenes.mean().rename(band_sel)

sentinel_mosaic.getInfo()

This next cell will export a Geotiff for each point of interest, by creating a bounding box around each point and then clipping the Sentinel layer to that area.

# crop_data = op.join(my_root_dir, "mosaic_s2")
"Colab notebooks/data/mosaic_s2"

Quiz

Let's extract seasonal Sentinel-2 instead of annual data

Run the Export (Optional)

 # This task will run in the background even if you close this notebook.
  # You can also check on the status of the task through the Javascript GEE interface
  # https://code.earthengine.google.com
  task.start()

# If you want to keep track of the export you can run this code
# However if run this, you will need to wait for it to finish before running additional code

while task.active():
  print('Polling for task (id: {}).'.format(task.id))
  time.sleep(15)

Visualizing Sentinel-2 imagery RGB in 2018 in our study area

# To make a map we first need some helper functions

# Define the URL format used for Earth Engine generated map tiles.
EE_TILES = 'https://earthengine.googleapis.com/map/{mapid}/{{z}}/{{x}}/{{y}}?token={token}'

#@title Mapdisplay: Display GEE objects using folium.
def Mapdisplay(center, dicc, Tiles="OpensTreetMap",zoom_start=8):
    '''
    :param center: Center of the map (Latitude and Longitude).
    :param dicc: Earth Engine Geometries or Tiles dictionary
    :param Tiles: Mapbox Bright,Mapbox Control Room,Stamen Terrain,Stamen Toner,stamenwatercolor,cartodbpositron.
    :zoom_start: Initial zoom level for the map.
    :return: A folium.Map object.
    '''
    mapViz = folium.Map(location=center,tiles=Tiles, zoom_start=zoom_start)
    for k,v in dicc.items():
      if gpd.geodataframe.GeoDataFrame == type(v):
        folium.GeoJson(
        data = v,
        name = k
          ).add_to(mapViz) 
      elif ee.image.Image in [type(x) for x in v.values()]:
        folium.TileLayer(
            tiles = v["tile_fetcher"].url_format,
            attr  = 'Google Earth Engine',
            overlay =True,
            name  = k
          ).add_to(mapViz)
      else:
        folium.GeoJson(
        data = v,
        name = k
          ).add_to(mapViz)
    mapViz.add_child(folium.LayerControl())
    return mapViz

s2_vis_params = {
  'bands': ['B4', 'B3', 'B2'],
  'min': 0,
  'max': 3000,
}

Mapdisplay(center=[y_map, x_map],
           dicc={'S2':sentinel_mosaic.getMapId(s2_vis_params)},
                #  'TrainingData':wheat_2018_disctrict}, 
           zoom_start=11)

Extract average temperature and precipitation from WorldClim

To enrich the attributes data for crop yield modeling, you can go to WorldClim to explore data including:

Minimum temperature
Maximum temperature
Average temperature
Precipitation
Solar radiation
Wind speed
Water vapor presure

But only average temperature and precipitation are available in GEE. More information of the dataset on GEE, please look over here.

# worldclim data on GEE
# https://developers.google.com/earth-engine/datasets/catalog/WORLDCLIM_V1_MONTHLY#bands
wcband_sel = ['tavg', 'prec']
worldclim = ee.ImageCollection("WORLDCLIM/V1/MONTHLY").select(wcband_sel)
worldclim_mosaic = worldclim.mean().rename(wcband_sel)

wheat_2018.head()

Quiz

Let's extract other dataset under WorldClim, e.g. max temperature, min temperature

Combine Sentinel-2 and Worldclim data

# Helper functions
def normalized_arr(ls_data):
  """normalized data 
  Args:
      np_array: data in a list
  Returns:
      x_scaled: rescaled/nomalized data in the value range between 0 and 1
  """
  arr = np.array(ls_data)
  #normalized the data 
  min_max_scaler = preprocessing.MinMaxScaler()
  x_scaled = min_max_scaler.fit_transform(arr)
  return x_scaled

def construct_indices(band1, band2):
  """compute band indices
  Args:
      band1, band2: specific band data 
  Returns:
      ind: index, e.g. ndvi
  """
  ind = (band1 - band2)/(band1 + band2)
  return ind

Sentinel-2 data reconstruction

wheat_2018_pts.head()

wheat_2018_gdf = wheat_2018_pts.geometry.buffer(0.05)
wheat_2018_gdf['geometry_x'] = wheat_2018_pts.geometry.buffer(0.05)
wheat_2018_gdf.head()

0    POLYGON ((80.24300 28.98000, 80.24276 28.97510...
1    POLYGON ((80.24300 28.98000, 80.24276 28.97510...
2    POLYGON ((80.24300 28.98000, 80.24276 28.97510...
3    POLYGON ((80.23800 28.98600, 80.23776 28.98110...
4    POLYGON ((80.23800 28.98600, 80.23776 28.98110...
dtype: geometry

# Reference the raster on disk.
bands = []
for i in range(len(wheat_2018)):
  raster = op.join(crop_data,'mosaic_s2',f'S2_mosaic_Terai_Nepal_2018_{i}.tif')
  with rasterio.open(raster, 'r') as src:
    frsc = src.read()
    print(frsc.mean(axis=(1, 2)))
    bands.append(frsc.mean(axis=(1, 2)))

band_arr = normalized_arr(bands)
bands_df = pd.DataFrame(band_arr, columns=['B2', 'B3', 'B4', 'B8', 'B11', 'B12'])

bands_df['ndvi'] = construct_indices(bands_df.B8, bands_df.B4)
bands_df['ndwi'] = construct_indices(bands_df.B8, bands_df.B3)
bands_df.head()

len(bands_df)

145

Adding Worldclim data

wcs = []
for i in range(len(wheat_2018)):
  raster = op.join(crop_data,'worldclim',f'wc_monthly_{i}.tif')
  with rasterio.open(raster, 'r') as src:
    frsc = src.read()
    print(frsc.mean(axis=(1, 2)))
    wcs.append(frsc.mean(axis=(1, 2)))

wc_arr = normalized_arr(wcs)
wcs_df = pd.DataFrame(wc_arr, columns=['tavg', 'prec'])
wcs_df.head()

len(wcs_df)

145

Concatenate Sentinel-2 and Worldclim to the wheat yield dataframe

dfs = [wheat_2018, bands_df, wcs_df, soil_df, indices_df]
wheat_2010_all = pd.concat(dfs, axis=1)
wheat_2010_all.head()

wheat_2010_all

len(wheat_2010_all)

145

Quiz

If you want to include both 'wheat_variety' and 'Treat' in the regressor, what should you do next?

# write you update code here

Crop yield modeling with XGBoost

Beginer friendly XGBoost regressor tutorial

More parameter search and tunning tutorial

Best ways to evaluate regression model: 3 Best metrics to evaluate Regression Model?

Train Model Using k-fold Cross Validation

X = wheat_2010_all[['PH', 'OM', 'B2', 'B3', 'B4', 'B8', 'B11', 'B12', 'ndvi', 'ndwi', 'tavg', 'prec']]
y= wheat_2010_all.wheat_yield

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=123)

Configure the Model

# learn more about the regressor parameter https://xgboost.readthedocs.io/en/latest/parameter.html#learning-task-parameters
xg_reg= xgb.XGBRegressor(objective = 'reg:squarederror' , 
                         colsample_bytree=0.4,
                         gamma=0,                 
                         learning_rate=0.07,
                         max_depth=3,
                         min_child_weight=1.5,
                         n_estimators=10000,                                                                    
                         reg_alpha=0.75,
                         reg_lambda=0.45,
                         subsample=0.6,
                         seed=42)

Fit the Model then Predict

xg_reg.fit(X_train,y_train)

preds = xg_reg.predict(X_test)

Evaluate the Model

rmse = np.sqrt(mean_squared_error(y_test, preds))
print("RMSE: %f" % (rmse))

RMSE: 7.058864

Visualize Metrics

Now we can plot y_test (test data) vs preds (predictions). Then we can investigage the variable importance to the model.

plt.scatter(preds, y_test, c = "green", marker = "s")
plt.title("Scatter plot for y_test and prediction")
plt.xlabel("Prediction")
plt.ylabel("True labels")
plt.show()

xgb.plot_importance(xg_reg)
plt.rcParams['figure.figsize'] = [10, 20]
plt.show()

Quiz

Try to use grid search and parameter search and tunning to tune the model

# write your code here

# provide script here
# this will take a bit of time to run, like 30mins
parameters_for_testing = {
 'colsample_bytree':[0.4,0.6,0.8],
 'gamma':[0,0.03,0.1,0.3],
 'min_child_weight':[1.5,6,10],
 'learning_rate':[0.1,0.07],
 'max_depth':[3,5],
 'n_estimators':[5000, 10000, 15000, 20000],
 'reg_alpha':[1e-5, 1e-2, 0.75],
 'reg_lambda':[1e-5, 1e-2, 0.45],
 'subsample':[0.6,0.95] 
}
 
xgb_model = xgb.XGBRegressor(learning_rate =0.1, n_estimators=1000, max_depth=5, min_child_weight=1, gamma=0, subsample=0.8, colsample_bytree=0.8, nthread=6, scale_pos_weight=1, seed=27)
gsearch1 = GridSearchCV(estimator = xgb_model, param_grid = parameters_for_testing, n_jobs=4,iid=False, verbose=10,scoring='neg_mean_squared_error')
gsearch1.fit(X_train,y_train)
print (gsearch1.grid_scores_)
print('best params')
print (gsearch1.best_params_)
print('best score')
print (gsearch1.best_score_)

		wheat_yield
		mean	min
District	wheat_variety
Banke	Banganga	51.673333	47.052
	Indian Hybrid 2285	38.859500	38.719
	NL 971	56.288250	49.360
	Unknown	48.200833	37.908
	Vijay	54.399500	53.343
Bardiya	Banganga	52.215875	49.540
	Local	48.159500	39.768
	NL 971	48.746857	46.184
	Tilottama	42.915500	42.361
	Unknown	45.468600	38.390
	Vijay	51.320375	40.525
Kailali	Aditya	48.239444	45.045
	Banganga	51.846667	44.290
	Local	48.487333	35.360
	NL 971	48.385000	37.670
	Vijay	48.782500	44.370
Kanchanpur	Banganga	55.053667	52.631
	Gautam	43.013500	37.003
	Local	50.484500	44.880
	NHS 1755	45.060000	40.440
	NL 971	45.857222	35.677
	Tilganga	50.620000	50.620
	Vijay	49.881000	38.920

		wheat_yield
		mean	min	max
District	Treat
Banke	D1	52.336600	47.052	57.303
	D2	49.431833	39.000	57.505
	D3	49.603667	37.908	60.985
Bardiya	D1	50.597000	46.184	56.418
	D2	47.107600	38.390	56.181
	D3	48.404067	40.346	57.723
Kailali	D1	49.673750	37.670	63.030
	D2	48.619375	35.360	58.530
	D3	47.573667	37.770	56.790
Kanchanpur	D1	48.156333	35.677	58.610
	D2	50.056583	37.003	62.500
	D3	47.920083	38.920	57.340

	Lat	Long	District	wheat_variety	Treat	wheat_yield	PH	OM	geometry
0	28.980	80.193	Kanchanpur	Banganga	D1	53.92	7.817228	1.688501	POINT (80.19300 28.98000)
1	28.980	80.193	Kanchanpur	Tilganga	D2	50.62	7.817228	1.688501	POINT (80.19300 28.98000)
2	28.980	80.193	Kanchanpur	NHS 1755	D3	40.44	7.817228	1.688501	POINT (80.19300 28.98000)
3	28.986	80.188	Kanchanpur	NL 971	D1	46.32	7.725938	1.644356	POINT (80.18800 28.98600)
4	28.986	80.188	Kanchanpur	Vijay	D2	62.50	7.725938	1.644356	POINT (80.18800 28.98600)

	Lat	Long	District	wheat_variety	Treat	wheat_yield	PH	OM	geometry
0	28.980	80.193	Kanchanpur	Banganga	D1	53.92	7.817228	1.688501	POINT (80.19300 28.98000)
1	28.980	80.193	Kanchanpur	Tilganga	D2	50.62	7.817228	1.688501	POINT (80.19300 28.98000)
2	28.980	80.193	Kanchanpur	NHS 1755	D3	40.44	7.817228	1.688501	POINT (80.19300 28.98000)
3	28.986	80.188	Kanchanpur	NL 971	D1	46.32	7.725938	1.644356	POINT (80.18800 28.98600)
4	28.986	80.188	Kanchanpur	Vijay	D2	62.50	7.725938	1.644356	POINT (80.18800 28.98600)

	B2	B3	B4	B8	B11	B12	ndvi	ndwi
0	0.246886	0.340987	0.397884	0.441628	0.281341	0.334612	0.052106	0.128595
1	0.246886	0.340987	0.397884	0.441628	0.281341	0.334612	0.052106	0.128595
2	0.246886	0.340987	0.397884	0.441628	0.281341	0.334612	0.052106	0.128595
3	0.223089	0.326231	0.384491	0.558959	0.269641	0.335698	0.184925	0.262913
4	0.223089	0.326231	0.384491	0.558959	0.269641	0.335698	0.184925	0.262913

	Lat	Long	District	wheat_variety	Treat	wheat_yield	PH	OM	geometry
0	28.98	80.193	Kanchanpur	Banganga	D1	53.92	7.817228	1.688501	POINT (80.19300 28.98000)
1	28.98	80.193	Kanchanpur	Tilganga	D2	50.62	7.817228	1.688501	POINT (80.19300 28.98000)
2	28.98	80.193	Kanchanpur	NHS 1755	D3	40.44	7.817228	1.688501	POINT (80.19300 28.98000)

	tavg	prec
0	0.598473	0.372607
1	0.598473	0.372607
2	0.598473	0.372607
3	0.538391	0.391679
4	0.538391	0.391679