【wrf-python】wrf domain 嵌套土地利用类型绘制

2023-12-15 10:36:49

参考自:?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()
    

文章来源:https://blog.csdn.net/A18040554844/article/details/135010359
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。