import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

import netCDF4 as nc
import pandas as pd
import os
import csv
from glob import glob
import xarray as xr
import matplotlib.ticker as ticker 

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib.colors import LinearSegmentedColormap, Normalize
from mpl_toolkits.axes_grid1 import make_axes_locatable

import geopandas as gpd
import matplotlib.colors as mcolors

# Step 1: Data Preparation
# Load your dataset into a pandas DataFrame
os.chdir('/Users/chenchenren/postdoc/paper/2N and water-US/Figure 3 Yiled under warming scenario/')

# Load data from an Excel file
excel_file = "scenario_management.xlsx"
sheet_name = "soybean"
data = pd.read_excel(excel_file, sheet_name=sheet_name)

data1=data[data['sample'] == 1].copy()
data1=data1.sort_values(by=['geoid', 'scenario'])
print(data1)

data2015=data1[data1['scenario'] == 2015].copy()
print(data2015)


data1['tmp'] = data1.groupby(['geoid','scenario'])['tmp'].transform('mean')
data1['pre'] = data1.groupby(['geoid','scenario'])['pre'].transform('mean')
data1['tmp2015'] = data1.groupby(['geoid','scenario'])['tmp2015'].transform('mean')
data1['pre2015'] = data1.groupby(['geoid','scenario'])['pre2015'].transform('mean')
data1=data1.drop_duplicates(subset=['geoid','scenario'])
data1=data1.drop(['model','lnyield'], axis=1)
print(data1)

       geoid  year       tmp       pre       har  irrigation1  lnyield_obs  \
26081   1001  2015  25.44350  5.702415   779.625     0.377088     7.509678   
26084   1001  2015  24.71859  5.866428   779.625     0.377088     7.509678   
26086   1001  2015  25.15585  5.960842   779.625     0.377088     7.509678   
26087   1001  2015  24.83175  7.541522   779.625     0.377088     7.509678   
26089   1001  2015  25.34826  6.200767   779.625     0.377088     7.509678   
...      ...   ...       ...       ...       ...          ...          ...   
43802  55141  2015  15.16129  7.028600  6700.507     0.000000     7.921128   
43804  55141  2015  16.77098  5.360460  6700.507     0.000000     7.921128   
43805  55141  2015  16.46397  6.299626  6700.507     0.000000     7.921128   
43807  55141  2015  16.07257  6.365017  6700.507     0.000000     7.921128   
43809  55141  2015  18.05135  5.488240  6700.507     0.000000     7.921128   

          lnfer  lnyield  group  sample  scenario   pre2015   model   tmp2015  
26081  2.457656      NaN      1       1       1.5  4.174539    gfdl  26.21985  
26084  2.457656      NaN      1       1       1.5  6.327019    ipsl  24.13167  
26086  2.457656      NaN      1       1       1.5  4.268701     mpi  24.93584  
26087  2.457656      NaN      1       1       1.5  8.572182     mri  25.04051  
26089  2.457656      NaN      1       1       1.5  9.683448  ukesm1  24.72551  
...         ...      ...    ...     ...       ...       ...     ...       ...  
43802  2.556201      NaN      2       1    2015.0  7.028600    gfdl  15.16129  
43804  2.556201      NaN      2       1    2015.0  5.360460    ipsl  16.77098  
43805  2.556201      NaN      2       1    2015.0  6.299626     mpi  16.46397  
43807  2.556201      NaN      2       1    2015.0  6.365017     mri  16.07257  
43809  2.556201      NaN      2       1    2015.0  5.488240  ukesm1  18.05135  

[26595 rows x 15 columns]
       geoid  year       tmp       pre       har  irrigation1  lnyield_obs  \
17218   1001  2015  24.93584  4.268701   779.625     0.377088     7.509678   
26082   1001  2015  26.21985  4.174539   779.625     0.377088     7.509678   
26083   1001  2015  24.13167  6.327019   779.625     0.377088     7.509678   
26088   1001  2015  25.04051  8.572182   779.625     0.377088     7.509678   
26090   1001  2015  24.72551  9.683448   779.625     0.377088     7.509678   
...      ...   ...       ...       ...       ...          ...          ...   
43802  55141  2015  15.16129  7.028600  6700.507     0.000000     7.921128   
43804  55141  2015  16.77098  5.360460  6700.507     0.000000     7.921128   
43805  55141  2015  16.46397  6.299626  6700.507     0.000000     7.921128   
43807  55141  2015  16.07257  6.365017  6700.507     0.000000     7.921128   
43809  55141  2015  18.05135  5.488240  6700.507     0.000000     7.921128   

          lnfer  lnyield  group  sample  scenario   pre2015   model   tmp2015  
17218  2.457656      NaN      1       1    2015.0  4.268701     mpi  24.93584  
26082  2.457656      NaN      1       1    2015.0  4.174539    gfdl  26.21985  
26083  2.457656      NaN      1       1    2015.0  6.327019    ipsl  24.13167  
26088  2.457656      NaN      1       1    2015.0  8.572182     mri  25.04051  
26090  2.457656      NaN      1       1    2015.0  9.683448  ukesm1  24.72551  
...         ...      ...    ...     ...       ...       ...     ...       ...  
43802  2.556201      NaN      2       1    2015.0  7.028600    gfdl  15.16129  
43804  2.556201      NaN      2       1    2015.0  5.360460    ipsl  16.77098  
43805  2.556201      NaN      2       1    2015.0  6.299626     mpi  16.46397  
43807  2.556201      NaN      2       1    2015.0  6.365017     mri  16.07257  
43809  2.556201      NaN      2       1    2015.0  5.488240  ukesm1  18.05135  

[8865 rows x 15 columns]
       geoid  year        tmp       pre        har  irrigation1  lnyield_obs  \
26081   1001  2015  25.099590  6.254395    779.625     0.377088     7.509678   
17216   1001  2015  26.784940  6.208599    779.625     0.377088     7.509678   
17218   1001  2015  25.010676  6.605178    779.625     0.377088     7.509678   
17222   1003  2015  26.261318  9.274925   6358.433     0.738441     7.927428   
17221   1003  2015  27.808318  8.926058   6358.433     0.738441     7.927428   
...      ...   ...        ...       ...        ...          ...          ...   
26072  55139  2015  20.012748  5.680057  14541.420     0.000000     8.028848   
26071  55139  2015  17.287318  5.693542  14541.420     0.000000     8.028848   
26077  55141  2015  17.398562  5.719637   6700.507     0.000000     7.921128   
26076  55141  2015  19.304302  6.116161   6700.507     0.000000     7.921128   
43802  55141  2015  16.504032  6.108389   6700.507     0.000000     7.921128   

          lnfer  group  sample  scenario   pre2015    tmp2015  
26081  2.457656      1       1       1.5  6.605178  25.010676  
17216  2.457656      1       1       3.0  6.605178  25.010676  
17218  2.457656      1       1    2015.0  6.605178  25.010676  
17222  1.292553      1       1       1.5  9.548671  26.127096  
17221  1.292553      1       1       3.0  9.548671  26.127096  
...         ...    ...     ...       ...       ...        ...  
26072  1.255133      2       1       3.0  5.693542  17.287318  
26071  1.255133      2       1    2015.0  5.693542  17.287318  
26077  2.556201      2       1       1.5  6.108389  16.504032  
26076  2.556201      2       1       3.0  6.108389  16.504032  
43802  2.556201      2       1    2015.0  6.108389  16.504032  

[5319 rows x 13 columns]

# Load data from an Excel file
excel_file = "scenario_management.xlsx"
sheet_name = "soybean"
data = pd.read_excel(excel_file, sheet_name=sheet_name)

# Create interaction terms
data['tmp_tmp_interaction'] = data['tmp'] ** 2
data['pre_pre_interaction'] = data['pre'] ** 2
data['irrigation12'] = data['irrigation1'] ** 2
data['lnfer_irrigation1_tmp_interaction'] = data['lnfer'] * data['irrigation1'] * data['tmp']
data['lnfer_irrigation12_tmp_tmp_interaction'] = data['lnfer'] * data['irrigation12'] * data['tmp'] ** 2
data['lnfer_irrigation1_pre_interaction'] = data['lnfer'] * data['irrigation1'] * data['pre']
data['lnfer_irrigation12_pre_pre_interaction'] = data['lnfer'] * data['irrigation12'] * data['pre'] ** 2

data['lnfer_tmp_interaction'] = data['lnfer'] * data['tmp']
data['lnfer_tmp_tmp_interaction'] = data['lnfer'] * data['tmp'] ** 2
data['lnfer_pre_interaction'] = data['lnfer'] * data['pre']
data['lnfer_pre_pre_interaction'] = data['lnfer'] * data['pre'] ** 2

# Define predictor variables and response variable
predictor_variables1 = ['lnfer', 'tmp', 'irrigation1', 'irrigation12', 'tmp_tmp_interaction',
                        'pre', 'pre_pre_interaction', 'lnfer_irrigation1_pre_interaction', 'lnfer_irrigation12_pre_pre_interaction',
                        'lnfer_irrigation1_tmp_interaction', 'lnfer_irrigation12_tmp_tmp_interaction']

predictor_variables2 = ['lnfer', 'tmp', 'tmp_tmp_interaction', 'irrigation1', 'irrigation12', 'pre', 'pre_pre_interaction',
                        'lnfer_tmp_interaction', 'lnfer_tmp_tmp_interaction',
                        'lnfer_pre_interaction', 'lnfer_pre_pre_interaction']

response_variable = 'lnyield'

# Function to perform regression analysis
def perform_regression(group_data, predictor_variables):
    X = sm.add_constant(group_data[predictor_variables])
    y = group_data[response_variable]
    return sm.OLS(y, X).fit()

