参考自：?WRF模拟区域绘制和局地放大 (qq.com)
from netCDF4 import Dataset
import numpy as np
from wrf import getvar,get_cartopy
import matplotlib.pyplot as plt
plt.rcParams['font.family']='Arial'
plt.rcParams['font.size']=12.5
import matplotlib.transforms as mtransforms
import cartopy.crs as ccrs
from cartopy.io.shapereader import BasicReader
from cartopy.mpl.gridliner import LATITUDE_FORMATTER,LONGITUDE_FORMATTER
import cmaps
from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
from mpl_toolkits.axes_grid1.inset_locator import TransformedBbox, BboxPatch, BboxConnector
## 子图连线函数
def mark_inset(parent_axes, inset_axes, loc1a=1, loc1b=1, loc2a=2, loc2b=2, **kwargs):
rect = TransformedBbox(inset_axes.viewLim, parent_axes.transData)
pp = BboxPatch(rect, fill=False, **kwargs)
parent_axes.add_patch(pp)
p1 = BboxConnector(inset_axes.bbox, rect, loc1=loc1a, loc2=loc1b, **kwargs)
inset_axes.add_patch(p1)
p1.set_clip_on(False)
p2 = BboxConnector(inset_axes.bbox, rect, loc1=loc2a, loc2=loc2b, **kwargs)
inset_axes.add_patch(p2)
p2.set_clip_on(False)
return pp, p1, p2
## 等高子图设置
def add_equal_axes(ax, loc, pad, width,projection):
'''
在原有的Axes旁新添一个等高或等宽的Axes并返回该对象.
Parameters
----------
ax : Axes or array_like of Axes
原有的Axes,也可以为一组Axes构成的数组.
loc : {'left', 'right', 'bottom', 'top'}
新Axes相对于旧Axes的位置.
pad : float
新Axes与旧Axes的间距.
width: float
当loc='left'或'right'时,width表示新Axes的宽度.
当loc='bottom'或'top'时,width表示新Axes的高度.
Returns
-------
ax_new : Axes
新Axes对象.
'''
# 无论ax是单个还是一组Axes,获取ax的大小位置.
axes = np.atleast_1d(ax).ravel()
bbox = mtransforms.Bbox.union([ax.get_position() for ax in axes])
# 决定新Axes的大小位置.
if loc == 'left':
x0_new = bbox.x0 - pad - width
x1_new = x0_new + width
y0_new = bbox.y0
y1_new = bbox.y1
elif loc == 'right':
x0_new = bbox.x1 + pad
x1_new = x0_new + width
y0_new = bbox.y0
y1_new = bbox.y1
elif loc == 'bottom':
x0_new = bbox.x0
x1_new = bbox.x1
y0_new = bbox.y0 - pad - width
y1_new = y0_new + width
elif loc == 'top':
x0_new = bbox.x0
x1_new = bbox.x1
y0_new = bbox.y1 + pad
y1_new = y0_new + width
# 创建新Axes.
fig = axes[0].get_figure()
bbox_new = mtransforms.Bbox.from_extents(x0_new, y0_new, x1_new, y1_new)
ax_new = fig.add_axes(bbox_new,projection=projection)
return ax_new
def LU_MODIS21(uniq=np.arange(1, 22)):
inds = uniq-1
C = np.array([
[0,.4,0], # 1 Evergreen Needleleaf Forest
[0,.4,.2], # 2 Evergreen Broadleaf Forest
[.2,.8,.2], # 3 Deciduous Needleleaf Forest
[.2,.8,.4], # 4 Deciduous Broadleaf Forest
[.2,.6,.2], # 5 Mixed Forests
[.3,.7,0], # 6 Closed Shrublands
[.82,.41,.12], # 7 Open Shurblands
[.74,.71,.41], # 8 Woody Savannas
[1,.84,.0], # 9 Savannas
[0,1,0], # 10 Grasslands
[0,1,1], # 11 Permanant Wetlands
[1,1,0], # 12 Croplands
[1,0,0], # 13 Urban and Built-up
[.7,.9,.3], # 14 Cropland/Natual Vegation Mosaic
[1,1,1], # 15 Snow and Ice
[.914,.914,.7], # 16 Barren or Sparsely Vegetated
[.5,.7,1], # 17 Water (like oceans)
[.86,.08,.23], # 18 Wooded Tundra
[.97,.5,.31], # 19 Mixed Tundra
[.91,.59,.48], # 20 Barren Tundra
[0,0,.88]]) # 21 Lake
cm = ListedColormap(C[inds])
labels = ['Evergreen Needleleaf Forest',
'Evergreen Broadleaf Forest',
'Deciduous Needleleaf Forest',
'Deciduous Broadleaf Forest',
'Mixed Forests',
'Closed Shrublands',
'Open Shrublands',
'Woody Savannas',
'Savannas',
'Grasslands',
'Permanent Wetlands',
'Croplands',
'Urban and Built-Up',
'Cropland/Natural Vegetation Mosaic',
'Snow and Ice',
'Barren or Sparsely Vegetated',
'Water',
'Wooded Tundra',
'Mixed Tundra',
'Barren Tundra',
'Lake']
return cm, np.array(labels)[inds]
def ld1(landuse):
# type 1:
cm, labels = LU_MODIS21()
ticks = [x+1.5 for x in range(len(labels))]
n_max = len(labels)
return to_np(landuse), ticks, labels, cm, n_max
def start(file_in_1, file_in_2):
ncfile_1 = Dataset(file_in_1)
ncfile_2 = Dataset(file_in_2)
landuse_1 = getvar(ncfile_1, 'LU_INDEX')
landuse_2 = getvar(ncfile_2, 'LU_INDEX')
lats_1, lons_1 = latlon_coords(landuse_1)
lats_2, lons_2 = latlon_coords(landuse_2)
cart_proj = get_cartopy(landuse_1)
# Create a figure
fig = plt.figure(figsize=(15,15),dpi=900)
# Set the GeoAxes to the projection used by WRF
ax = plt.axes(projection=cart_proj)
# Use the data extent
ax.set_xlim(cartopy_xlim(landuse_1))
ax.set_ylim(cartopy_ylim(landuse_1))
#-------------------------------------------------绘图部分d01------------------------------------------
plt.tight_layout() #自动调整子图参数
landuse_new_1, ticks, labels, cm, n_max = ld1(landuse_1)
#create具有非规则矩形网格的伪色彩图。
im_1 = ax.pcolormesh(to_np(lons_1), to_np(lats_1), landuse_new_1, vmin=1, vmax=n_max+1,
cmap=cm, alpha=0.8, transform=ccrs.PlateCarree())
# #设置图例bar
# cbar = fig.colorbar(im_1, shrink=0.98, pad = 0.08, ax=ax_inset)
# cbar.set_ticks(ticks)
# cbar.ax.set_yticklabels(labels)
# *must* call draw in order to get the axis boundary used to add ticks:
fig.canvas.draw()
# Define gridline locations and draw the lines using cartopy's built-in gridliner:
xticks_1=[80,90,100,110,120] #add xticks
yticks_1=[30,35,40,45,50]
# font3={'family':'SimHei','size':25,'color':'k'}
# Add the gridlines
ax.gridlines(xlocs=xticks_1, ylocs=yticks_1, draw_labels=False, linewidth=0.9, color='grey', alpha=0.6, linestyle='--')
# 设置经纬度和刻度
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)
#设置坐标轴字体大小
ax.text(0.005,0.945,'D01',fontsize=22,fontweight='bold',transform=ax.transAxes)
ax.tick_params(axis="x", labelsize=22)
ax.tick_params(axis="y", labelsize=22)
countries = BasicReader("F:/world_adm0_Project.shp")
provinces = BasicReader("F:/chinaProvince.shp")
provinces1 = BasicReader("F:/gnnq.shp")
#叠加shp
ax.add_geometries(provinces.geometries(), linewidth=1.5, edgecolor='blue', crs=ccrs.PlateCarree(),
facecolor='none')
ax.add_geometries(countries.geometries(), linewidth=1.5, edgecolor='black', crs=ccrs.PlateCarree(),
facecolor='none')
#g1.rotate_labels = False
# ------------------------------在d01的模拟区域上框出d02的模拟区域范围--------------------------------------------------
ax.plot([lon_1[0, 0], lon_1[-1, 0]], [lat_1[0, 0], lat_1[-1, 0]], color='red', lw=2, transform=ccrs.PlateCarree())
ax.plot([lon_1[-1, 0], lon_1[-1, -1]], [lat_1[-1, 0], lat_1[-1, -1]], color='red', lw=2, transform=ccrs.PlateCarree())
ax.plot([lon_1[-1, -1], lon_1[0, -1]], [lat_1[-1, -1], lat_1[0, -1]], color='red', lw=2, transform=ccrs.PlateCarree())
ax.plot([lon_1[0, -1], lon_1[0, 0]], [lat_1[0, -1], lat_1[0, 0]], color='red', lw=2, transform=ccrs.PlateCarree())
#----------------------------------------------绘图部分d02--------------------------------------------------------------
ax_inset = add_equal_axes(ax, loc='right', pad=0.05, width=0.5, projection=cart_proj)
landuse_new_2, ticks, labels, cm, n_max = ld1(landuse_2)
#create具有非规则矩形网格的伪色彩图。
im_2 = ax_inset.pcolormesh(to_np(lons_2), to_np(lats_2), landuse_new_2, vmin=1, vmax=n_max+1,
cmap=cm, alpha=0.8, transform=ccrs.PlateCarree())
#叠加shp数据
fig.canvas.draw()
ax_inset.add_geometries(provinces1.geometries(), linewidth=1, edgecolor='black', crs=ccrs.PlateCarree(),
facecolor='none')
# 添加经纬度网格和刻度
# Define gridline locations and draw the lines using cartopy's built-in gridliner:
# ax_inset.gridlines(xlocs=xticks, ylocs=yticks)
xticks_2=[100,102,104]
yticks_2=[36,37,38,39,40]
# Use the data extent
# ax_inset.set_xlim(cartopy_xlim(landuse_2))
# ax_inset.set_ylim(cartopy_ylim(landuse_2))
ax_inset.gridlines(xlocs=xticks_2, ylocs=yticks_2,
draw_labels=False, linewidth=0.8, color='grey', alpha=0.5, linestyle='--')
# Label the end-points of the gridlines using the custom tick makers:
ax_inset.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax_inset.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_xticks(ax_inset, xticks_2)
lambert_yticks(ax_inset, yticks_2)
ax_inset.text(0.005,0.999,'D02',fontsize=22,fontweight='bold',transform=ax_inset.transAxes)
ax_inset.tick_params(axis="x", labelsize=20)
ax_inset.tick_params(axis="y", labelsize=20)
ax_inset.gridlines(xlocs=xticks_2, ylocs=yticks_2, draw_labels=False, linewidth=0.9, color='grey', alpha=0.6, linestyle='--')
#叠加shp
ax_inset.add_geometries(provinces1.geometries(), linewidth=1.5, edgecolor='blue', crs=ccrs.PlateCarree(),
facecolor='none')
#添加子图连线
mark_inset(ax, ax_inset, loc1a=2, loc1b=1, loc2a=3, loc2b=4, fc="none", ec='k', lw=0.85)
#设置图例bar
cbar = fig.colorbar(im_1, shrink=0.98, pad = 0.05, ax=ax_inset)
cbar.set_ticks(ticks)
cbar.ax.set_yticklabels(labels)
pic_name = "F:/pythonplot/lu_type_modis_1.png"
fig.savefig(pic_name, bbox_inches='tight', dpi=900)
plt.close()
print('土地利用绘制完毕')
def main():
file_in_1 = "F:/pythonplot/geo_em.d01.nc"
file_in_2 = "F:/pythonplot/geo_em.d02.nc"
start(file_in_1, file_in_2)
if __name__ == '__main__':
main()