<jupyter_start><jupyter_text>Brazil soybeans crop model using Gro APIIn this notebook, we walk through the process of building a simple "crop model" using data in Gro. The accompanying Gro web app display for [Brazil Soybeans](https://app.gro-intelligence.com/displays/25894) provides an overview of the data series we will use. The crop model we will create is a very basic  model that uses crop-production-weighted NDVI to forecast  yield. PreliminaryFirst we construct a `CropModel` object. We also save the entity ids for the `item` [soybeans](https://app.gro-intelligence.com/dictionary/items/270), the `region` [Brazil](https://app.gro-intelligence.com/dictionary/regions/1029), the `metric` [yield](https://app.gro-intelligence.com/dictionary/metrics/170037) and yield units (sources report yield in different units, we specify tonnes per hectare so they'll be converted to a consistent unit for model training). We will use these throughout the notebook to simplify the retrieval, storage and manipulation of data.<jupyter_code>import os
from groclient.crop_model import CropModel

model = CropModel('api.gro-intelligence.com', os.environ['GROAPI_TOKEN'])
country_id = model.search_for_entity('regions', "brazil")
crop_id =  model.search_for_entity('items', "soybeans")
yield_id = model.search_for_entity('metrics', "yield mass/area")
TONNES_PER_HECTARE_UNIT_ID = 61
# Quiet logs for demonstration purposes:
import warnings
warnings.filterwarnings('ignore')<jupyter_output><empty_output><jupyter_text>Historical yieldsTo get data series for national level soybeans yield for Brazil, we set the entities (item, metric and region), and get the available data series.  Note that there are multiple data series in this case because there are a few different sources that report the Brazil soybeans yields: [FAO](https://app.gro-intelligence.com/dictionary/sources/2), [PS&D](https://app.gro-intelligence.com/dictionary/sources/14), [IGC](https://app.gro-intelligence.com/dictionary/sources/19), and [IBGE](https://app.gro-intelligence.com/dictionary/sources/114).<jupyter_code>yield_entities = {}
yield_entities['item_id']  = crop_id
yield_entities['region_id'] = country_id
yield_entities['metric_id'] = yield_id

data_series_list = model.get_data_series(**yield_entities)
print("There are {} data series for {}.".format(len(data_series_list), yield_entities))<jupyter_output>There are 6 data series for {'item_id': 270, 'region_id': 1029, 'metric_id': 170037}.<jupyter_text>Each source covers a slightly different time range. Here we simply take all of them as training data, to have maximum coverage. We exclude source 32, the [Gro Yield Model](https://app.gro-intelligence.com/dictionary/sources/32), since that is itself a forecast from Gro's more sophisticated automated yield models, as discussed in the conclusion below.<jupyter_code>GRO_YM_SOURCE_ID = 32
for data_series in data_series_list:
    if data_series['source_id'] == GRO_YM_SOURCE_ID:
        continue
    data_series['unit_id'] = TONNES_PER_HECTARE_UNIT_ID
    print("source_id {}: {} to {}".format(
        data_series['source_id'], data_series['start_date'], data_series['end_date']))
    model.add_single_data_series(data_series)
%matplotlib inline
import pandas
import matplotlib.pyplot as plt

df = model.get_df()
yield_df = df.loc[(df.metric_id == yield_entities['metric_id']) & \
                  (df.region_id == yield_entities['region_id'])]
yield_df.end_date = pandas.to_datetime(yield_df.end_date)
yield_df.set_index('end_date')
fig, axes = plt.subplots(1, 1)
for source_id, group in yield_df.groupby('source_id'):
    group.plot(x='end_date', y='value', ax=axes, x_compat=True,
               label=model.lookup('sources', source_id)['name'])
plt.ylabel("Yield (t/ha)")
plt.show()<jupyter_output><empty_output><jupyter_text>Production "weight" of provincesWe compute the "weight" of each province, based on the full history of production quantity per province. The weights are normalized to add up to 1.0. A weight of 0.2 means we expect that province to account for 20% of the country's production.<jupyter_code>provinces = model.get_provinces("brazil")
weights = model.compute_weights("soybeans", "Production Quantity mass", provinces)
max_weight = max(filter(lambda x: x >0, weights))
print ('Total weight: {}'.format(sum(filter(lambda x: x >0, weights))))
print ('Maximum weight: {}'.format(max_weight))
import geopandas
from shapely.geometry import shape

gdf = geopandas.GeoDataFrame([{
       'region_id': province['id'], 
       'geometry': shape(model.get_geojson(province['id'])['geometries'][0]),
       'production_weight': weight if weight > 0 else 0
    } for (province, weight) in zip(provinces, weights)])
gdf.set_index('region_id')

fig, axes = plt.subplots(1, 1)
gdf.plot(column='production_weight', 
         cmap='Blues', ax=axes, vmin=0, vmax=max_weight, legend=True)
axes.set_title('Crop Production Weight')
axes.xaxis.set_visible(False)
axes.yaxis.set_visible(False)
plt.show()<jupyter_output><empty_output><jupyter_text>NDVI crop-weighted, province levelThe main signal we will use to model crop yields is [NDVI](https://app.gro-intelligence.com/dictionary/items/321), which represents vegetation biomass per pixel, and thus is a good physical proxy for yield (production mass per unit of area). NDVI represents the biomass  of all vegetation, not just soybeans. To make it a more accurate signal, we compute the crop-weighted version of NDVI, using [compute_crop_weighted_series](https://developers.gro-intelligence.com/api.htmlapi.client.crop_model.CropModel.compute_crop_weighted_series).<jupyter_code>soy_ndvi_df = model.compute_crop_weighted_series(
    'soybeans', 'Production Quantity mass',
    'Vegetation NDVI', 'Vegetation Indices index',
    provinces)<jupyter_output><empty_output><jupyter_text>Now let's look at NDVI by province, on a particular date. First, we select NDVI for all provinces from the model's dataframe, and join it with the geopandas dataframe that has the provinces geometry.   Second, we select the crop weighted NDVI for all provinces from the `soy_ndvi_df` dataframe and do the same.<jupyter_code>entities = {}
entities['item_id'] =  model.search_for_entity('items', "Vegetation NDVI")
entities['metric_id'] = model.search_for_entity('metrics', "Vegetation Indices index")

EVAL_DATE = u'2019-02-25T00:00:00.000Z'
df = model.get_df()
ndvi_df = df.loc[(df.metric_id == entities['metric_id']) & \
                 (df.end_date == EVAL_DATE),:]
gdf1 = gdf.join(ndvi_df.set_index('region_id'), on='region_id')

cwndvi_df = soy_ndvi_df.loc[(soy_ndvi_df.metric_id == entities['metric_id']) & \
                            (soy_ndvi_df.end_date == EVAL_DATE), :]
gdf2 = gdf.join(cwndvi_df.set_index('region_id'), on='region_id')

fig, axes = plt.subplots(1, 2)
gdf1.plot(column='value', 
          cmap='Greens', ax=axes[0], vmin=0.0, vmax=1.0)
axes[0].set_title('NDVI')
axes[0].xaxis.set_visible(False)
axes[0].yaxis.set_visible(False)
gdf2.plot(column='value', 
          cmap='Greens', ax=axes[1], vmin=0.0, vmax=max_weight)
axes[1].set_title('Crop-weighted NDVI')
axes[1].xaxis.set_visible(False)
axes[1].yaxis.set_visible(False)
plt.show()<jupyter_output><empty_output><jupyter_text>NDVI, country levelNDVI is available as a data series at all region levels, so we just add it to the model at the country level.<jupyter_code>entities['region_id'] = country_id
for data_series in model.get_data_series(**entities):
    model.add_single_data_series(data_series)
df = model.get_df()
df.end_date = pandas.to_datetime(df.end_date, utc=True)
country_ndvi_df = df.loc[(df.metric_id == entities['metric_id']) & \
                         (df.region_id == entities['region_id'])]<jupyter_output><empty_output><jupyter_text>For crop weighted NDVI, we need to aggregate it up from provinces. Since the crop-weighted NDVI series was built with weights normalized to 1.0, the sum of province level values is the national level value, for each date.<jupyter_code>soy_ndvi_df.end_date = pandas.to_datetime(soy_ndvi_df.end_date, utc=True)
soy_ndvi_df.set_index('end_date')
country_ndvi_df.set_index('end_date')

country_ndvi_df = country_ndvi_df.merge(pandas.DataFrame([{
    'end_date': end_date,
    'crop_weighted_value': value
} for (end_date, value) in zip(soy_ndvi_df.end_date, 
                               soy_ndvi_df.groupby(['end_date']).sum().value)]))

fig, axes = plt.subplots(1, 1)
country_ndvi_df.plot(x='end_date', y=['crop_weighted_value', 'value'], ax=axes)
plt.ylabel('NDVI')
plt.show()<jupyter_output><empty_output><jupyter_text>Relationship between  yield and crop-weighted NDVIOn any day during the season, we can look at crop weighted NDVI as a rough aproximation of the amount of soybeans being grown per unit of area. Let's see if this would work as a predictor of the final yield of the crop for the current year. For this example, we pick February 25th of each year, which, as we can see from the crop calendar for this crop on our [display](https://app.gro-intelligence.com/displays/25894), is the end of the growth period, and just before harvesting starts.  This date is late enough in the season (end of growth) that the weighted NDVI should represent the final amount of the crop that will be produced, but it is early enough (before harvest) for the forecast to be useful.  Let's look at a scatter plot of historical values of  crop weighted NDVI on Feb 25 and final yields for the same year.<jupyter_code>PREDICTION_DATE_MONTH = 2
PREDICTION_DATE_DAY = 25
tmp_df = country_ndvi_df.loc[(country_ndvi_df.metric_id == entities['metric_id']) & \
                             (country_ndvi_df.region_id == country_id) & \
                             (country_ndvi_df.end_date.map(lambda x: x.month) == PREDICTION_DATE_MONTH) &\
                             (country_ndvi_df.end_date.map(lambda x: x.day) == PREDICTION_DATE_DAY),:]
ym_df = pandas.DataFrame([{'year': end_date.year, 'yield': value}
                          for (end_date, value) in zip(yield_df.end_date, 
                                                   yield_df.value)])
cw_df = pandas.DataFrame([{'year': end_date.year, 'crop_weighted_ndvi': value}
                          for (end_date, value) in zip(tmp_df.end_date,
                                                       tmp_df.crop_weighted_value)])
ym_df = ym_df.merge(cw_df)
ym_df.plot.scatter(x='crop_weighted_ndvi', y = 'yield')
plt.show()<jupyter_output><empty_output><jupyter_text>Indeed there seems to be a relationship between these two series, which means weighted NDVI on Feb 25 can be a reasonable predictor of final yield. Linear modelAn obvious next step is to do a linear regression, where we model `y = a * x + b` where x is the crop weighted NDVI.<jupyter_code>from sklearn import linear_model

lm = linear_model.LinearRegression()
X = ym_df.loc[:,['crop_weighted_ndvi']]
lm.fit(X, ym_df['yield'])
lm_df =pandas.DataFrame([{'year': year, 'lm_yield': value}
                         for (year, value) in zip(ym_df.year, lm.predict(X))])
plt.scatter(ym_df['yield'], (lm_df['lm_yield']))
plt.xlabel("Actual yield")
plt.ylabel("Predicted yield (t/ha)")
plt.plot([min(ym_df['yield']), max(ym_df['yield'])], 
         [min(ym_df['yield']), max(ym_df['yield'])], color='red')
plt.show()<jupyter_output><empty_output><jupyter_text>That is not a great fit, especially at higher yields, which is not surprising as we're using crop weights based on production rather than planted area.  Ideally to get the national yield, we would weight regional NDVI by  area.  Since production equals area multiplied by yield,  weighting by production instead of area will tend to over-wheight the high yield regions. Thus the model will be less accurate in years when the difference between regions is relatively large.  Quadratic+trend modelA polynomial of degree 2 should help us capture this second order effect, and provide a better fit. So we try to fit  `y = a * x^2 + b * x + c` where x is the  crop weighted NDVI. Further, as there appears to be a longer term trend of improvement in yields, we add the year as another dependent variable to capture the long term trend. So the model becomes `y = a * x^2 + b * x + c + d * t`, where t is time (year).<jupyter_code>X2 = pandas.DataFrame(X.crop_weighted_ndvi.map(lambda x: x*x))
X2.rename(columns={'crop_weighted_ndvi': 'crop_weighted_ndvi2'}, inplace=True)
X2 = pandas.DataFrame([X2.crop_weighted_ndvi2, X.crop_weighted_ndvi, ym_df.year]).transpose()
qm = linear_model.LinearRegression()
qm.fit(X2, ym_df['yield'])
qm_df =pandas.DataFrame([{'year': year, 'qm_yield': value}
                         for (year, value) in zip(ym_df.year, qm.predict(X2))])
plt.scatter(ym_df['yield'], (qm_df['qm_yield']))
plt.plot([min(ym_df['yield']), max(ym_df['yield'])], 
         [min(ym_df['yield']), max(ym_df['yield'])], color='red')
plt.xlabel("Actual yield (t/ha)")
plt.ylabel("Predicted yield (t/ha)")
plt.show()<jupyter_output><empty_output><jupyter_text>ResultsThe second model seems to provide a better looking fit. Looking at the standard regression statistics to confirm that:<jupyter_code>from sklearn.metrics import mean_squared_error, r2_score
print('''Linear Model:
      Mean Squared Error: {}
      R2: {}'''.format(mean_squared_error(ym_df['yield'], lm.predict(X)),
                      r2_score(ym_df['yield'], lm.predict(X))))
print('''Quadtratic Model:
      Mean Squared Error: {}
      R2: {}'''.format(mean_squared_error(ym_df['yield'], qm.predict(X2)),
                          r2_score(ym_df['yield'], qm.predict(X2))))<jupyter_output>Linear Model:
      Mean Squared Error: 0.0985408856247
      R2: 0.037707087065
Quadtratic Model:
      Mean Squared Error: 0.0298944694076
      R2: 0.708068018016<jupyter_text>Indeed the quadratic+trend model- has a root-mean-square error of less than 0.03 t/ha, i.e. under 1% - explains about 71% of the variance in yieldsTo get more perspective, let's see how well the models' fit reality as a function of time. As we can see in the following figure, the quadratic model does indeed outperform the linear model, especially in the high yield years.<jupyter_code>backtest_df = ym_df.merge(lm_df).merge(qm_df)
backtest_df.year = pandas.to_datetime(backtest_df.year, format='%Y')
fig, axes = plt.subplots(1, 1)
backtest_df.plot(x='year', y = ['yield', 'lm_yield', 'qm_yield'], ax=axes)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11664700575948928, 'fraction_numerical': 0.014155853167719104, 'mean_word_length': 4.811290322580645, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 54, 'https://': 12, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '202128', 'license': 'MIT', 'n_tokens_mistral': 4953, 'n_tokens_neox': 4621, 'n_words': 1332, 'path': 'api/client/samples/crop_models/brazil_soybeans.ipynb', 'repo_name': 'eric-gro/api-client', 'size': '14411'}