# Create a list to store data rows
data_rows = []

sample_data = data[data['sample'] == 0].copy()
print(sample_data)
project_data = data1
print(project_data)

# Iterate through 'results' and extract values
for group in sample_data['group'].unique():
    group_data = sample_data[sample_data['group'] == group].copy()
    predictor_variables = predictor_variables1 if group == 1 else predictor_variables2
    reg = perform_regression(group_data, predictor_variables)

    for group_proj in project_data['group'].unique():
        if group == group_proj:
            # Calculate y_temp_base based on regression parameters
            tmp_values = project_data[project_data['group'] == group_proj]['tmp'].values
            pre_values = project_data[project_data['group'] == group_proj]['pre'].values
            tmp2015 = project_data[project_data['group'] == group_proj]['tmp2015'].values
            pre2015 = project_data[project_data['group'] == group_proj]['pre2015'].values
            lnfer_values = project_data[project_data['group'] == group_proj]['lnfer'].values
            irrigation1_values = project_data[project_data['group'] == group_proj]['irrigation1'].values
            
            if group_proj == 1:
                print(reg.params['tmp_tmp_interaction'])
                print(reg.params['const'])
                y_temp_base = (
                    reg.params['tmp_tmp_interaction'] * tmp_values**2 +
                    reg.params['tmp'] * tmp_values +
                    reg.params['pre_pre_interaction'] * pre_values**2 +
                    reg.params['pre'] * pre_values +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['irrigation1'] * irrigation1_values +
                    reg.params['irrigation12'] * irrigation1_values**2 +
                    reg.params['lnfer_irrigation1_tmp_interaction'] * irrigation1_values * tmp_values * lnfer_values +
                    reg.params['lnfer_irrigation12_tmp_tmp_interaction'] * irrigation1_values**2 * lnfer_values * tmp_values**2 +
                    reg.params['lnfer_irrigation1_pre_interaction'] * irrigation1_values * pre_values * lnfer_values +
                    reg.params['lnfer_irrigation12_pre_pre_interaction'] * irrigation1_values**2 * lnfer_values * pre_values**2 +
                    reg.params['const']
                )
                y_tem = ( # temperature induced
                    reg.params['tmp_tmp_interaction'] * tmp_values**2 +
                    reg.params['tmp'] * tmp_values +
                    reg.params['pre_pre_interaction'] * pre2015**2 +
                    reg.params['pre'] * pre2015 +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['irrigation1'] * irrigation1_values +
                    reg.params['irrigation12'] * irrigation1_values**2 +
                    reg.params['lnfer_irrigation1_tmp_interaction'] * irrigation1_values * tmp_values * lnfer_values +
                    reg.params['lnfer_irrigation12_tmp_tmp_interaction'] * irrigation1_values**2 * lnfer_values * tmp_values**2 +
                    reg.params['lnfer_irrigation1_pre_interaction'] * irrigation1_values * pre2015 * lnfer_values +
                    reg.params['lnfer_irrigation12_pre_pre_interaction'] * irrigation1_values**2 * lnfer_values * pre2015**2 +
                    reg.params['const']
                )
                y_pre = ( # pre induced
                    reg.params['tmp_tmp_interaction'] * tmp2015**2 +
                    reg.params['tmp'] * tmp2015 +
                    reg.params['pre_pre_interaction'] * pre_values**2 +
                    reg.params['pre'] * pre_values +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['irrigation1'] * irrigation1_values +
                    reg.params['irrigation12'] * irrigation1_values**2 +
                    reg.params['lnfer_irrigation1_tmp_interaction'] * irrigation1_values * tmp2015 * lnfer_values +
                    reg.params['lnfer_irrigation12_tmp_tmp_interaction'] * irrigation1_values**2 * lnfer_values * tmp2015**2 +
                    reg.params['lnfer_irrigation1_pre_interaction'] * irrigation1_values * pre_values * lnfer_values +
                    reg.params['lnfer_irrigation12_pre_pre_interaction'] * irrigation1_values**2 * lnfer_values * pre_values**2 +
                    reg.params['const']
                )
                print(len(y_temp_base))
                print(len(y_tem))
                print(len(y_pre))
            else:
                print(reg.params['tmp_tmp_interaction'])
                print(reg.params['const'])
                y_temp_base = (
                    reg.params['tmp_tmp_interaction'] * tmp_values**2 +
                    reg.params['tmp'] * tmp_values +
                    reg.params['pre_pre_interaction'] * pre_values**2 +
                    reg.params['pre'] * pre_values +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['lnfer_tmp_interaction'] * tmp_values * lnfer_values +
                    reg.params['lnfer_tmp_tmp_interaction'] * lnfer_values * tmp_values**2 +
                    reg.params['lnfer_pre_interaction'] * pre_values * lnfer_values +
                    reg.params['lnfer_pre_pre_interaction'] * lnfer_values * pre_values**2 +
                    reg.params['const']
                )
                y_tem = ( # tem induced
                    reg.params['tmp_tmp_interaction'] * tmp_values**2 +
                    reg.params['tmp'] * tmp_values +
                    reg.params['pre_pre_interaction'] * pre2015**2 +
                    reg.params['pre'] * pre2015 +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['lnfer_tmp_interaction'] * tmp_values * lnfer_values +
                    reg.params['lnfer_tmp_tmp_interaction'] * lnfer_values * tmp_values**2 +
                    reg.params['lnfer_pre_interaction'] * pre2015 * lnfer_values +
                    reg.params['lnfer_pre_pre_interaction'] * lnfer_values * pre2015**2 +
                    reg.params['const']
                )
                y_pre = ( # pre induced
                    reg.params['tmp_tmp_interaction'] * tmp2015**2 +
                    reg.params['tmp'] * tmp2015 +
                    reg.params['pre_pre_interaction'] * pre_values**2 +
                    reg.params['pre'] * pre_values +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['lnfer_tmp_interaction'] * tmp2015 * lnfer_values +
                    reg.params['lnfer_tmp_tmp_interaction'] * lnfer_values * tmp2015**2 +
                    reg.params['lnfer_pre_interaction'] * pre_values * lnfer_values +
                    reg.params['lnfer_pre_pre_interaction'] * lnfer_values * pre_values**2 +
                    reg.params['const']
                )
                print(len(y_temp_base))
                print(len(y_tem))
                print(len(y_pre))
            
            # Extract 'geoid', 'scenario', and 'model' from the original data
            geoid_values = project_data[project_data['group'] == group_proj]['geoid'].values
            scenario_values = project_data[project_data['group'] == group_proj]['scenario'].values
            lnyield_obs_values= project_data[project_data['group'] == group_proj]['lnyield_obs'].values
            har_values= project_data[project_data['group'] == group_proj]['har'].values
            
            for tmp, pre, lnfer, irrigation1,y_temp_base_val, y_tem_val, y_pre_val, geoid, scenario, lnyield_obs,har in zip(
                tmp_values, pre_values, lnfer_values, irrigation1_values,
                 y_temp_base, y_tem, y_pre, geoid_values, scenario_values, lnyield_obs_values,har_values
            ):
                data_rows.append([tmp, pre, lnfer, irrigation1, y_temp_base_val, y_tem_val, y_pre_val, geoid, scenario, 
                                  lnyield_obs,har])

# Create a DataFrame from the data
results_base = pd.DataFrame(data_rows, columns=['tmp', 'pre', 'lnfer', 'irrigation1',
                                                'y_temp_base', 'y_tem', 'y_pre', 'geoid', 'scenario', 
                                                'lnyield_obs', 'har'])

results_base_t = results_base[(results_base['scenario'] == 1.5) | (results_base['scenario'] == 3)]

# Rename the 'y_temp_base' column to 'y_temp_base2015'
results_base_t.rename(columns={'y_temp_base': 'y_temp_base_t'}, inplace=True)
results_base_t = results_base_t[['geoid',  'scenario', 'y_temp_base_t', 'y_tem', 'y_pre', 'lnyield_obs', 'har']]

print(results_base_t)

       geoid  year       tmp       pre      har  irrigation1  lnyield_obs  \
0       1001  2009  24.14066  8.808839   1900.0     0.718851     7.609614   
1       1001  2010  25.67548  5.438990   1900.0     0.830087     7.035731   
2       1001  2011  24.54482  5.270162   1500.0     0.595086     7.466514   
3       1001  2013  24.03659  7.514143   2400.0     0.001241     7.926855   
4       1003  2008  24.72723  8.675528  19000.0     0.001980     7.934111   
...      ...   ...       ...       ...      ...          ...          ...   
17211  55141  2016  17.48690  7.304833  18200.0     0.000000     8.134343   
17212  55141  2017  16.59024  6.467205  19597.5     0.000000     7.978877   
17213  55141  2018  17.26666  8.373576  23000.0     0.000000     7.838077   
17214  55141  2019  16.23859  8.735889  10100.0     0.000000     7.912185   
17215  55141  2020  16.20263  6.832954  18500.0     0.000000     8.173033   

          lnfer   lnyield  group  ...  pre_pre_interaction  irrigation12  \
0      2.022537  7.900063      1  ...            77.595645      0.516747   
1      1.910005  7.661077      1  ...            29.582612      0.689044   
2      2.025641  7.785866      1  ...            27.774608      0.354128   
3      2.906756  8.084917      1  ...            56.462345      0.000002   
4      1.240602  7.807824      1  ...            75.264786      0.000004   
...         ...       ...    ...  ...                  ...           ...   
17211  2.860158  8.062482      2  ...            53.360585      0.000000   
17212  2.633171  7.905756      2  ...            41.824741      0.000000   
17213  2.639034  8.039181      2  ...            70.116775      0.000000   
17214  2.671287  7.918515      2  ...            76.315757      0.000000   
17215  2.625827  7.873821      2  ...            46.689260      0.000000   

       lnfer_irrigation1_tmp_interaction  \
