<jupyter_start><jupyter_text>Figures 4 & 5. Validate CSDL against NASS county statistics Here we create Figures 4 & 5 from the paper, which summarize the agreement between CSDL and NASS statistics at the county level. CDL and NASS agreement is also shown for comparision.<jupyter_code># if needed, change these to the directories containing your data and the name of your data file
nass_dir = '../data/nass_counties/'
cdl_dir = '../data/cdl_counties/'
csdl_dir = '../data/csdl_counties/'<jupyter_output><empty_output><jupyter_text>Step 1. Load Datasets  1.1. NASS county-level statistics First, we load the NASS county-level statistics from 1999-2018.<jupyter_code>nass_corn = pd.read_csv(os.path.join(nass_dir, 'NASS_cropAreaCorn_1999to2018_raw.csv'))
nass_soy = pd.read_csv(os.path.join(nass_dir, 'NASS_cropAreaSoy_1999to2018_raw.csv'))

cols = ['state_fips_code', 'county_code', 'year', 'state_alpha', 
        'class_desc', 'short_desc','statisticcat_desc', 'commodity_desc',
        'util_practice_desc', 'Value']
nass_soy = nass_soy[cols]
nass_corn = nass_corn[cols]
nass = pd.concat([nass_corn, nass_soy])
print(nass.shape)
nass.head()
# Add the unique county FIPS code: stateFIPS+countyFIPS
nass['county_code'] = nass['county_code'].map(int).apply(lambda x: '{0:0>3}'.format(x))
nass['state_fips_code'] = nass['state_fips_code'].apply(lambda x: '{0:0>2}'.format(x))
nass['fips'] = (nass['state_fips_code'].map(str)+nass['county_code'].map(str)).map(int)

nass['commodity_desc'] = nass['commodity_desc'].str.title()
nass = nass.rename(columns={"commodity_desc":"crop", "state_alpha":"state", "Value":"Nass_Area_acres"})
nass['Nass_Area_acres'] = nass['Nass_Area_acres'].str.replace(',', '').astype('float')
nass["Nass_Area_m2"] = nass["Nass_Area_acres"] * 4046.86
nass = nass[['fips', 'year', 'state', 'crop', 'Nass_Area_m2']]
nass.head()<jupyter_output><empty_output><jupyter_text>1.2. Aggregated CDL at county level Next, we load CDL aggregated to the county level across the 13 states from 1999-2018. This data was exported in Google Earth Engine. A few steps of reorganizing the dataframe allows us to compare it to the NASS data.<jupyter_code>cdl = pd.DataFrame()
for filename in sorted(os.listdir(cdl_dir)):
    if (filename.endswith('.csv')):
        temp = pd.read_csv(os.path.join(cdl_dir, filename)).drop(['.geo','system:index'],axis=1)
        cdl = pd.concat([temp, cdl], sort=True)
cdl = cdl[cdl['masterid'] != 0] # drop dummy feature
print(cdl.shape)
cdl.head()
# compute CDL coverage by county
classes = list(set(cdl.columns.tolist()) - set(['Year', 'area_m2', 'masterid', 'COUNTYFP', 'STATEFP']))
other = list(set(classes) - set(['1', '5']))
cdl['cdl_coverage'] = cdl[classes].sum(axis=1)
cdl['Other'] = cdl[other].sum(axis=1)
cdl = cdl.drop(other, axis=1)
cdl.head()
maxcoverage = cdl[cdl['Year'] == 2018][['masterid', 'cdl_coverage']].rename(
    {'cdl_coverage': '2018cdl_coverage'}, axis=1)
