# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
import numpy as np
import pandas as pd
import os
from sklearn import linear_model

import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
# %matplotlib inline
# -

# # Figure 7. Validate CSDL against ARMS crop rotation statistics

# if needed, change this to the directory containing your data and the name of your data file
arms_dir = '../data/arms_rotation/'
csdl_dir = '../data/csdl_rotation/'

# Agricultural Resource Management Survey (ARMS) crop rotation data was downloaded from the "Tailored reports: Crop Production Practices" found at https://www.ers.usda.gov/data-products/arms-farm-financial-and-crop-production-practices/.

# load ARMS data
arms = pd.read_csv(os.path.join(arms_dir, 'ARMS_priorCrop.csv'))
arms = arms.rename({'corn_estimate': 'arms_corn', 'soy_estimate': 'arms_soy', 'state':'STATE_ABBR'}, axis=1)
print(arms.shape)
arms.head()

# Crop rotation implied by CSDL maps was computed and exported in Google Earth Engine.

# load 2001, 2005 and 2010 data for CSDL
csdl_2001 = pd.read_csv(os.path.join(csdl_dir, '2001_CSDL_cropsYearPrior_corn.csv'))
csdl_2005 = pd.read_csv(os.path.join(csdl_dir, '2005_CSDL_cropsYearPrior_corn.csv'))
csdl_2010 = pd.read_csv(os.path.join(csdl_dir, '2010_CSDL_cropsYearPrior_corn.csv'))
csdl = pd.concat([csdl_2001, csdl_2005, csdl_2010], axis=0)
csdl = csdl.drop(['system:index', 'DRAWSEQ', 'SUB_REGION', 'groups', '.geo'], axis=1)
csdl = csdl.rename({'1': 'csdl_corn','5': 'csdl_soy', '0': 'csdl_other'}, axis=1)
print(csdl.shape) # 13states * 3years = 39
csdl.head()

# +
# merge datasets
df = arms.merge(csdl, on=['STATE_ABBR', 'year'], how='left')

df['arms_tot'] = df['arms_corn'] + df['arms_soy']
df['arms_corn_corn'] = df['arms_corn'] / df['arms_tot']
df['arms_soy_corn'] = df['arms_soy'] / df['arms_tot']

df['csdl_tot'] = df['csdl_corn'] + df['csdl_soy']
df['csdl_corn_corn'] = df['csdl_corn'] / df['csdl_tot']
df['csdl_soy_corn'] = df['csdl_soy'] / df['csdl_tot']

df['csdl_ratio_soy_corn'] = df['csdl_soy'] / df['csdl_corn']
df['arms_ratio_soy_corn'] = df['arms_soy'] / df['arms_corn']

df.head()

# +
fig, ax = plt.subplots(nrows=1, ncols=3,figsize=(15,5))
fsize = 15

NUM_COLORS = 13
cm = plt.get_cmap('tab20')
lm = linear_model.LinearRegression()
bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.0)

for i,(id,group) in enumerate(df.groupby('year')):
    
    ax[i].set_prop_cycle(color=[cm(1.*i/NUM_COLORS) for i in range(NUM_COLORS)])
    ax[i].set_title(id,fontsize=fsize)
    ax[i].set_xlabel('ARMS', fontsize=fsize)
    
    for id2,group2 in group.groupby('STATE_NAME'):
        ax[i].scatter(group2['arms_soy_corn'],group2['csdl_soy_corn'],label=id2)

    group= group.dropna(subset=['arms_soy_corn', 'csdl_soy_corn'])
    model = lm.fit(group[['arms_soy_corn']], group['csdl_soy_corn'])
    R2 = lm.score(group[['arms_soy_corn']], group['csdl_soy_corn'])

    predictions = lm.predict(group[['arms_soy_corn']])
    ax[i].plot(group[['arms_soy_corn']], predictions, linewidth=1)
    ax[i].text(0.05, 0.9, '$r^2$={0:.3f}'.format(R2), ha="left", va="center", 
               size=fsize, bbox=bbox_props, transform=ax[i].transAxes)

    # plot 1:1 line
    ax[i].plot([0,1], [0,1], 'k--', alpha=0.75, zorder=0)
    ax[i].set_aspect('equal')
    ax[i].set_xlim([0,1])
    ax[i].set_ylim([0,1])
    ax[i].grid(True)
    ax[i].legend().set_visible(False)
    
    ax[i].set_yticklabels('', rotation=0)

ax[0].set_ylabel('CSDL',fontsize=fsize)
ylabels = xlabels = np.arange(0,10, 2)/10
ax[0].set_yticklabels(ylabels,fontsize=fsize)
ax[0].set_xticklabels(xlabels,fontsize=fsize)
ax[1].set_xticklabels(xlabels,fontsize=fsize)
ax[2].set_xticklabels(xlabels,fontsize=fsize)


# Create the legend
colors=[cm(1.*i/NUM_COLORS) for i in range(NUM_COLORS)]
lines = [Line2D([0], [0], color=c, linewidth=3, linestyle='-') for c in colors]
mylabels = ['Fit','1:1','Illinois', 'Indiana', 'Iowa', 'Kansas', 'Kentucky', 'Michigan', 'Minnesota',
 'Missouri', 'Nebraska', 'North Dakota', 'Ohio', 'South Dakota', 'Wisconsin']

fig.legend(lines,     # The line objects
           labels=mylabels,   # The labels for each line
           loc="lower right",   # Position of legend
           bbox_to_anchor=(0.95, 0.03),
           fontsize='x-large')

fig.tight_layout(rect=[0,0,0.8,1])
fig.show()

StarcoderdataJupyterScriptsDedupFiltered