0                              35.098172   
1                              40.707692   
2                              29.587092   
3                               0.086679   
4                               0.060734   
...                                  ...   
17211                           0.000000   
17212                           0.000000   
17213                           0.000000   
17214                           0.000000   
17215                           0.000000   

      lnfer_irrigation12_tmp_tmp_interaction  \
0                                 609.077444   
1                                 867.597822   
2                                 432.157531   
3                                   0.002585   
4                                   0.002973   
...                                      ...   
17211                               0.000000   
17212                               0.000000   
17213                               0.000000   
17214                               0.000000   
17215                               0.000000   

       lnfer_irrigation1_pre_interaction  \
0                              12.807195   
1                               8.623353   
2                               6.352818   
3                               0.027097   
4                               0.021308   
...                                  ...   
17211                           0.000000   
17212                           0.000000   
17213                           0.000000   
17214                           0.000000   
17215                           0.000000   

       lnfer_irrigation12_pre_pre_interaction  lnfer_tmp_interaction  \
0                                   81.098269              48.825378   
1                                   38.932997              49.040295   
2                                   19.923715              49.718994   
3                                    0.000253              69.868502   
4                                    0.000366              30.676651   
...                                       ...                    ...   
17211                                0.000000              50.015297   
17212                                0.000000              43.684939   
17213                                0.000000              45.567303   
17214                                0.000000              43.377934   
17215                                0.000000              42.545303   

       lnfer_tmp_tmp_interaction  lnfer_pre_interaction  \
0                    1178.676851              17.816203   
1                    1259.133118              10.388498   
2                    1220.343752              10.675456   
3                    1679.400541              21.841780   
4                     758.548605              10.762877   
...                          ...                    ...   
17211                 874.612496              20.892977   
17212                 724.743620              17.029257   
17213                 786.795125              22.098152   
17214                 704.396491              23.336067   
17215                 689.345808              17.942155   

       lnfer_pre_pre_interaction  
0                     156.940062  
1                      56.502937  
2                      56.261384  
3                     164.122260  
4                      93.373644  
...                          ...  
17211                 152.619705  
17212                 110.131694  
17213                 185.040553  
17214                 203.861289  
17215                 122.597920  

[17216 rows x 26 columns]
       geoid  year        tmp       pre        har  irrigation1  lnyield_obs  \
26081   1001  2015  25.099590  6.254395    779.625     0.377088     7.509678   
17216   1001  2015  26.784940  6.208599    779.625     0.377088     7.509678   
17218   1001  2015  25.010676  6.605178    779.625     0.377088     7.509678   
17222   1003  2015  26.261318  9.274925   6358.433     0.738441     7.927428   
17221   1003  2015  27.808318  8.926058   6358.433     0.738441     7.927428   
...      ...   ...        ...       ...        ...          ...          ...   
26072  55139  2015  20.012748  5.680057  14541.420     0.000000     8.028848   
26071  55139  2015  17.287318  5.693542  14541.420     0.000000     8.028848   
26077  55141  2015  17.398562  5.719637   6700.507     0.000000     7.921128   
26076  55141  2015  19.304302  6.116161   6700.507     0.000000     7.921128   
43802  55141  2015  16.504032  6.108389   6700.507     0.000000     7.921128   

          lnfer  group  sample  scenario   pre2015    tmp2015  
26081  2.457656      1       1       1.5  6.605178  25.010676  
17216  2.457656      1       1       3.0  6.605178  25.010676  
17218  2.457656      1       1    2015.0  6.605178  25.010676  
17222  1.292553      1       1       1.5  9.548671  26.127096  
17221  1.292553      1       1       3.0  9.548671  26.127096  
...         ...    ...     ...       ...       ...        ...  
26072  1.255133      2       1       3.0  5.693542  17.287318  
26071  1.255133      2       1    2015.0  5.693542  17.287318  
26077  2.556201      2       1       1.5  6.108389  16.504032  
26076  2.556201      2       1       3.0  6.108389  16.504032  
43802  2.556201      2       1    2015.0  6.108389  16.504032  

[5319 rows x 13 columns]
-0.013398123665888008
0.9408977132139006
3816
3816
3816
-0.00627099550076526
4.94855569752093
1503
1503
1503
      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
3      1003       1.5       7.590898  7.581375  7.609421     7.927428   
4      1003       3.0       7.351790  7.333040  7.618647     7.927428   
6      1005       1.5       7.807509  7.813926  7.814709     8.161582   
...     ...       ...            ...       ...       ...          ...   
5311  55137       3.0       8.006494  8.008905  7.852551     7.941248   
5313  55139       1.5       7.878360  7.896531  7.831343     8.028848   
5314  55139       3.0       7.962289  7.963008  7.848795     8.028848   
5316  55141       1.5       7.946194  7.971672  7.845699     7.921128   
5317  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  
0       779.625  
1       779.625  
3      6358.433  
4      6358.433  
6       526.500  
...         ...  
5311   5335.875  
5313  14541.420  
5314  14541.420  
5316   6700.507  
5317   6700.507  

[3546 rows x 7 columns]

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1921580644.py:173: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  results_base_t.rename(columns={'y_temp_base': 'y_temp_base_t'}, inplace=True)

# Load data from an Excel file
excel_file = "scenario_management.xlsx"
sheet_name = "soybean"
data = pd.read_excel(excel_file, sheet_name=sheet_name)

# Create interaction terms
data['tmp_tmp_interaction'] = data['tmp'] ** 2
data['pre_pre_interaction'] = data['pre'] ** 2
data['irrigation12'] = data['irrigation1'] ** 2
data['lnfer_irrigation1_tmp_interaction'] = data['lnfer'] * data['irrigation1'] * data['tmp']
data['lnfer_irrigation12_tmp_tmp_interaction'] = data['lnfer'] * data['irrigation12'] * data['tmp'] ** 2
data['lnfer_irrigation1_pre_interaction'] = data['lnfer'] * data['irrigation1'] * data['pre']
data['lnfer_irrigation12_pre_pre_interaction'] = data['lnfer'] * data['irrigation12'] * data['pre'] ** 2

data['lnfer_tmp_interaction'] = data['lnfer'] * data['tmp']
data['lnfer_tmp_tmp_interaction'] = data['lnfer'] * data['tmp'] ** 2
data['lnfer_pre_interaction'] = data['lnfer'] * data['pre']
data['lnfer_pre_pre_interaction'] = data['lnfer'] * data['pre'] ** 2

# Define predictor variables and response variable
predictor_variables1 = ['lnfer', 'tmp', 'irrigation1', 'irrigation12', 'tmp_tmp_interaction',
                        'pre', 'pre_pre_interaction', 'lnfer_irrigation1_pre_interaction', 'lnfer_irrigation12_pre_pre_interaction',
                        'lnfer_irrigation1_tmp_interaction', 'lnfer_irrigation12_tmp_tmp_interaction']

predictor_variables2 = ['lnfer', 'tmp', 'tmp_tmp_interaction', 'irrigation1', 'irrigation12', 'pre', 'pre_pre_interaction',
                        'lnfer_tmp_interaction', 'lnfer_tmp_tmp_interaction',
                        'lnfer_pre_interaction', 'lnfer_pre_pre_interaction']

response_variable = 'lnyield'

# Function to perform regression analysis
def perform_regression(group_data, predictor_variables):
    X = sm.add_constant(group_data[predictor_variables])
    y = group_data[response_variable]
    return sm.OLS(y, X).fit()

# Create a list to store data rows
data_rows = []

sample_data = data[data['sample'] == 0].copy()
project_data = data1

# Iterate through 'results' and extract values
for group in sample_data['group'].unique():
    group_data = sample_data[sample_data['group'] == group].copy()
    predictor_variables = predictor_variables1 if group == 1 else predictor_variables2
    reg = perform_regression(group_data, predictor_variables)

    for group_proj in project_data['group'].unique():
        if group == group_proj:
            # Calculate y_temp_base based on regression parameters
            tmp_values = project_data[project_data['group'] == group_proj]['tmp'].values
            pre_values = project_data[project_data['group'] == group_proj]['pre'].values
            lnfer_values = project_data[project_data['group'] == group_proj]['lnfer'].values
            irrigation1_values = project_data[project_data['group'] == group_proj]['irrigation1'].values
            
            if group_proj == 1:
                print(reg.params['tmp_tmp_interaction'])
                print(reg.params['const'])
                y_temp_base = (
                    reg.params['tmp_tmp_interaction'] * tmp_values**2 +
                    reg.params['tmp'] * tmp_values +
                    reg.params['pre_pre_interaction'] * pre_values**2 +
                    reg.params['pre'] * pre_values +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['irrigation1'] * irrigation1_values +
                    reg.params['irrigation12'] * irrigation1_values**2 +
                    reg.params['lnfer_irrigation1_tmp_interaction'] * irrigation1_values * tmp_values * lnfer_values +
                    reg.params['lnfer_irrigation12_tmp_tmp_interaction'] * irrigation1_values**2 * lnfer_values * tmp_values**2 +
                    reg.params['lnfer_irrigation1_pre_interaction'] * irrigation1_values * pre_values * lnfer_values +
                    reg.params['lnfer_irrigation12_pre_pre_interaction'] * irrigation1_values**2 * lnfer_values * pre_values**2 +
                    reg.params['const']
                )
                print(len(y_temp_base))
            else:
                print(reg.params['tmp_tmp_interaction'])
                y_temp_base = (
                    reg.params['tmp_tmp_interaction'] * tmp_values**2 +
                    reg.params['tmp'] * tmp_values +
                    reg.params['pre_pre_interaction'] * pre_values**2 +
                    reg.params['pre'] * pre_values +
                    reg.params['lnfer'] * lnfer_values +
                    reg.params['lnfer_tmp_interaction'] * tmp_values * lnfer_values +
                    reg.params['lnfer_tmp_tmp_interaction'] * lnfer_values * tmp_values**2 +
                    reg.params['lnfer_pre_interaction'] * pre_values * lnfer_values +
                    reg.params['lnfer_pre_pre_interaction'] * lnfer_values * pre_values**2 +
                    reg.params['const']
                )
                print(len(y_temp_base))
            
            # Extract 'geoid', 'scenario', and 'model' from the original data
            geoid_values = project_data[project_data['group'] == group_proj]['geoid'].values
            scenario_values = project_data[project_data['group'] == group_proj]['scenario'].values
            
            for tmp, pre, lnfer, irrigation1,y_temp_base_val, geoid, scenario in zip(
                tmp_values, pre_values, lnfer_values, irrigation1_values,
                 y_temp_base, geoid_values, scenario_values
            ):
                data_rows.append([tmp, pre, lnfer, irrigation1, y_temp_base_val, geoid, scenario])