cdl = cdl.merge(maxcoverage, on='masterid')
cdl['CDL_perccov'] = cdl['cdl_coverage'] / cdl['2018cdl_coverage'] * 100
cdl_key = {'1': "Corn", '5': "Soybeans", 'Year': "year"}  
cdl = cdl.rename(cdl_key, axis=1)
cdl.head()
crops = ["Corn", "Soybeans", "Other"]
cdl = pd.melt(cdl, id_vars=['masterid', 'year', 'COUNTYFP', 'STATEFP', 'CDL_perccov'], value_vars=crops, value_name='CDL_Area_m2')
cdl = cdl.rename(columns={"variable": "crop", "masterid": "fips"})
cdl.head()<jupyter_output><empty_output><jupyter_text>1.3. Aggregated CSDL at county level Lastly, we load CSDL aggregated to the county level, also data that was exported from GEE.<jupyter_code>csdl = pd.DataFrame()
for filename in sorted(os.listdir(csdl_dir)):
    if (filename.endswith('.csv')):
        temp = pd.read_csv(os.path.join(csdl_dir, filename))
        csdl = pd.concat([temp, csdl], sort=True)
csdl = csdl.dropna(subset=['masterid'])
csdl[['COUNTYFP', 'STATEFP', 'masterid', 'year']] = csdl[['COUNTYFP', 'STATEFP', 'masterid', 'year']].astype(int)
csdl = csdl[['0', '1', '5', 'COUNTYFP', 'STATEFP', 'masterid', 'year']]
print(csdl.shape)
csdl.head()
csdl = csdl.merge(maxcoverage, on='masterid')
csdl['CSDL_coverage'] = csdl[['0','1','5']].sum(axis=1)
csdl['CSDL_perccov'] = csdl['CSDL_coverage'] / csdl['2018cdl_coverage']*100
# note that CDL covers lakes, rivers; CSDL does not
csdl_key = {'0':"Other",'1': "Corn", '5': "Soybeans"} 
csdl = csdl.rename(csdl_key,axis=1)
csdl.head()
crops = ["Corn", "Soybeans", "Other"]
csdl = pd.melt(csdl, id_vars=['masterid', 'year', 'COUNTYFP', 'STATEFP', 'CSDL_perccov'], value_vars=crops, value_name='CSDL_Area_m2')
csdl = csdl.rename(columns={"variable": "crop", "masterid": "fips"})
csdl.head()<jupyter_output><empty_output><jupyter_text>1.4. Join CDL, CSDL and NASS<jupyter_code>df = nass.merge(cdl, on=['year', 'fips', 'crop'], how='left').merge(
    csdl, on=['year', 'fips', 'crop', 'COUNTYFP', 'STATEFP'], how='left')

# convert m2 to km2: 1e6
df['Nass_Area_km2'] = df['Nass_Area_m2'] / 1e6
df['CDL_Area_km2'] = df['CDL_Area_m2'] / 1e6
df['CSDL_Area_km2'] = df['CSDL_Area_m2'] / 1e6

print(df.shape)
df.head()<jupyter_output>(39986, 14)<jupyter_text>Step 2. Generate Figure 4: comparing  NASS-CSDL<jupyter_code>recent_years = np.arange(2008, 2019)
past_years = np.arange(1999, 2008)
df['year_group'] = '1999-2007'
df.loc[df['year'].isin(recent_years),'year_group'] = '2008-2018'

df_sub = df[df['CDL_perccov'] > 95]
print(df.shape)
print(df_sub.shape)
color_mapCSDL = { "Corn": 'orange',  "Soybeans": 'green'}
color_mapCDL = { "Corn": 'brown',  "Soybeans": 'darkseagreen'}
fsize = 15

lm = linear_model.LinearRegression()
bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.0)

fig, ax = plt.subplots(nrows=2, ncols=4,figsize=(16,8))

