In [5]:
# Load your data
excel_file = "sensitivity.xlsx"
sheet_name = "maize"
data = pd.read_excel(excel_file, sheet_name=sheet_name)
# Create interaction terms
data['tmp_tmp_interaction'] = data['tmp'] * data['tmp']
data['irrigation12'] = data['irrigation1'] * data['irrigation1']
data['lnfer_irrigation1_tmp_interaction'] = data['lnfer'] * data['irrigation1'] * data['tmp']
data['lnfer_irrigation12_tmp_tmp_interaction'] = data['lnfer']* data['irrigation1'] * data['irrigation1'] * data['tmp'] * data['tmp']
data['lnfer_tmp_interaction'] = data['lnfer'] * data['tmp']
data['lnfer_tmp_tmp_interaction'] = data['lnfer'] * data['tmp'] * data['tmp']
# Define predictor variables and response variable
predictor_variables1 = ['lnfer', 'tmp', 'irrigation1', 'irrigation12', 'tmp_tmp_interaction',
'lnfer_irrigation1_tmp_interaction', 'lnfer_irrigation12_tmp_tmp_interaction']
predictor_variables2 = ['lnfer', 'tmp', 'tmp_tmp_interaction', 'irrigation1',
'lnfer_tmp_interaction', 'lnfer_tmp_tmp_interaction']
response_variable = 'lnyield'
# Get unique group values
group_values = data['group'].unique()
# Initialize an empty list to store data for each group
group_dfs = []
for group_idx, group in enumerate(group_values):
# Filter data for the current group
group_data = data[data['group'] == group].copy() # Make a copy to avoid warnings
# Select predictor variables based on the group
predictor_variables = predictor_variables1 if group == 1 else predictor_variables2
# Perform regression analysis for the current group
X = sm.add_constant(group_data[predictor_variables])
y = group_data[response_variable]
reg = sm.OLS(y, X).fit()
group_data = data[data['group'] == group].copy()
if group == 1:
group_data['y_temp'] = (reg.params['tmp_tmp_interaction'] * group_data['tmp']*2 +
reg.params['tmp'] +
reg.params['lnfer_irrigation1_tmp_interaction'] * group_data['lnfer'] * group_data['irrigation1'] +
reg.params['lnfer_irrigation12_tmp_tmp_interaction'] * group_data['irrigation12'] * group_data['lnfer'] * group_data['tmp']*2)*100
else:
group_data['y_temp'] = (reg.params['tmp_tmp_interaction'] * group_data['tmp']*2 +
reg.params['tmp'] + reg.params['lnfer_tmp_interaction'] * group_data['lnfer'] +
reg.params['lnfer_tmp_tmp_interaction'] * group_data['lnfer'] * group_data['tmp']*2
)*100
group_dfs.append(group_data)
# Concatenate the dataframes for all groups
data = pd.concat(group_dfs)
# Filter data for the year 2020
year_2020_data = data[data['year'] == 2020]
# Convert the 'geoid' column in 'year_2020_data' to 'object' data type
year_2020_data['geoid'] = year_2020_data['geoid'].astype(str)
# 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(str)
# Rename the 'geoid' column to 'GEOID' in year_2020_data
year_2020_data.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Now, merge 'year_2020_data' with county boundaries based on 'GEOID' attribute
merged_data = gdf_county.merge(year_2020_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]
# Create a custom colormap with white color for 0 values and normalize it to the data range
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.RdYlGn(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-11, vmax=11)
filtered_data.loc[filtered_data['y_temp'] > 11, 'y_temp'] = 11
filtered_data.loc[filtered_data['y_temp'] < -11, 'y_temp'] = 11
# 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 = [-10, -5, 0, 5, 10]
# Adjust the height of the colorbar by setting the shrink parameter
cbar = plt.colorbar(sc, ax=ax, label='Yield change (%)', orientation='vertical', pad=-0.02, shrink=0.8, 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 (%)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16) # Adjust tick label fontsize
filtered_data.plot(column='y_temp', 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.27, 0.22, '(a) Maize', 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('maize_tmp.jpg', format='jpg', dpi=300, bbox_inches='tight', pad_inches=0.1)
# Show the plot
plt.show()
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_5814/68807675.py:60: 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
year_2020_data['geoid'] = year_2020_data['geoid'].astype(str)
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_5814/68807675.py:71: 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
year_2020_data.rename(columns={'geoid': 'GEOID'}, inplace=True)