# Create a DataFrame from the data
results_base = pd.DataFrame(data_rows, columns=['tmp','pre','lnfer','irrigation1',
                                               'y_temp_base', 'geoid', 'scenario'])

results_base = results_base[results_base['scenario'] == 2015]

# Rename the 'y_temp_base' column to 'y_temp_base2015'
results_base.rename(columns={'y_temp_base': 'y_temp_base2015'}, inplace=True)

# Keep only the desired columns in results_2015
results_base2015 = results_base[['geoid',  'y_temp_base2015']]

print(results_base2015)

-0.013398123665888008
0.9408977132139006
3816
-0.00627099550076526
1503
      geoid  y_temp_base2015
2      1001         7.871762
5      1003         7.599898
8      1005         7.821125
11     1009         8.285447
14     1011         8.201724
...     ...              ...
5306  55133         7.901357
5309  55135         7.816405
5312  55137         7.854962
5315  55139         7.849514
5318  55141         7.871178

[1773 rows x 2 columns]

data = results_base_t.merge(results_base2015, left_on=['geoid'], right_on=['geoid'], how='left')
data = data.sort_values(by=['geoid', 'scenario'])
print(data)

data['yield_warming']=np.exp(data['y_temp_base_t']-data['y_temp_base2015']+data['lnyield_obs'])
data['gap_yield_warming']=(data['yield_warming']-np.exp(data['lnyield_obs']))/1000 #tonne/ha
data['gap_production_warming']=data['gap_yield_warming']*data['har']/10000  #10000 tonne
print(data)

data['yield_tem']=np.exp(data['y_tem']-data['y_temp_base2015']+data['lnyield_obs'])
data['gap_yield_tem']=(data['yield_tem']-np.exp(data['lnyield_obs']))/1000 #tonne/ha
data['gap_production_tem']=data['gap_yield_tem']*data['har']/10000  #10000 tonne

data['yield_pre']=np.exp(data['y_pre']-data['y_temp_base2015']+data['lnyield_obs'])
data['gap_yield_pre']=(data['yield_pre']-np.exp(data['lnyield_obs']))/1000 #tonne/ha
data['gap_production_pre']=data['gap_yield_pre']*data['har']/10000  #10000 tonne

data['yield_obs']=(np.exp(data['lnyield_obs']))/1000 #tonne/ha
data['production_obs']=data['yield_obs']*data['har']/10000  #10000 tonne
print(data)
alldata=data

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
...     ...       ...            ...       ...       ...          ...   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  
0       779.625         7.871762  
1       779.625         7.871762  
2      6358.433         7.599898  
3      6358.433         7.599898  
4       526.500         7.821125  
...         ...              ...  
3541   5335.875         7.854962  
3542  14541.420         7.849514  
3543  14541.420         7.849514  
3544   6700.507         7.871178  
3545   6700.507         7.871178  

[3546 rows x 8 columns]
      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
...     ...       ...            ...       ...       ...          ...   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
1       779.625         7.871762    1423.263948          -0.402362   
2      6358.433         7.599898    2747.448365          -0.024839   
3      6358.433         7.599898    2163.149263          -0.609138   
4       526.500         7.821125    3456.341787          -0.047383   
...         ...              ...            ...                ...   
3541   5335.875         7.854962    3270.766621           0.459900   
3542  14541.420         7.849514    3157.999025           0.089794   
3543  14541.420         7.849514    3434.487053           0.366282   
3544   6700.507         7.871178    2969.485302           0.214609   
3545   6700.507         7.871178    3563.285992           0.808409   

      gap_production_warming  
0                  -0.003556  
1                  -0.031369  
2                  -0.015794  
3                  -0.387316  
4                  -0.002495  
...                      ...  
3541                0.245397  
3542                0.130573  
3543                0.532626  
3544                0.143799  
3545                0.541675  

[3546 rows x 11 columns]
      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
...     ...       ...            ...       ...       ...          ...   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
1       779.625         7.871762    1423.263948          -0.402362   
2      6358.433         7.599898    2747.448365          -0.024839   
3      6358.433         7.599898    2163.149263          -0.609138   
4       526.500         7.821125    3456.341787          -0.047383   
...         ...              ...            ...                ...   
3541   5335.875         7.854962    3270.766621           0.459900   
3542  14541.420         7.849514    3157.999025           0.089794   
3543  14541.420         7.849514    3434.487053           0.366282   
3544   6700.507         7.871178    2969.485302           0.214609   
3545   6700.507         7.871178    3563.285992           0.808409   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
0                  -0.003556  1808.246423      -0.017379           -0.001355   
1                  -0.031369  1449.169067      -0.376457           -0.029349   
2                  -0.015794  2721.408325      -0.050879           -0.032351   
3                  -0.387316  2122.968625      -0.649319           -0.412865   
4                  -0.002495  3478.590658      -0.025134           -0.001323   
...                      ...          ...            ...                 ...   
3541                0.245397  3278.662976       0.467797            0.249610   
3542                0.130573  3215.907385       0.147702            0.214780   
3543                0.532626  3436.956451       0.368751            0.536217   
3544                0.143799  3046.116440       0.291240            0.195145   
3545                0.541675  3561.607831       0.806731            0.540551   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
0     1797.117618      -0.028508           -0.002223   1.825626   
1     1792.990999      -0.032635           -0.002544   1.825626   
2     2798.814184       0.026527            0.016867   2.772287   
3     2824.757355       0.052470            0.033363   2.772287   
4     3481.315485      -0.022410           -0.001180   3.503725   
...           ...            ...                 ...        ...   
3541  2804.096571      -0.006770           -0.003612   2.810866   
3542  3012.956357      -0.055249           -0.080339   3.068205   
3543  3066.000609      -0.002204           -0.003206   3.068205   
3544  2685.572371      -0.069304           -0.046437   2.754877   
3545  2756.174840       0.001298            0.000870   2.754877   

      production_obs  
0           0.142330  
1           0.142330  
2           1.762740  
3           1.762740  
4           0.184471  
...              ...  
3541        1.499843  
3542        4.461606  
3543        4.461606  
3544        1.845907  
3545        1.845907  

[3546 rows x 19 columns]

# warming scenario
# Filter rows where 'scenario' is equal to 1.5
data = alldata[alldata['scenario'] == 1.5]
print(data)