for i,(id,group) in enumerate(df_sub[df_sub['year_group']=='1999-2007'].groupby(['crop'])):
            
    ax[i,0].scatter(group['CDL_Area_km2'],group['Nass_Area_km2'], color = color_mapCDL[id], alpha=0.3, label=id)
    ax[i,1].scatter(group['CSDL_Area_km2'],group['Nass_Area_km2'], color = color_mapCSDL[id], alpha=0.3, label=id)

    group1= group.dropna(subset=['CDL_Area_km2','Nass_Area_km2'])
    R2 = r2_score(group1['CDL_Area_km2'].values,group1['Nass_Area_km2'].values)
    ax[i,0].text(0.05, 0.9, '$R^2$={0:.3f}'.format(R2), ha="left", va="center", size=fsize, bbox=bbox_props, transform=ax[i,0].transAxes)

    group2= group.dropna(subset=['CSDL_Area_km2','Nass_Area_km2'])
    R2 = r2_score(group2['CSDL_Area_km2'].values,group2['Nass_Area_km2'].values)
    ax[i,1].text(0.05, 0.9, '$R^2$={0:.3f}'.format(R2), ha="left", va="center", size=fsize, bbox=bbox_props, transform=ax[i,1].transAxes)

    lims = [0,2.6e+9/1e6] # km2
    # now plot both limits against each other
    ax[i,0].plot(lims, lims, 'k--', alpha=0.75, zorder=0)
    ax[i,1].plot(lims, lims, 'k--', alpha=0.75, zorder=0)
    ax[i,0].set_xlim(lims)
    ax[i,1].set_xlim(lims)
    ax[i,0].set_ylim(lims)
    ax[i,1].set_ylim(lims)  
    ax[i,0].set_aspect('equal')
    ax[i,1].set_aspect('equal')
    ax[0,0].set_xticklabels('', rotation=90)
    ax[0,1].set_xticklabels('', rotation=90)
    ax[i,1].set_yticklabels('', rotation=90)

for i,(id,group) in enumerate(df_sub[df_sub['year_group']=='2008-2018'].groupby(['crop'])):
       
    ax[i,2].scatter(group['CDL_Area_km2'],group['Nass_Area_km2'], color = color_mapCDL[id], alpha=0.3, label=id)
    ax[i,3].scatter(group['CSDL_Area_km2'],group['Nass_Area_km2'], color = color_mapCSDL[id], alpha=0.3, label=id)
    
    # set background to gray and alpha=0.2
    ax[i,2].set_facecolor((0.5019607843137255, 0.5019607843137255, 0.5019607843137255, 0.2))
    ax[i,3].set_facecolor((0.5019607843137255, 0.5019607843137255, 0.5019607843137255, 0.2))
    
    group1= group.dropna(subset=['CDL_Area_km2','Nass_Area_km2'])
    R2 = r2_score(group1['CDL_Area_km2'].values, group1['Nass_Area_km2'].values)
    ax[i,2].text(0.05, 0.9, '$R^2$={0:.3f}'.format(R2), ha="left", va="center", size=fsize, bbox=bbox_props, transform=ax[i,2].transAxes)

    group2= group.dropna(subset=['CSDL_Area_km2','Nass_Area_km2'])
    R2 = r2_score(group2['CSDL_Area_km2'].values, group2['Nass_Area_km2'].values)
    ax[i,3].text(0.05, 0.9, '$R^2$={0:.3f}'.format(R2), ha="left", va="center", size=fsize, bbox=bbox_props, transform=ax[i,3].transAxes)

    lims = [0,2.6e+9/1e6] # area in km2
    # now plot both limits against each other
    ax[i,2].plot(lims, lims, 'k--', alpha=0.75, zorder=0)
    ax[i,3].plot(lims, lims, 'k--', alpha=0.75, zorder=0)
    ax[i,2].set_xlim(lims)
    ax[i,3].set_xlim(lims)
    ax[i,2].set_ylim(lims)
    ax[i,3].set_ylim(lims)    
    ax[i,2].set_aspect('equal')
    ax[i,3].set_aspect('equal')
    ax[0,2].set_xticklabels('', rotation=90)
    ax[0,3].set_xticklabels('', rotation=90)
    ax[i,2].set_yticklabels('', rotation=90)
    ax[i,3].set_yticklabels('', rotation=90)

ax[0,0].set_ylabel('NASS county area [$km^2$]',fontsize=fsize)    
ax[1,0].set_ylabel('NASS county area [$km^2$]',fontsize=fsize)   
ax[1,0].set_xlabel('CDL county area [$km^2$]',fontsize=fsize)
ax[1,1].set_xlabel('CSDL county area [$km^2$]',fontsize=fsize)
ax[1,2].set_xlabel('CDL county area [$km^2$]',fontsize=fsize)
ax[1,3].set_xlabel('CSDL county area [$km^2$]',fontsize=fsize)
ax[0,0].set_title('1999-2007',fontsize=fsize) 
ax[0,1].set_title('1999-2007',fontsize=fsize) 
ax[0,2].set_title('2008-2018',fontsize=fsize) 
ax[0,3].set_title('2008-2018',fontsize=fsize) 

