import pandas as pd
import numpy as np
import statsmodels.api as sm
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 tifffile
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
import matplotlib.lines as mlines
from geotiff import GeoTiff
from rasterio import open as rasopen
import matplotlib.cm as cm
from PIL import Image
from matplotlib.font_manager import FontProperties
/Users/chenchenren/anaconda3/lib/python3.11/site-packages/pandas/core/arrays/masked.py:60: UserWarning: Pandas requires version '1.3.6' or newer of 'bottleneck' (version '1.3.5' currently installed). from pandas.core import (
os.chdir('/Users/chenchenren/postdoc/paper/2N and water-US/Figure 7 Cost and Benefit/')
excel_file = "maize_cost_benefit_county.csv"
data = pd.read_csv(excel_file)
data15 = data[data['scenario'] == 1.5]
data15['geoid'] = data15['geoid'].astype(float)
print(data15)
data3 = data[data['scenario'] == 3.0]
data3['geoid'] = data3['geoid'].astype(float)
print(data3)
geoid scenario production_benefit_irri production_benefit_lnfer \ 1 48393.0 1.5 -0.251293 0.170940 3 48195.0 1.5 -1.242042 4.431190 5 48233.0 1.5 -1.542894 1.454517 7 48163.0 1.5 4.037799 0.732635 9 48129.0 1.5 0.142622 0.132830 ... ... ... ... ... 4181 53021.0 1.5 0.101028 0.100647 4182 53005.0 1.5 0.060051 0.060051 4184 41059.0 1.5 0.000352 0.000561 4186 41045.0 1.5 0.430206 0.430206 4188 41049.0 1.5 0.076315 0.076315 production_benefit_mana irri_cost_irri irri_cost_mana fer_cost_fer \ 1 0.184301 0.000000 0.000000 0.205265 3 4.829888 1.522186 0.000000 3.092772 5 1.604419 0.000000 0.000000 1.804322 7 0.842611 0.232803 0.000000 1.988732 9 0.142632 0.014755 0.000000 0.038230 ... ... ... ... ... 4181 0.101065 0.013235 -0.561538 0.056264 4182 0.060051 0.000000 0.000000 0.000000 4184 0.000561 0.000029 -0.004488 0.000469 4186 0.430206 0.000000 0.000000 0.000000 4188 0.076315 0.000000 0.000000 0.000000 fer_cost_mana net_benefits_irri net_benefits_fer net_benefits_mana \ 1 0.205265 -0.251293 -0.034325 -0.020963 3 3.092772 -2.764228 1.338418 1.737116 5 1.804322 -1.542894 -0.349805 -0.199903 7 1.988732 3.804996 -1.256097 -1.146121 9 0.038230 0.127867 0.094601 0.104403 ... ... ... ... ... 4181 3.656216 0.087793 0.044383 -2.993612 4182 0.000000 0.060051 0.060051 0.060051 4184 0.046783 0.000323 0.000093 -0.041733 4186 0.000000 0.430206 0.430206 0.430206 4188 0.000000 0.076315 0.076315 0.076315 state1 production_obs 1 Texas 1.271143e+07 3 Texas 2.893905e+08 5 Texas 8.392635e+07 7 Texas 3.026130e+07 9 Texas 1.042043e+07 ... ... ... 4181 Washington 1.248886e+08 4182 Washington 1.089175e+08 4184 Oregon 3.779036e+07 4186 Oregon 8.652065e+07 4188 Oregon 3.742098e+07 [2095 rows x 14 columns] geoid scenario production_benefit_irri production_benefit_lnfer \ 0 48393.0 3.0 0.089430 0.568409 2 48195.0 3.0 5.779283 12.817246 4 48233.0 3.0 0.546965 3.889810 6 48163.0 3.0 6.631563 1.920470 8 48129.0 3.0 0.517039 0.407964 ... ... ... ... ... 4180 53021.0 3.0 0.143144 0.143453 4183 53005.0 3.0 0.089902 0.089902 4185 41059.0 3.0 0.016020 0.016020 4187 41045.0 3.0 0.361436 0.361436 4189 41049.0 3.0 0.162405 0.162405 production_benefit_mana irri_cost_irri irri_cost_mana fer_cost_fer \ 0 0.748957 0.000000 0.000000 0.937968 2 16.837568 1.522186 0.000000 12.006741 4 5.185463 0.000000 0.000000 6.506429 6 2.919224 0.232803 0.000000 9.964947 8 0.517105 0.032833 0.000000 0.156157 ... ... ... ... ... 4180 0.144304 0.017942 -0.561538 0.080147 4183 0.089902 0.000000 0.000000 0.000000 4185 0.016020 0.000000 0.000000 0.000000 4187 0.361436 0.000000 0.000000 0.000000 4189 0.162405 0.000000 0.000000 0.000000 fer_cost_mana net_benefits_irri net_benefits_fer net_benefits_mana \ 0 0.937968 0.089430 -0.369560 -0.189011 2 12.006741 4.257097 0.810505 4.830827 4 6.506429 0.546965 -2.616619 -1.320966 6 9.964947 6.398760 -8.044477 -7.045723 8 0.156157 0.484206 0.251807 0.360948 ... ... ... ... ... 4180 3.888411 0.125201 0.063306 -3.182569 4183 0.000000 0.089902 0.089902 0.089902 4185 0.000000 0.016020 0.016020 0.016020 4187 0.000000 0.361436 0.361436 0.361436 4189 0.000000 0.162405 0.162405 0.162405 state1 production_obs 0 Texas 1.271143e+07 2 Texas 2.893905e+08 4 Texas 8.392635e+07 6 Texas 3.026130e+07 8 Texas 1.042043e+07 ... ... ... 4180 Washington 1.248886e+08 4183 Washington 1.089175e+08 4185 Oregon 3.779036e+07 4187 Oregon 8.652065e+07 4189 Oregon 3.742098e+07 [2095 rows x 14 columns]
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/1149997662.py:5: 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 data15['geoid'] = data15['geoid'].astype(float) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/1149997662.py:9: 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 data3['geoid'] = data3['geoid'].astype(float)
positive_irri=data15[data15['net_benefits_irri']>=0]
print(np.sum(positive_irri['net_benefits_irri']))
positive_fer=data15[data15['net_benefits_fer']>=0]
print(np.sum(positive_fer['net_benefits_fer']))
positive_mana=data15[data15['net_benefits_mana']>=0]
print(np.sum(positive_mana['net_benefits_mana']))
positive_irri=data3[data3['net_benefits_irri']>=0]
print(np.sum(positive_irri['net_benefits_irri']))
positive_fer=data3[data3['net_benefits_fer']>=0]
print(np.sum(positive_fer['net_benefits_fer']))
positive_mana=data3[data3['net_benefits_mana']>=0]
print(np.sum(positive_mana['net_benefits_mana']))
1099.500194206254 923.8257195169043 1132.558445295366 1548.1302568316698 1203.1421144872845 1519.6743318418678
print(np.min(data15['net_benefits_irri']),np.max(data15['net_benefits_irri']),np.mean(data15['net_benefits_irri']))
print(np.percentile(data15['net_benefits_irri'], 5),np.percentile(data15['net_benefits_irri'], 95))
p95 = np.percentile(data15['net_benefits_irri'], 95)
data15.loc[(data15['net_benefits_irri'] > 4), 'net_benefits_irri'] = 4
data15.loc[(data15['net_benefits_irri'] < -4), 'net_benefits_irri'] = -4
data15['geoid'] = data15['geoid'].astype(float)
-66.80826169726707 31.622774542706605 -0.44507942642921616 -3.790673424685539 2.383502395564208
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/1154970319.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 data15['geoid'] = data15['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' data15 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data15' Data15Frame
data15.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data15' Data15Frame to 'object' data15 type
data15['GEOID'] = data15['GEOID'].astype(float)
# Merge 'data15' with county boundaries based on 'GEOID' attribute
merged_data15 = gdf_county.merge(data15, on='GEOID', how='left')
# Filter the data15 to include only geometries within the specified longitude and latitude range
filtered_data15 = merged_data15.cx[-130:-60, 25:50]
# Plot the county boundaries
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')
# Create a ScalarMappable without plotting the data15
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)
cbar.ax.set_ylabel('Net benefits (Million $)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
# Plot the mean_gap_irrigation values
filtered_data15.plot(column='net_benefits_irri', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=16)
font_prop.set_weight('bold')
# Additional plot configurations
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.28, 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_profit_irri_15.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_12671/2558518357.py:9: 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 data15.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/2558518357.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 data15['GEOID'] = data15['GEOID'].astype(float)
print(np.min(data3['net_benefits_irri']),np.max(data3['net_benefits_irri']),np.mean(data3['net_benefits_irri']))
print(np.percentile(data3['net_benefits_irri'], 5),np.percentile(data3['net_benefits_irri'], 95))
p95 = np.percentile(data3['net_benefits_irri'], 95)
data3.loc[(data3['net_benefits_irri'] > 4), 'net_benefits_irri'] = 4
data3.loc[(data3['net_benefits_irri'] < -4), 'net_benefits_irri'] = -4
-65.73237750743054 43.81009787180551 -0.5582196616058241 -6.366208228430618 2.876500034682916
# 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' data3 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data3' Data3Frame
data3.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data3' Data3Frame to 'object' data3 type
data3['GEOID'] = data3['GEOID'].astype(float)
# Merge 'data3' with county boundaries based on 'GEOID' attribute
merged_data3 = gdf_county.merge(data3, on='GEOID', how='left')
# Filter the data3 to include only geometries within the specified longitude and latitude range
filtered_data3 = merged_data3.cx[-130:-60, 25:50]
# Plot the county boundaries
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')
# Create a ScalarMappable without plotting the data3
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)
cbar.ax.set_ylabel('Net benefits (Million $)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
# Plot the mean_gap_irrigation values
filtered_data3.plot(column='net_benefits_irri', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=16)
font_prop.set_weight('bold')
# Additional plot configurations
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.28, 0.22, '(g) 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_profit_irri_3.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_12671/2348013683.py:9: 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 data3.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/2348013683.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 data3['GEOID'] = data3['GEOID'].astype(float)
print(np.min(data15['net_benefits_fer']),np.max(data15['net_benefits_fer']),np.mean(data15['net_benefits_fer']))
print(np.percentile(data15['net_benefits_fer'], 5),np.percentile(data15['net_benefits_fer'], 95))
p95 = np.percentile(data15['net_benefits_fer'], 95)
data15.loc[(data15['net_benefits_fer'] > 4), 'net_benefits_fer'] = 4
data15.loc[(data15['net_benefits_fer'] < -4), 'net_benefits_fer'] = -4
-3.842688271223061 10.89293678486294 0.4090441391634506 -0.16621207528483145 2.395941635927279
# 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' data15 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data15' Data15Frame
data15.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data15' Data15Frame to 'object' data15 type
data15['GEOID'] = data15['GEOID'].astype(float)
# Merge 'data15' with county boundaries based on 'GEOID' attribute
merged_data15 = gdf_county.merge(data15, on='GEOID', how='left')
# Filter the data15 to include only geometries within the specified longitude and latitude range
filtered_data15 = merged_data15.cx[-130:-60, 25:50]
# Plot the county boundaries
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')
# Create a ScalarMappable without plotting the data15
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)
cbar.ax.set_ylabel('Net benefits (Million $)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
# Plot the mean_gap_fergation values
filtered_data15.plot(column='net_benefits_fer', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=16)
font_prop.set_weight('bold')
# Additional plot configurations
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.28, 0.22, '(b) 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_profit_fer_15.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_12671/825270477.py:9: 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 data15.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/825270477.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 data15['GEOID'] = data15['GEOID'].astype(float)
print(np.min(data3['net_benefits_fer']),np.max(data3['net_benefits_fer']),np.mean(data3['net_benefits_fer']))
print(np.percentile(data3['net_benefits_fer'], 5),np.percentile(data3['net_benefits_fer'], 95))
p95 = np.percentile(data3['net_benefits_fer'], 95)
data3.loc[(data3['net_benefits_fer'] > 4), 'net_benefits_fer'] = 4
data3.loc[(data3['net_benefits_fer'] < -4), 'net_benefits_fer'] = -4
-23.43803642654318 12.858336180578352 0.4078927623917264 -0.7532070938520621 3.3531496412994906
# 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' data3 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data3' Data3Frame
data3.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data3' Data3Frame to 'object' data3 type
data3['GEOID'] = data3['GEOID'].astype(float)
# Merge 'data3' with county boundaries based on 'GEOID' attribute
merged_data3 = gdf_county.merge(data3, on='GEOID', how='left')
# Filter the data3 to include only geometries within the specified longitude and latitude range
filtered_data3 = merged_data3.cx[-130:-60, 25:50]
# Plot the county boundaries
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')
# Create a ScalarMappable without plotting the data3
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)
cbar.ax.set_ylabel('Net benefits (Million $)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
# Plot the mean_gap_fergation values
filtered_data3.plot(column='net_benefits_fer', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=16)
font_prop.set_weight('bold')
# Additional plot configurations
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.28, 0.22, '(h) 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_profit_fer_3.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_12671/825294603.py:9: 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 data3.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/825294603.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 data3['GEOID'] = data3['GEOID'].astype(float)
print(np.min(data15['net_benefits_mana']),np.max(data15['net_benefits_mana']),np.mean(data15['net_benefits_mana']))
print(np.percentile(data15['net_benefits_mana'], 5),np.percentile(data15['net_benefits_mana'], 95))
p95 = np.percentile(data15['net_benefits_mana'], 95)
data15.loc[(data15['net_benefits_mana'] > 4), 'net_benefits_mana'] = 4
data15.loc[(data15['net_benefits_mana'] < -4), 'net_benefits_mana'] = -4
-35.47509712649492 54.41902102067451 -0.010490560809896953 -3.625147354166761 2.6448254312054247
# 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' data15 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data15' Data15Frame
data15.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data15' Data15Frame to 'object' data15 type
data15['GEOID'] = data15['GEOID'].astype(float)
# Merge 'data15' with county boundaries based on 'GEOID' attribute
merged_data15 = gdf_county.merge(data15, on='GEOID', how='left')
# Filter the data15 to include only geometries within the specified longitude and latitude range
filtered_data15 = merged_data15.cx[-130:-60, 25:50]
# Plot the county boundaries
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')
# Create a ScalarMappable without plotting the data15
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)
cbar.ax.set_ylabel('Net benefits (Million $)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
# Plot the mean_gap_managation values
filtered_data15.plot(column='net_benefits_mana', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=16)
font_prop.set_weight('bold')
# Additional plot configurations
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.15, 0.15, '(a) Maize under\n1.5°C warming', transform=plt.gcf().transFigure,
fontsize=24, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')
# Save the entire figure as a JPG with DPI 300
plt.savefig('maize_profit_mana_15.svg', format='svg',
dpi=300, bbox_inches='tight', pad_inches=0.1, transparent=True)
# Show the plot
plt.show()
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/2281804117.py:9: 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 data15.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/2281804117.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 data15['GEOID'] = data15['GEOID'].astype(float)
print(np.min(data3['net_benefits_mana']),np.max(data3['net_benefits_mana']),np.mean(data3['net_benefits_mana']))
print(np.percentile(data3['net_benefits_mana'], 5),np.percentile(data3['net_benefits_mana'], 95))
p95 = np.percentile(data3['net_benefits_mana'], 95)
data3.loc[(data3['net_benefits_mana'] > 4), 'net_benefits_mana'] = 4
data3.loc[(data3['net_benefits_mana'] < -4), 'net_benefits_mana'] = -4
-37.66545605638162 43.76440028674214 -0.04210149335560882 -4.115095414596415 3.827956170396017
# 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' data3 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data3' Data3Frame
data3.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data3' Data3Frame to 'object' data3 type
data3['GEOID'] = data3['GEOID'].astype(float)
# Merge 'data3' with county boundaries based on 'GEOID' attribute
merged_data3 = gdf_county.merge(data3, on='GEOID', how='left')
# Filter the data3 to include only geometries within the specified longitude and latitude range
filtered_data3 = merged_data3.cx[-130:-60, 25:50]
# Plot the county boundaries
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')
# Create a ScalarMappable without plotting the data3
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $)', orientation='vertical', pad=0.02, shrink=1, ticks=custom_ticks)
cbar.ax.set_ylabel('Net benefits (Million $)', fontsize=16, fontfamily='Arial')
cbar.ax.tick_params(labelsize=16)
cbar.ax.set_yticklabels(custom_ticks, fontsize=16, fontfamily='Arial')
# Plot the mean_gap_managation values
filtered_data3.plot(column='net_benefits_mana', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=16)
font_prop.set_weight('bold')
# Additional plot configurations
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.15, 0.15, '(b) Maize under\n3°C warming', transform=plt.gcf().transFigure,
fontsize=24, fontweight='bold', fontfamily='Arial', ha='left', va='bottom')
# Save the entire figure as a JPG with DPI 300
plt.savefig('maize_profit_mana_3.svg', format='svg',
dpi=300, bbox_inches='tight', pad_inches=0.1, transparent=True)
# Show the plot
plt.show()
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/1786408533.py:9: 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 data3.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/1786408533.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 data3['GEOID'] = data3['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' data3 type
gdf_county['GEOID'] = gdf_county['GEOID'].astype(float)
# Rename the 'geoid' column to 'GEOID' in the 'data3' Data3Frame
data3.rename(columns={'geoid': 'GEOID'}, inplace=True)
# Convert the 'GEOID' column in the 'data3' Data3Frame to 'object' data3 type
data3['GEOID'] = data3['GEOID'].astype(float)
# Merge 'data3' with county boundaries based on 'GEOID' attribute
merged_data3 = gdf_county.merge(data3, on='GEOID', how='left')
# Filter the data3 to include only geometries within the specified longitude and latitude range
filtered_data3 = merged_data3.cx[-130:-60, 25:50]
# Plot the county boundaries
fig, ax = plt.subplots(figsize=(20, 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')
# Create a ScalarMappable without plotting the data3
cmap = mcolors.ListedColormap(['white'] + list(plt.cm.PiYG(np.linspace(0, 1, 255))))
norm = mcolors.Normalize(vmin=-4.2, vmax=4.2)
sc = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sc.set_array([])
# Add colorbar
custom_ticks = [-4, -2, 0, 2, 4]
cbar = plt.colorbar(sc, ax=ax, label='Net benefits (Million $/yr)', orientation='horizontal', pad=0.02, shrink=0.8, ticks=custom_ticks)
cbar.ax.set_xlabel('Cost-competitiveness (Million $/yr)', fontsize=40, fontfamily='Arial')
cbar.ax.tick_params(labelsize=40) # Adjust tick label fontsize
cbar.set_ticks(custom_ticks)
# Plot the mean_gap_managation values
filtered_data3.plot(column='net_benefits_mana', cmap=cmap, ax=ax, norm=norm)
# Set font properties for the legend
font_prop = FontProperties(family='Arial', size=40)
font_prop.set_weight('bold')
# Additional plot configurations
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('profit_legend.svg', format='svg',
dpi=300, bbox_inches='tight', pad_inches=0.1, transparent=True)
# Show the plot
plt.show()
/var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/3453793490.py:9: 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 data3.rename(columns={'geoid': 'GEOID'}, inplace=True) /var/folders/vd/0_phd7hx2n51y4412862zww00000gp/T/ipykernel_12671/3453793490.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 data3['GEOID'] = data3['GEOID'].astype(float)