print(np.min(data['gap_yield_tem']),np.max(data['gap_yield_tem']),np.mean(data['gap_yield_tem']))
print(np.percentile(data['gap_yield_tem'], 5),np.percentile(data['gap_yield_tem'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_yield_tem'] > 0.82, 'gap_yield_tem'] = 0.82
data.loc[data['gap_yield_tem'] < -0.82, 'gap_yield_tem'] = -0.82

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.RdYlBu(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-0.85, vmax=0.85)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-0.8, -0.4, 0, 0.4, 0.8]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Yield change (tonnes/ha/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Yield change (tonnes/ha/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
filtered_data.plot(column='gap_yield_tem', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(c) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_yield_15_tem.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
6      1009       1.5       8.252038  8.263953  8.273532     7.840701   
8      1011       1.5       8.183765  8.196448  8.189040     7.984614   
...     ...       ...            ...       ...       ...          ...   
3536  55133       1.5       7.939076  7.945049  7.895383     7.974899   
3538  55135       1.5       7.859488  7.881852  7.794041     7.957613   
3540  55137       1.5       7.893858  7.916888  7.831933     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
2      6358.433         7.599898    2747.448365          -0.024839   
4       526.500         7.821125    3456.341787          -0.047383   
6      1107.736         8.285447    2458.463020          -0.083523   
8       870.750         8.201724    2883.197758          -0.052246   
...         ...              ...            ...                ...   
3536   7837.872         7.901357    3018.810254           0.111746   
3538   9909.051         7.816405    2983.031280           0.125787   
3540   5335.875         7.854962    2922.351719           0.111485   
3542  14541.420         7.849514    3157.999025           0.089794   
3544   6700.507         7.871178    2969.485302           0.214609   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
0                  -0.003556  1808.246423      -0.017379           -0.001355   
2                  -0.015794  2721.408325      -0.050879           -0.032351   
4                  -0.002495  3478.590658      -0.025134           -0.001323   
6                  -0.009252  2487.931418      -0.054055           -0.005988   
8                  -0.004549  2919.998466      -0.015446           -0.001345   
...                      ...          ...            ...                 ...   
3536                0.087585  3036.896480       0.129832            0.101761   
3538                0.124643  3050.496205       0.193252            0.191494   
3540                0.059487  2990.432916       0.179567            0.095815   
3542                0.130573  3215.907385       0.147702            0.214780   
3544                0.143799  3046.116440       0.291240            0.195145   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
0     1797.117618      -0.028508           -0.002223   1.825626   
2     2798.814184       0.026527            0.016867   2.772287   
4     3481.315485      -0.022410           -0.001180   3.503725   
6     2511.877490      -0.030109           -0.003335   2.541986   
8     2898.448637      -0.036995           -0.003221   2.935444   
...           ...            ...                 ...        ...   
3536  2889.751233      -0.017313           -0.013570   2.907064   
3538  2794.053601      -0.063191           -0.062616   2.857245   
3540  2746.873159      -0.063993           -0.034146   2.810866   
3542  3012.956357      -0.055249           -0.080339   3.068205   
3544  2685.572371      -0.069304           -0.046437   2.754877   

      production_obs  
0           0.142330  
2           1.762740  
4           0.184471  
6           0.281585  
8           0.255604  
...              ...  
3536        2.278520  
3538        2.831258  
3540        1.499843  
3542        4.461606  
3544        1.845907  

[1773 rows x 19 columns]
-0.2705061645449869 0.3930858603063348 0.019464671523886125
-0.12008486829578888 0.2107192864307359

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/2378038445.py:13: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/2378038445.py:24: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/2378038445.py:27: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

# warming scenario
# Filter rows where 'scenario' is equal to 3.0
data = alldata[alldata['scenario'] == 3]

print(np.min(data['gap_yield_tem']),np.max(data['gap_yield_tem']),np.mean(data['gap_yield_tem']))
print(np.percentile(data['gap_yield_tem'], 5),np.percentile(data['gap_yield_tem'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_yield_tem'] > 0.82, 'gap_yield_tem'] = 0.82
data.loc[data['gap_yield_tem'] < -0.82, 'gap_yield_tem'] = -0.82

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.RdYlBu(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-0.85, vmax=0.85)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-0.8, -0.4, 0, 0.4, 0.8]
# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Yield change (tonnes/ha/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Yield change (tonnes/ha/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
filtered_data.plot(column='gap_yield_tem', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(g) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_yield_3_tem.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

-1.0246078643134875 1.424309730498316 -0.09531852705028614
-0.5732467023482948 0.6114727730280216

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1039305278.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1039305278.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1039305278.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

# warming scenario
# Filter rows where 'scenario' is equal to 1.5
data = alldata[alldata['scenario'] == 1.5]

print(np.min(data['gap_yield_pre']),np.max(data['gap_yield_pre']),np.mean(data['gap_yield_pre']))
print(np.percentile(data['gap_yield_pre'], 5),np.percentile(data['gap_yield_pre'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_yield_pre'] > 0.82, 'gap_yield_pre'] = 0.82
data.loc[data['gap_yield_pre'] < -0.82, 'gap_yield_pre'] = -0.82

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.RdYlBu(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-0.85, vmax=0.85)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-0.8, -0.4, 0, 0.4, 0.8]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Yield change (tonnes/ha/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Yield change (tonnes/ha/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
filtered_data.plot(column='gap_yield_pre', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(d) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_yield_15_pre.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

-0.22096733063125384 0.26111832478073305 -0.013044942645089862
-0.13051479600755855 0.1141240435924736

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1402444114.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1402444114.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1402444114.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

# warming scenario
# Filter rows where 'scenario' is equal to 3.0
data = alldata[alldata['scenario'] == 3]

print(np.min(data['gap_yield_pre']),np.max(data['gap_yield_pre']),np.mean(data['gap_yield_pre']))
print(np.percentile(data['gap_yield_pre'], 5),np.percentile(data['gap_yield_pre'], 95))

data.loc[data['gap_yield_pre'] > 0.82, 'gap_yield_pre'] = 0.82
data.loc[data['gap_yield_pre'] < -0.82, 'gap_yield_pre'] = -0.82

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.RdYlBu(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-0.85, vmax=0.85)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-0.8, -0.4, 0, 0.4, 0.8]
# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Yield change (tonnes/ha/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Yield change (tonnes/ha/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
filtered_data.plot(column='gap_yield_pre', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(h) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_yield_3_pre.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

-0.29892407330188825 0.23823591853341441 -0.014247440048533452
-0.13660592311274003 0.12175341699056687

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/31486314.py:11: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/31486314.py:22: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/31486314.py:25: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 1.5]

print(data)
print(np.min(data['gap_production_tem']),np.max(data['gap_production_tem']),np.mean(data['gap_production_tem']))
print(np.percentile(data['gap_production_tem'], 5),np.percentile(data['gap_production_tem'], 95))
print(np.sum(data['gap_production_tem']))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_tem'] > 2.1, 'gap_production_tem'] = 2.1
data.loc[data['gap_production_tem'] < -2.1, 'gap_production_tem'] = -2.1

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-2.2, vmax=2.2)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-2, -1, 0, 1, 2]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change\n($10^{4}$ tonnes/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change\n($10^{4}$ tonnes/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_tem', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(c) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_15_tem.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
6      1009       1.5       8.252038  8.263953  8.273532     7.840701   
8      1011       1.5       8.183765  8.196448  8.189040     7.984614   
...     ...       ...            ...       ...       ...          ...   
3536  55133       1.5       7.939076  7.945049  7.895383     7.974899   
3538  55135       1.5       7.859488  7.881852  7.794041     7.957613   
3540  55137       1.5       7.893858  7.916888  7.831933     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
2      6358.433         7.599898    2747.448365          -0.024839   
4       526.500         7.821125    3456.341787          -0.047383   
6      1107.736         8.285447    2458.463020          -0.083523   
8       870.750         8.201724    2883.197758          -0.052246   
...         ...              ...            ...                ...   
3536   7837.872         7.901357    3018.810254           0.111746   
3538   9909.051         7.816405    2983.031280           0.125787   
3540   5335.875         7.854962    2922.351719           0.111485   
3542  14541.420         7.849514    3157.999025           0.089794   
3544   6700.507         7.871178    2969.485302           0.214609   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
0                  -0.003556  1808.246423      -0.017379           -0.001355   
2                  -0.015794  2721.408325      -0.050879           -0.032351   
4                  -0.002495  3478.590658      -0.025134           -0.001323   
6                  -0.009252  2487.931418      -0.054055           -0.005988   
8                  -0.004549  2919.998466      -0.015446           -0.001345   
...                      ...          ...            ...                 ...   
3536                0.087585  3036.896480       0.129832            0.101761   
3538                0.124643  3050.496205       0.193252            0.191494   
3540                0.059487  2990.432916       0.179567            0.095815   
3542                0.130573  3215.907385       0.147702            0.214780   
3544                0.143799  3046.116440       0.291240            0.195145   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
0     1797.117618      -0.028508           -0.002223   1.825626   
2     2798.814184       0.026527            0.016867   2.772287   
4     3481.315485      -0.022410           -0.001180   3.503725   
6     2511.877490      -0.030109           -0.003335   2.541986   
8     2898.448637      -0.036995           -0.003221   2.935444   
...           ...            ...                 ...        ...   
3536  2889.751233      -0.017313           -0.013570   2.907064   
3538  2794.053601      -0.063191           -0.062616   2.857245   
3540  2746.873159      -0.063993           -0.034146   2.810866   
3542  3012.956357      -0.055249           -0.080339   3.068205   
3544  2685.572371      -0.069304           -0.046437   2.754877   

      production_obs  
0           0.142330  
2           1.762740  
4           0.184471  
6           0.281585  
8           0.255604  
...              ...  
3536        2.278520  
3538        2.831258  
3540        1.499843  
3542        4.461606  
3544        1.845907  

[1773 rows x 19 columns]
-1.7232100621426445 2.360531492628775 0.06837798599747648
-0.28157389572941444 0.738641855720369
121.2341691735258

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3711736066.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3711736066.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3711736066.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 3]

print(data)
print(np.min(data['gap_production_tem']),np.max(data['gap_production_tem']),np.mean(data['gap_production_tem']))
print(np.percentile(data['gap_production_tem'], 5),np.percentile(data['gap_production_tem'], 95))
print(np.sum(data['gap_production_tem']))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_tem'] > 2.1, 'gap_production_tem'] = 2.1
data.loc[data['gap_production_tem'] < -2.1, 'gap_production_tem'] = -2.1

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-2.2, vmax=2.2)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-2, -1, 0, 1, 2]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change\n($10^{4}$ tonnes/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change\n($10^{4}$ tonnes/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_tem', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(g) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_3_tem.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
5      1005       3.0       7.574791  7.580497  7.815420     8.161582   
7      1009       3.0       8.067909  8.085865  8.267491     7.840701   
9      1011       3.0       7.957275  7.967675  8.191323     7.984614   
...     ...       ...            ...       ...       ...          ...   
3537  55133       3.0       8.016329  8.009911  7.907775     7.974899   
3539  55135       3.0       7.978732  7.985286  7.809852     7.957613   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
1       779.625         7.871762    1423.263948          -0.402362   
3      6358.433         7.599898    2163.149263          -0.609138   
5       526.500         7.821125    2738.726589          -0.764999   
7      1107.736         8.285447    2045.019309          -0.496967   
9       870.750         8.201724    2298.851794          -0.636592   
...         ...              ...            ...                ...   
3537   7837.872         7.901357    3261.266121           0.354202   
3539   9909.051         7.816405    3360.818522           0.503574   
3541   5335.875         7.854962    3270.766621           0.459900   
3543  14541.420         7.849514    3434.487053           0.366282   
3545   6700.507         7.871178    3563.285992           0.808409   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
1                  -0.031369  1449.169067      -0.376457           -0.029349   
3                  -0.387316  2122.968625      -0.649319           -0.412865   
5                  -0.040277  2754.397273      -0.749328           -0.039452   
7                  -0.055051  2082.072183      -0.459914           -0.050946   
9                  -0.055431  2322.886600      -0.612557           -0.053338   
...                      ...          ...            ...                 ...   
3537                0.277619  3240.402785       0.333339            0.261266   
3539                0.498994  3382.916708       0.525672            0.520891   
3541                0.245397  3278.662976       0.467797            0.249610   
3543                0.532626  3436.956451       0.368751            0.536217   
3545                0.541675  3561.607831       0.806731            0.540551   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
1     1792.990999      -0.032635           -0.002544   1.825626   
3     2824.757355       0.052470            0.033363   2.772287   
5     3483.791254      -0.019934           -0.001050   3.503725   
7     2496.748569      -0.045238           -0.005011   2.541986   
9     2905.071094      -0.030373           -0.002645   2.935444   
...           ...            ...                 ...        ...   
3537  2925.781379       0.018717            0.014670   2.907064   
3539  2838.580219      -0.018664           -0.018495   2.857245   
3541  2804.096571      -0.006770           -0.003612   2.810866   
3543  3066.000609      -0.002204           -0.003206   3.068205   
3545  2756.174840       0.001298            0.000870   2.754877   

      production_obs  
1           0.142330  
3           1.762740  
5           0.184471  
7           0.281585  
9           0.255604  
...              ...  
3537        2.278520  
3539        2.831258  
3541        1.499843  
3543        4.461606  
3545        1.845907  

[1773 rows x 19 columns]
-7.1569150862130995 9.991046513367879 -0.034810631570730985
-1.503981910801795 1.7494962914269376
-61.71924977490603

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/4156552768.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/4156552768.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/4156552768.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 1.5]

print(data)
print(np.min(data['gap_production_pre']),np.max(data['gap_production_pre']),np.mean(data['gap_production_pre']))
print(np.percentile(data['gap_production_pre'], 5),np.percentile(data['gap_production_pre'], 95))
print(np.sum(data['gap_production_pre']))

p95 = np.percentile(data['gap_production_pre'], 95)
# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_pre'] > 2.1, 'gap_production_pre'] = 2.1
data.loc[data['gap_production_pre'] < -2.1, 'gap_production_pre'] = -2.1

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-2.2, vmax=2.2)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-2, -1, 0, 1, 2]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change\n($10^{4}$ tonnes/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change\n($10^{4}$ tonnes/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_pre', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(d) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_15_pre.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
6      1009       1.5       8.252038  8.263953  8.273532     7.840701   
8      1011       1.5       8.183765  8.196448  8.189040     7.984614   
...     ...       ...            ...       ...       ...          ...   
3536  55133       1.5       7.939076  7.945049  7.895383     7.974899   
3538  55135       1.5       7.859488  7.881852  7.794041     7.957613   
3540  55137       1.5       7.893858  7.916888  7.831933     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
2      6358.433         7.599898    2747.448365          -0.024839   
4       526.500         7.821125    3456.341787          -0.047383   
6      1107.736         8.285447    2458.463020          -0.083523   
8       870.750         8.201724    2883.197758          -0.052246   
...         ...              ...            ...                ...   
3536   7837.872         7.901357    3018.810254           0.111746   
3538   9909.051         7.816405    2983.031280           0.125787   
3540   5335.875         7.854962    2922.351719           0.111485   
3542  14541.420         7.849514    3157.999025           0.089794   
3544   6700.507         7.871178    2969.485302           0.214609   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
0                  -0.003556  1808.246423      -0.017379           -0.001355   
2                  -0.015794  2721.408325      -0.050879           -0.032351   
4                  -0.002495  3478.590658      -0.025134           -0.001323   
6                  -0.009252  2487.931418      -0.054055           -0.005988   
8                  -0.004549  2919.998466      -0.015446           -0.001345   
...                      ...          ...            ...                 ...   
3536                0.087585  3036.896480       0.129832            0.101761   
3538                0.124643  3050.496205       0.193252            0.191494   
3540                0.059487  2990.432916       0.179567            0.095815   
3542                0.130573  3215.907385       0.147702            0.214780   
3544                0.143799  3046.116440       0.291240            0.195145   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
0     1797.117618      -0.028508           -0.002223   1.825626   
2     2798.814184       0.026527            0.016867   2.772287   
4     3481.315485      -0.022410           -0.001180   3.503725   
6     2511.877490      -0.030109           -0.003335   2.541986   
8     2898.448637      -0.036995           -0.003221   2.935444   
...           ...            ...                 ...        ...   
3536  2889.751233      -0.017313           -0.013570   2.907064   
3538  2794.053601      -0.063191           -0.062616   2.857245   
3540  2746.873159      -0.063993           -0.034146   2.810866   
3542  3012.956357      -0.055249           -0.080339   3.068205   
3544  2685.572371      -0.069304           -0.046437   2.754877   

      production_obs  
0           0.142330  
2           1.762740  
4           0.184471  
6           0.281585  
8           0.255604  
...              ...  
3536        2.278520  
3538        2.831258  
3540        1.499843  
3542        4.461606  
3544        1.845907  

[1773 rows x 19 columns]
-1.9153194352568128 1.098382750376755 -0.043688019711214915
-0.45890603154174003 0.21905820475850707
-77.45885894798404

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3150142858.py:13: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3150142858.py:24: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3150142858.py:27: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 3]

print(data)
print(np.min(data['gap_production_pre']),np.max(data['gap_production_pre']),np.mean(data['gap_production_pre']))
print(np.percentile(data['gap_production_pre'], 5),np.percentile(data['gap_production_pre'], 95))
print(np.sum(data['gap_production_pre']))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_pre'] > 2.1, 'gap_production_pre'] = 2.1
data.loc[data['gap_production_pre'] < -2.1, 'gap_production_pre'] = -2.1

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-2.2, vmax=2.2)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-2, -1, 0, 1, 2]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change\n($10^{4}$ tonnes/yr)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change\n($10^{4}$ tonnes/yr)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_pre', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(h) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_3_pre.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
5      1005       3.0       7.574791  7.580497  7.815420     8.161582   
7      1009       3.0       8.067909  8.085865  8.267491     7.840701   
9      1011       3.0       7.957275  7.967675  8.191323     7.984614   
...     ...       ...            ...       ...       ...          ...   
3537  55133       3.0       8.016329  8.009911  7.907775     7.974899   
3539  55135       3.0       7.978732  7.985286  7.809852     7.957613   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
1       779.625         7.871762    1423.263948          -0.402362   
3      6358.433         7.599898    2163.149263          -0.609138   
5       526.500         7.821125    2738.726589          -0.764999   
7      1107.736         8.285447    2045.019309          -0.496967   
9       870.750         8.201724    2298.851794          -0.636592   
...         ...              ...            ...                ...   
3537   7837.872         7.901357    3261.266121           0.354202   
3539   9909.051         7.816405    3360.818522           0.503574   
3541   5335.875         7.854962    3270.766621           0.459900   
3543  14541.420         7.849514    3434.487053           0.366282   
3545   6700.507         7.871178    3563.285992           0.808409   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
1                  -0.031369  1449.169067      -0.376457           -0.029349   
3                  -0.387316  2122.968625      -0.649319           -0.412865   
5                  -0.040277  2754.397273      -0.749328           -0.039452   
7                  -0.055051  2082.072183      -0.459914           -0.050946   
9                  -0.055431  2322.886600      -0.612557           -0.053338   
...                      ...          ...            ...                 ...   
3537                0.277619  3240.402785       0.333339            0.261266   
3539                0.498994  3382.916708       0.525672            0.520891   
3541                0.245397  3278.662976       0.467797            0.249610   
3543                0.532626  3436.956451       0.368751            0.536217   
3545                0.541675  3561.607831       0.806731            0.540551   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
1     1792.990999      -0.032635           -0.002544   1.825626   
3     2824.757355       0.052470            0.033363   2.772287   
5     3483.791254      -0.019934           -0.001050   3.503725   
7     2496.748569      -0.045238           -0.005011   2.541986   
9     2905.071094      -0.030373           -0.002645   2.935444   
...           ...            ...                 ...        ...   
3537  2925.781379       0.018717            0.014670   2.907064   
3539  2838.580219      -0.018664           -0.018495   2.857245   
3541  2804.096571      -0.006770           -0.003612   2.810866   
3543  3066.000609      -0.002204           -0.003206   3.068205   
3545  2756.174840       0.001298            0.000870   2.754877   

      production_obs  
1           0.142330  
3           1.762740  
5           0.184471  
7           0.281585  
9           0.255604  
...              ...  
3537        2.278520  
3539        2.831258  
3541        1.499843  
3543        4.461606  
3545        1.845907  

[1773 rows x 19 columns]
-1.6554353565582052 1.0917603521739567 -0.038505560139215286
-0.386414719208073 0.22641246435388723
-68.2703581268287

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3329200732.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3329200732.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3329200732.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

Percentage

data = alldata[alldata['scenario'] == 1.5]

print(data)
data['gap_production_tem_p']=data['gap_production_tem']/data['production_obs']*100 # change in percentage
print(np.min(data['gap_production_tem_p']),np.max(data['gap_production_tem_p']),np.mean(data['gap_production_tem_p']))
print(np.percentile(data['gap_production_tem_p'], 5),np.percentile(data['gap_production_tem_p'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_tem_p'] > 20.3, 'gap_production_tem_p'] = 20.3
data.loc[data['gap_production_tem_p'] < -20.3, 'gap_production_tem_p'] = -20.3

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-20.5, vmax=20.5)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-20, -10, 0, 10, 20]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change (%)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change (%)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_tem_p', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(c) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_15_tem_p.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
6      1009       1.5       8.252038  8.263953  8.273532     7.840701   
8      1011       1.5       8.183765  8.196448  8.189040     7.984614   
...     ...       ...            ...       ...       ...          ...   
3536  55133       1.5       7.939076  7.945049  7.895383     7.974899   
3538  55135       1.5       7.859488  7.881852  7.794041     7.957613   
3540  55137       1.5       7.893858  7.916888  7.831933     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
2      6358.433         7.599898    2747.448365          -0.024839   
4       526.500         7.821125    3456.341787          -0.047383   
6      1107.736         8.285447    2458.463020          -0.083523   
8       870.750         8.201724    2883.197758          -0.052246   
...         ...              ...            ...                ...   
3536   7837.872         7.901357    3018.810254           0.111746   
3538   9909.051         7.816405    2983.031280           0.125787   
3540   5335.875         7.854962    2922.351719           0.111485   
3542  14541.420         7.849514    3157.999025           0.089794   
3544   6700.507         7.871178    2969.485302           0.214609   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
0                  -0.003556  1808.246423      -0.017379           -0.001355   
2                  -0.015794  2721.408325      -0.050879           -0.032351   
4                  -0.002495  3478.590658      -0.025134           -0.001323   
6                  -0.009252  2487.931418      -0.054055           -0.005988   
8                  -0.004549  2919.998466      -0.015446           -0.001345   
...                      ...          ...            ...                 ...   
3536                0.087585  3036.896480       0.129832            0.101761   
3538                0.124643  3050.496205       0.193252            0.191494   
3540                0.059487  2990.432916       0.179567            0.095815   
3542                0.130573  3215.907385       0.147702            0.214780   
3544                0.143799  3046.116440       0.291240            0.195145   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
0     1797.117618      -0.028508           -0.002223   1.825626   
2     2798.814184       0.026527            0.016867   2.772287   
4     3481.315485      -0.022410           -0.001180   3.503725   
6     2511.877490      -0.030109           -0.003335   2.541986   
8     2898.448637      -0.036995           -0.003221   2.935444   
...           ...            ...                 ...        ...   
3536  2889.751233      -0.017313           -0.013570   2.907064   
3538  2794.053601      -0.063191           -0.062616   2.857245   
3540  2746.873159      -0.063993           -0.034146   2.810866   
3542  3012.956357      -0.055249           -0.080339   3.068205   
3544  2685.572371      -0.069304           -0.046437   2.754877   

      production_obs  
0           0.142330  
2           1.762740  
4           0.184471  
6           0.281585  
8           0.255604  
...              ...  
3536        2.278520  
3538        2.831258  
3540        1.499843  
3542        4.461606  
3544        1.845907  

[1773 rows x 19 columns]
-9.505131950352416 15.42906908839307 0.5406694107684533
-5.430827420656857 7.057769469297233

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1975420342.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['gap_production_tem_p']=data['gap_production_tem']/data['production_obs']*100 # change in percentage
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1975420342.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1975420342.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1975420342.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 3]

print(data)
data['gap_production_tem_p']=data['gap_production_tem']/data['production_obs']*100 # change in percentage
print(np.min(data['gap_production_tem_p']),np.max(data['gap_production_tem_p']),np.mean(data['gap_production_tem_p']))
print(np.percentile(data['gap_production_tem_p'], 5),np.percentile(data['gap_production_tem_p'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_tem_p'] > 20.3, 'gap_production_tem_p'] = 20.3
data.loc[data['gap_production_tem_p'] < -20.3, 'gap_production_tem_p'] = -20.3

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-20.5, vmax=20.5)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-20, -10, 0, 10, 20]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change (%)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change (%)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_tem_p', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(g) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_3_tem_p.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
5      1005       3.0       7.574791  7.580497  7.815420     8.161582   
7      1009       3.0       8.067909  8.085865  8.267491     7.840701   
9      1011       3.0       7.957275  7.967675  8.191323     7.984614   
...     ...       ...            ...       ...       ...          ...   
3537  55133       3.0       8.016329  8.009911  7.907775     7.974899   
3539  55135       3.0       7.978732  7.985286  7.809852     7.957613   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
1       779.625         7.871762    1423.263948          -0.402362   
3      6358.433         7.599898    2163.149263          -0.609138   
5       526.500         7.821125    2738.726589          -0.764999   
7      1107.736         8.285447    2045.019309          -0.496967   
9       870.750         8.201724    2298.851794          -0.636592   
...         ...              ...            ...                ...   
3537   7837.872         7.901357    3261.266121           0.354202   
3539   9909.051         7.816405    3360.818522           0.503574   
3541   5335.875         7.854962    3270.766621           0.459900   
3543  14541.420         7.849514    3434.487053           0.366282   
3545   6700.507         7.871178    3563.285992           0.808409   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
1                  -0.031369  1449.169067      -0.376457           -0.029349   
3                  -0.387316  2122.968625      -0.649319           -0.412865   
5                  -0.040277  2754.397273      -0.749328           -0.039452   
7                  -0.055051  2082.072183      -0.459914           -0.050946   
9                  -0.055431  2322.886600      -0.612557           -0.053338   
...                      ...          ...            ...                 ...   
3537                0.277619  3240.402785       0.333339            0.261266   
3539                0.498994  3382.916708       0.525672            0.520891   
3541                0.245397  3278.662976       0.467797            0.249610   
3543                0.532626  3436.956451       0.368751            0.536217   
3545                0.541675  3561.607831       0.806731            0.540551   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
1     1792.990999      -0.032635           -0.002544   1.825626   
3     2824.757355       0.052470            0.033363   2.772287   
5     3483.791254      -0.019934           -0.001050   3.503725   
7     2496.748569      -0.045238           -0.005011   2.541986   
9     2905.071094      -0.030373           -0.002645   2.935444   
...           ...            ...                 ...        ...   
3537  2925.781379       0.018717            0.014670   2.907064   
3539  2838.580219      -0.018664           -0.018495   2.857245   
3541  2804.096571      -0.006770           -0.003612   2.810866   
3543  3066.000609      -0.002204           -0.003206   3.068205   
3545  2756.174840       0.001298            0.000870   2.754877   

      production_obs  
1           0.142330  
3           1.762740  
5           0.184471  
7           0.281585  
9           0.255604  
...              ...  
3537        2.278520  
3539        2.831258  
3541        1.499843  
3543        4.461606  
3545        1.845907  

[1773 rows x 19 columns]
-30.556039700742744 51.456837077503124 -4.046315695932895
-22.196113893756138 24.00336699857401

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3641423086.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['gap_production_tem_p']=data['gap_production_tem']/data['production_obs']*100 # change in percentage
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3641423086.py:12: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3641423086.py:23: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3641423086.py:26: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 1.5]

print(data)
data['gap_production_pre_p']=data['gap_production_pre']/data['production_obs']*100 # change in percentage

print(np.min(data['gap_production_pre_p']),np.max(data['gap_production_pre_p']),np.mean(data['gap_production_pre_p']))
print(np.percentile(data['gap_production_pre_p'], 5),np.percentile(data['gap_production_pre_p'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_pre_p'] > 20.3, 'gap_production_pre_p'] = 20.3
data.loc[data['gap_production_pre_p'] < -20.3, 'gap_production_pre_p'] = -20.3

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-20.5, vmax=20.5)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-20, -10, 0, 10, 20]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change (%)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change (%)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_pre_p', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(d) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_15_pre_p.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
0      1001       1.5       7.846458  7.862197  7.856023     7.509678   
2      1003       1.5       7.590898  7.581375  7.609421     7.927428   
4      1005       1.5       7.807509  7.813926  7.814709     8.161582   
6      1009       1.5       8.252038  8.263953  8.273532     7.840701   
8      1011       1.5       8.183765  8.196448  8.189040     7.984614   
...     ...       ...            ...       ...       ...          ...   
3536  55133       1.5       7.939076  7.945049  7.895383     7.974899   
3538  55135       1.5       7.859488  7.881852  7.794041     7.957613   
3540  55137       1.5       7.893858  7.916888  7.831933     7.941248   
3542  55139       1.5       7.878360  7.896531  7.831343     8.028848   
3544  55141       1.5       7.946194  7.971672  7.845699     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
0       779.625         7.871762    1780.009827          -0.045616   
2      6358.433         7.599898    2747.448365          -0.024839   
4       526.500         7.821125    3456.341787          -0.047383   
6      1107.736         8.285447    2458.463020          -0.083523   
8       870.750         8.201724    2883.197758          -0.052246   
...         ...              ...            ...                ...   
3536   7837.872         7.901357    3018.810254           0.111746   
3538   9909.051         7.816405    2983.031280           0.125787   
3540   5335.875         7.854962    2922.351719           0.111485   
3542  14541.420         7.849514    3157.999025           0.089794   
3544   6700.507         7.871178    2969.485302           0.214609   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
0                  -0.003556  1808.246423      -0.017379           -0.001355   
2                  -0.015794  2721.408325      -0.050879           -0.032351   
4                  -0.002495  3478.590658      -0.025134           -0.001323   
6                  -0.009252  2487.931418      -0.054055           -0.005988   
8                  -0.004549  2919.998466      -0.015446           -0.001345   
...                      ...          ...            ...                 ...   
3536                0.087585  3036.896480       0.129832            0.101761   
3538                0.124643  3050.496205       0.193252            0.191494   
3540                0.059487  2990.432916       0.179567            0.095815   
3542                0.130573  3215.907385       0.147702            0.214780   
3544                0.143799  3046.116440       0.291240            0.195145   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
0     1797.117618      -0.028508           -0.002223   1.825626   
2     2798.814184       0.026527            0.016867   2.772287   
4     3481.315485      -0.022410           -0.001180   3.503725   
6     2511.877490      -0.030109           -0.003335   2.541986   
8     2898.448637      -0.036995           -0.003221   2.935444   
...           ...            ...                 ...        ...   
3536  2889.751233      -0.017313           -0.013570   2.907064   
3538  2794.053601      -0.063191           -0.062616   2.857245   
3540  2746.873159      -0.063993           -0.034146   2.810866   
3542  3012.956357      -0.055249           -0.080339   3.068205   
3544  2685.572371      -0.069304           -0.046437   2.754877   

      production_obs  
0           0.142330  
2           1.762740  
4           0.184471  
6           0.281585  
8           0.255604  
...              ...  
3536        2.278520  
3538        2.831258  
3540        1.499843  
3542        4.461606  
3544        1.845907  

[1773 rows x 19 columns]
-9.942069219302383 7.9464744349771745 -0.5290393162376932
-4.943241780596144 4.090930799241238

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/705259917.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['gap_production_pre_p']=data['gap_production_pre']/data['production_obs']*100 # change in percentage
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/705259917.py:13: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/705259917.py:24: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/705259917.py:27: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 3.0]

print(data)
data['gap_production_pre_p']=data['gap_production_pre']/data['production_obs']*100 # change in percentage

print(np.min(data['gap_production_pre_p']),np.max(data['gap_production_pre_p']),np.mean(data['gap_production_pre_p']))
print(np.percentile(data['gap_production_pre_p'], 5),np.percentile(data['gap_production_pre_p'], 95))

# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_pre_p'] > 20.3, 'gap_production_pre_p'] = 20.3
data.loc[data['gap_production_pre_p'] < -20.3, 'gap_production_pre_p'] = -20.3

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-20.5, vmax=20.5)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-20, -10, 0, 10, 20]

# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change (%)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)

# Set the font size of the colorbar label
# Set the colorbar label text, font size, and padding
cbar.ax.set_ylabel('Production change (%)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)  # Adjust tick label fontsize
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')

filtered_data.plot(column='gap_production_pre_p', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
plt.text(0.25, 0.2, '(h) Soybean', transform=plt.gcf().transFigure,
         fontsize=18, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')

# Save the entire figure as a JPG with DPI 300
plt.savefig('soybean_all_production_3_pre_p.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

      geoid  scenario  y_temp_base_t     y_tem     y_pre  lnyield_obs  \
1      1001       3.0       7.622792  7.640830  7.853724     7.509678   
3      1003       3.0       7.351790  7.333040  7.618647     7.927428   
5      1005       3.0       7.574791  7.580497  7.815420     8.161582   
7      1009       3.0       8.067909  8.085865  8.267491     7.840701   
9      1011       3.0       7.957275  7.967675  8.191323     7.984614   
...     ...       ...            ...       ...       ...          ...   
3537  55133       3.0       8.016329  8.009911  7.907775     7.974899   
3539  55135       3.0       7.978732  7.985286  7.809852     7.957613   
3541  55137       3.0       8.006494  8.008905  7.852551     7.941248   
3543  55139       3.0       7.962289  7.963008  7.848795     8.028848   
3545  55141       3.0       8.128488  8.128017  7.871649     7.921128   

            har  y_temp_base2015  yield_warming  gap_yield_warming  \
1       779.625         7.871762    1423.263948          -0.402362   
3      6358.433         7.599898    2163.149263          -0.609138   
5       526.500         7.821125    2738.726589          -0.764999   
7      1107.736         8.285447    2045.019309          -0.496967   
9       870.750         8.201724    2298.851794          -0.636592   
...         ...              ...            ...                ...   
3537   7837.872         7.901357    3261.266121           0.354202   
3539   9909.051         7.816405    3360.818522           0.503574   
3541   5335.875         7.854962    3270.766621           0.459900   
3543  14541.420         7.849514    3434.487053           0.366282   
3545   6700.507         7.871178    3563.285992           0.808409   

      gap_production_warming    yield_tem  gap_yield_tem  gap_production_tem  \
1                  -0.031369  1449.169067      -0.376457           -0.029349   
3                  -0.387316  2122.968625      -0.649319           -0.412865   
5                  -0.040277  2754.397273      -0.749328           -0.039452   
7                  -0.055051  2082.072183      -0.459914           -0.050946   
9                  -0.055431  2322.886600      -0.612557           -0.053338   
...                      ...          ...            ...                 ...   
3537                0.277619  3240.402785       0.333339            0.261266   
3539                0.498994  3382.916708       0.525672            0.520891   
3541                0.245397  3278.662976       0.467797            0.249610   
3543                0.532626  3436.956451       0.368751            0.536217   
3545                0.541675  3561.607831       0.806731            0.540551   

        yield_pre  gap_yield_pre  gap_production_pre  yield_obs  \
1     1792.990999      -0.032635           -0.002544   1.825626   
3     2824.757355       0.052470            0.033363   2.772287   
5     3483.791254      -0.019934           -0.001050   3.503725   
7     2496.748569      -0.045238           -0.005011   2.541986   
9     2905.071094      -0.030373           -0.002645   2.935444   
...           ...            ...                 ...        ...   
3537  2925.781379       0.018717            0.014670   2.907064   
3539  2838.580219      -0.018664           -0.018495   2.857245   
3541  2804.096571      -0.006770           -0.003612   2.810866   
3543  3066.000609      -0.002204           -0.003206   3.068205   
3545  2756.174840       0.001298            0.000870   2.754877   

      production_obs  
1           0.142330  
3           1.762740  
5           0.184471  
7           0.281585  
9           0.255604  
...              ...  
3537        2.278520  
3539        2.831258  
3541        1.499843  
3543        4.461606  
3545        1.845907  

[1773 rows x 19 columns]
-11.74440851900627 7.533661361393064 -0.6096099246111667
-5.852965407030569 4.41937139595061

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1548580543.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['gap_production_pre_p']=data['gap_production_pre']/data['production_obs']*100 # change in percentage
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1548580543.py:13: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1548580543.py:24: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/1548580543.py:27: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)

data = alldata[alldata['scenario'] == 3]

# Replace values in 'gap_yield' greater than p90 with p90
data['gap_production_pre_p']=data['gap_production_pre']/data['production_obs']*100 # change in percentage
data.loc[data['gap_production_pre_p'] > 20.3, 'gap_production_pre_p'] = 20.3
data.loc[data['gap_production_pre_p'] < -20.3, 'gap_production_pre_p'] = -20.3

data['geoid'] = data['geoid'].astype(float)


# Load US county boundaries shapefile using geopandas
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_county_20m/surface_shapefile.shp'
gdf_county = gpd.read_file(shapefile_path)

# Ensure 'GEOID' column in the shapefile has the 'object' data type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)

# Rename the 'geoid' column to 'GEOID' in year_2020_data
data.rename(columns={'geoid': 'GEOID'}, inplace=True)

# Convert the 'GEOID' column in the 'data' DataFrame to 'object' data type
data['GEOID'] = data['GEOID'].astype(float)
# Now, merge 'data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(data, on='GEOID', how='left')

# Filter the data to include only geometries within the specified longitude and latitude range
filtered_data = merged_data.cx[-130:-60, 25:50]

cmap = mcolors.ListedColormap(['white'] + list(plt.cm.Spectral(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-20.5, vmax=20.5)

# Plot both county boundaries and 'y_temp' values on the same map within the specified range
fig, ax = plt.subplots(figsize=(10, 4))
shapefile_path = '/Users/chenchenren/postdoc/paper/2N and water-US/Data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp'
gdf_state = gpd.read_file(shapefile_path)
gdf_state.rename(columns={'GEOID': 'GEOID_state'}, inplace=True)
gdf_state = gdf_state.cx[-130:-60, 25:50]
gdf_state.boundary.plot(ax=ax, linewidth=0.3, color='grey')  # Plot county boundaries in grey

sc = plt.scatter([], [], c=[], cmap=cmap, norm=norm)  # Create an empty ScalarMappable for the colorbar

# Define custom colorbar ticks
custom_ticks = [-20, -10, 0, 10, 20]

# Adjust the height of the colorbar by setting the shrink parameter
# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Production change (%)', orientation='horizontal', pad=0.02, shrink=0.8, ticks=custom_ticks)
# Set the colorbar label text, font size, and padding
cbar.ax.set_xlabel('Production change (%)', fontsize=20, fontfamily='Arial')
cbar.ax.tick_params(labelsize=20)  # Adjust tick label fontsize
cbar.set_ticks(custom_ticks)


filtered_data.plot(column='gap_production_pre', cmap=cmap, ax=ax, norm=norm)

# Add labels, titles, or other plot configurations as needed
ax.set_xticks([])  # Remove x-axis ticks
ax.set_yticks([])  # Remove y-axis ticks
ax.set_frame_on(False)  # Turn off the frame
# Save the entire figure as a JPG with DPI 300
plt.savefig('production_p_soybean_legend.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)

# Show the plot
plt.show()

/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3522885677.py:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['gap_production_pre_p']=data['gap_production_pre']/data['production_obs']*100 # change in percentage
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3522885677.py:8: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['geoid'] = data['geoid'].astype(float)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3522885677.py:19: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data.rename(columns={'geoid': 'GEOID'}, inplace=True)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_7700/3522885677.py:22: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  data['GEOID'] = data['GEOID'].astype(float)