# Create the legend
legend_elements = [ Line2D([0], [0], color='k', linewidth=1, linestyle='--', label='1:1'),
                    Line2D([0], [0], marker='o', color='w', markerfacecolor='brown', markersize=15, label='Corn CDL'),
                    Line2D([0], [0], marker='o', color='w', markerfacecolor='orange', markersize=15, label='Corn CSDL'),
                    Line2D([0], [0], marker='o', color='w', markerfacecolor='darkseagreen', markersize=15, label='Soybeans CDL'),
                    Line2D([0], [0], marker='o', color='w', markerfacecolor='green', markersize=15, label='Soybeans CSDL'),
                    mpatches.Patch(facecolor='gray', edgecolor='gray', alpha = 0.2, label='NASS informed by CDL')]
fig.legend(handles=legend_elements,
           loc='lower center',  
           bbox_to_anchor=(0.5,-0.01),
           fontsize = 'x-large',
           ncol=7)

fig.tight_layout(rect=[0,0.05,1,1]) # legend on the bottom<jupyter_output><empty_output><jupyter_text>Step 3. Generate Figure 5: comparing NASS and CSDL over time by state<jupyter_code>def corrfun(df, col1, col2):
    
    df = df.dropna(subset=[col2])
    
    if df.shape[0] != 0:
        r2 = r2_score(df[col1].values, df[col2].values)
        mse = mean_squared_error(df[col1].values, df[col2].values)
    else:
        r2 = np.nan
        mse = np.nan
        
    return pd.Series({'R': r2, 'mse': mse, 'Ncounties': df.shape[0]}) # return R2 (coef of determination)
totArea = df.groupby(['year', 'state', 'crop'])['Nass_Area_m2', 'CDL_Area_m2', 'CSDL_Area_m2'].sum().reset_index()
print(totArea.shape) # 20years * 13states * 2commodity = 520 rows
corr_cdl = df[df['CDL_perccov'] > 90].groupby(['state','year','crop']).apply(corrfun,'Nass_Area_m2','CDL_Area_m2').reset_index().rename(
    {'R': 'R_NASS_CDL', 'mse': 'mse_NASS_CDL', 'Ncounties': 'Ncounties_CDL'}, axis=1)
corr_csdl = df.groupby(['state','year','crop']).apply(corrfun, 'Nass_Area_m2', 'CSDL_Area_m2').reset_index().rename(
    {'R': 'R_NASS_CSDL', 'mse': 'mse_NASS_CSDL', 'Ncounties': 'Ncounties_CSDL'}, axis=1)
corr = corr_cdl.merge(corr_csdl, on=['state','year','crop'], how='outer')
corr = corr.merge(totArea, on=['state','year','crop'], how='left')
print(corr.shape)

abbr_to_state = {'IL':'Illinois', 'IA':'Iowa', 'IN':'Indiana', 'NE':'Nebraska', 'ND':'North Dakota',
                 'SD':'South Dakota', 'MN':'Minnesota', 'WI':'Wisconsin', 'MI':'Michigan',
                 'KS':'Kansas','KY':'Kentucky', 'OH':'Ohio', 'MO':'Missouri'}
corr['state_abbrs'] = corr['state']
corr['state_name'] = corr['state'].replace(abbr_to_state)
corr.head()
state_geosorted = ['North Dakota', 'Minnesota', 'Wisconsin', 'Michigan',
 'South Dakota','Iowa', 'Illinois', 'Indiana',
 'Nebraska', 'Missouri', 'Kentucky', 'Ohio',
 'Kansas']

color_mapCSDL = { "Corn": 'orange',  "Soybeans": 'green'}
color_mapCDL = { "Corn": 'brown',  "Soybeans": 'darkseagreen'}
marker_map ={ "Corn": 'o',  "Soybeans": '^'}
fsize = 15
bbox_props = dict(boxstyle="round", fc="w", ec="0.5", alpha=0.0)

