In [5]:
data = alldata[alldata['scenario'] == 1.5]
# 5% percentile
data['gap_production_warming'] = data.groupby(['geoid'])['gap_production_warming'].transform(lambda x: np.percentile(x, 5))
data=data.drop_duplicates(subset=['geoid'])
print(data)
print(np.min(data['gap_production_warming']),np.max(data['gap_production_warming']),np.mean(data['gap_production_warming']))
print(np.percentile(data['gap_production_warming'], 5),np.percentile(data['gap_production_warming'], 95))
# Replace values in 'gap_yield' greater than p90 with p90
data.loc[data['gap_production_warming'] > 2.1, 'gap_production_warming'] = 2.1
data.loc[data['gap_production_warming'] < -2.1, 'gap_production_warming'] = -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}$ t)', 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('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_warming', 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_p5.png', format='png',
dpi=300, bbox_inches='tight', pad_inches=0.1, transparent=True)
# Show the plot
plt.show()
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_61041/2113936764.py:3: 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_warming'] = data.groupby(['geoid'])['gap_production_warming'].transform(lambda x: np.percentile(x, 5))
geoid scenario id_cof tmp pre lnfer irrigation1 \
0 1001 1.5 1 25.099590 6.254395 2.457656 0.377088
2000 1003 1.5 1 26.261318 9.274925 1.292553 0.738441
4000 1005 1.5 1 25.408514 6.727813 2.114730 0.028216
6000 1009 1.5 1 23.976584 6.790346 4.361742 0.056342
8000 1011 1.5 1 25.246212 6.448847 4.858469 0.030586
... ... ... ... ... ... ... ...
3536000 55133 1.5 1 18.514920 5.639188 1.516143 NaN
3538000 55135 1.5 1 17.319042 5.592127 1.466300 NaN
3540000 55137 1.5 1 17.765668 5.555818 1.553310 NaN
3542000 55139 1.5 1 18.099080 5.372444 1.255133 NaN
3544000 55141 1.5 1 17.398562 5.719637 2.556201 NaN
y_temp_base_t y_temp_base2015 lnyield_obs har \
0 7.474178 7.499539 7.509678 779.625
2000 7.197338 7.207396 7.927428 6358.433
4000 7.428942 7.443278 8.161582 526.500
6000 7.889409 7.925261 7.840701 1107.736
8000 7.808706 7.827184 7.984614 870.750
... ... ... ... ...
3536000 7.954243 7.917639 7.974899 7837.872
3538000 7.876120 7.832488 7.957613 9909.051
3540000 7.910403 7.871683 7.941248 5335.875
3542000 7.894186 7.865316 8.028848 14541.420
3544000 7.966143 7.894034 7.921128 6700.507
gap_yield_warming gap_production_warming
0 -0.045718 -0.003727
2000 -0.027743 -0.019340
4000 -0.049872 -0.002605
6000 -0.089522 -0.009786
8000 -0.053741 -0.004705
... ... ...
3536000 0.108381 0.083619
3538000 0.127429 0.116701
3540000 0.110971 0.055787
3542000 0.089869 0.121146
3544000 0.205987 0.135636
[1773 rows x 13 columns]
-3.3181771942742713 3.053110777122164 0.005905536688759566
-0.43960474399223093 0.5036116026289954
