有人能告诉我如何修改下面的脚本来使用Geopandas添加美国的州边界吗?我试过了,但似乎不能让它工作。此外,如果有一个特定的国家边界shapefile,任何人都可以推荐,请让我知道。否则,这个图看起来和我想要的完全一样。注意,我不想使用map.drawstates()因为它也能画出其他国家的州。我只想要美国的州界。
from netCDF4 import Dataset as NetCDFFile
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as mcolors
import geopandas as gpd
# Read the U.S. state borders shapefile
us_states = gpd.read_file('C:/Users/scapa/Downloads/tl_2012_us_state/tl_2012_us_state.shp')
nc = NetCDFFile('C:/Users/scapa/Downloads/nclwAXbfFo6UE.nc')
lat = nc.variables['lat'][:]
lon = nc.variables['lon'][:]
prate = nc.variables['VAR'][:]
# Determine the extent of the shapefile
state_bounds = us_states.total_bounds # Get the bounding box of the state borders
# Calculate the map extent
map_extent = [state_bounds[0], state_bounds[2], state_bounds[1], state_bounds[3]]
fig, ax = plt.subplots(figsize=(10, 8))
map = Basemap(llcrnrlon=230., llcrnrlat=10., urcrnrlon=305., urcrnrlat=55.)
lons, lats = np.meshgrid(lon, lat)
x, y = map(lons, lats)
# Set the levels and colors
levels = np.arange(-3.6, 3.6, 0.4)
cmap = plt.get_cmap('BrBG', len(levels) - 1)
# Create a custom colormap with white for values between -0.4 and 0.4
cmaplist = [(1, 1, 1, 1) if -0.4 <= val <= 0.4 else cmap((val - levels.min()) / (levels.max() - levels.min())) for val in levels]
custom_cmap = mcolors.LinearSegmentedColormap.from_list('custom', cmaplist, len(levels))
# Add U.S. state borders
us_states.boundary.plot(ax=map, linewidth=0.5)
cs = map.contourf(x, y, prate, levels, cmap=custom_cmap)
cbar = plt.colorbar(cs, orientation='horizontal', ticks=levels, pad=0.07)
cbar.set_label('Surface Precipitation Rate Anomalies (mm/day)')
cbar.ax.tick_params(labelsize=8)
map.drawcoastlines()
map.drawparallels(np.arange(10, 60, 10), labels=[1, 1, 1, 1], linewidth=0.5)
map.drawmeridians(np.arange(240, 310, 20), labels=[1, 1, 0, 1], linewidth=0.5)
plt.title('El Niño Winter Precipitation Anomalies', fontsize=15, y=1.05)
plt.text(0.5, 1.02, 'Winters of 1957-58, 1965-66, 1972-73, 1982-83, 1987-88, 1991-92, 1997-98, 2009-10, 2015-16', horizontalalignment='center', fontsize=9, fontstyle='italic',transform=ax.transAxes)
plt.savefig('era5.precip.anom.states.png', dpi=600)
plt.show()
字符串
的数据
1条答案
按热度按时间dz6r00yl1#
事实证明,Cartopy是解决我的问题的最佳方案。下面是上面代码的一些关键添加/更改:
字符串