fig, ax = plt.subplots(nrows=4, ncols=4,figsize=(16,12))

for i, id in enumerate(state_geosorted):
    
    group = corr[corr['state_name']==id]

    thisAx = ax[int(np.floor(i/4)), i%4]
    thisAx.text(0.05, 0.1, id, ha="left", va="center", size=18, bbox=bbox_props, transform=thisAx.transAxes)
    
    ylabels = np.arange(0, 12, 2)/10
    thisAx.set_yticks(ylabels)
    thisAx.set_yticklabels('', rotation=0)
    
    xlabels = np.arange(1999, 2019, 2)
    thisAx.set_xticks(xlabels)
    thisAx.set_xticklabels('', rotation=90)
    thisAx.fill_between(np.arange(2008, 2019), 0, 1.1, alpha = 0.2, color='gray')

    crops = ['Soybeans','Corn']
    for l,id2 in enumerate(crops):
        group2 = group[group['crop']==id2]
        group2 = group2.sort_values('year')
        
        thisAx.plot(group2['year'],group2['R_NASS_CDL'],color=color_mapCDL[id2], alpha=1, marker=marker_map[id2])
        thisAx.plot(group2['year'],group2['R_NASS_CSDL'],color=color_mapCSDL[id2], alpha=1, marker=marker_map[id2])

    thisAx.set_xlim([1999,2018])
    thisAx.set_ylim([0.0, 1.1])
    thisAx.grid(True)

ax[0,0].set_ylabel('$R^{2}$',fontsize=fsize)
ax[1,0].set_ylabel('$R^{2}$',fontsize=fsize)
ax[2,0].set_ylabel('$R^{2}$',fontsize=fsize)
ax[3,0].set_ylabel('$R^{2}$',fontsize=fsize)
ax[0,0].set_yticklabels(ylabels, rotation=0,fontsize=fsize)
ax[1,0].set_yticklabels(ylabels, rotation=0,fontsize=fsize)
ax[2,0].set_yticklabels(ylabels, rotation=0,fontsize=fsize)
ax[3,0].set_yticklabels(ylabels, rotation=0,fontsize=fsize)
ax[3,0].set_xticklabels(xlabels, rotation=90,fontsize=fsize)
ax[-1,-1].axis('off')
ax[-1,-2].axis('off')
ax[-1,-3].axis('off')

# Create the legend manually
legend_elements = [ Line2D([0], [0], color='brown',  linewidth=3, marker='o', linestyle='-', label='Corn CDL'),
                    Line2D([0], [0], color='orange',  linewidth=3, marker='o', linestyle='-', label='Corn CSDL'),
                    Line2D([0], [0], color='darkseagreen',  linewidth=3, marker='^', linestyle='-', label='Soybeans CDL'),
                    Line2D([0], [0], color='green',  linewidth=3, marker='^', linestyle='-', label='Soybeans CSDL')]
fig.legend(handles=legend_elements,
           loc="lower right",   # Position of legend
           bbox_to_anchor=(0.6, 0.05),
           fontsize = 'xx-large')

bckground_patch = mpatches.Patch(color='gray', alpha = 0.2, label='Training years')
fig.legend(handles=[bckground_patch],
           loc="lower right",
           bbox_to_anchor=(0.8, 0.1),
           fontsize = 'xx-large')

fig.tight_layout()
ax[2,1].set_xticklabels(xlabels, rotation=90,fontsize=fsize)
ax[2,2].set_xticklabels(xlabels, rotation=90,fontsize=fsize)
ax[2,3].set_xticklabels(xlabels, rotation=90,fontsize=fsize)
fig.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20411799643634923, 'fraction_numerical': 0.05022107833432324, 'mean_word_length': 5.509450171821306, 'pattern_counts': {'":': 17, '<': 28, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '263637', 'license': 'MIT', 'n_tokens_mistral': 6710, 'n_tokens_neox': 6207, 'n_words': 988, 'path': 'validate_map/1_NASSvsCSDLvsCDL.ipynb', 'repo_name': 'LobellLab/csdl', 'size': '15153'}