In [1]:
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 (
In [2]:
os.chdir('/Users/chenchenren/postdoc/paper/2N and water-US/Figure 7 Cost and Benefit/')
In [3]:
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)
In [4]:
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
In [5]:
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)
In [6]:
# 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)
In [7]:
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
In [8]:
# 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)
In [9]:
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
In [10]:
# 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)
In [11]:
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
In [12]:
# 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)
In [13]:
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
In [21]:
# 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)
In [15]:
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
In [20]:
# 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)
In [23]:
# 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)
In [ ]: