下垫面类型对于WRF的地表过程十分重要,而在我们研究WRF的地表过程之前,需要对输入的土地利用类型进行一些绘制,以便后续的修改。
土地利用分类图的特点主要在于色标的设置,而在Github中,已有人根据WRF的不同土地利用数据设置了色标与标签,详见:landuse_colormao
在这里,我将以北极地区为例,绘制北极地区的WRF下垫面数据,我使用的MODIS21这类。
首先加载包,并进行函数定义:

## Custom Terrain ColorMaps by Brian
import os
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.path as mpath
import cmaps
import matplotlib.pyplot as plt###引入库包
import numpy as np
import numpy.ma as ma
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from netCDF4 import Dataset
try:
    import pykdtree.kdtree
    _IS_PYKDTREE = True
except ImportError:
    import scipy.spatial
    _IS_PYKDTREE = False
from wrf import getvar, interplevel, vertcross,vinterp, ALL_TIMES, CoordPair, xy_to_ll, ll_to_xy, to_np, get_cartopy, latlon_coords, cartopy_xlim, cartopy_ylim
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import pandas as pd
import datetime as dt
import time
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import numpy as np
# Four Different Land Use Categories:
#   LU_MODIS20
#   LU_MODIS21     includes lake category
#   LU_USGS
#   LU_NLCD


## Land Use Colormap
## ! represents categories not in my Utah domain
def LU_MODIS21(): #
    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)
    
    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, labels

LU_MODIS21函数定义了绘制的土地类型、标签与对应色标,我们使用时,直接引用LU_MODIS21()即可返回。
那么有一个小问题,就是当绘制范围内的数据并不包含所有土地利用类型,我们应该怎么做?
其实很简单,只要将你拥有的土地类型数据提取出来,将原本函数中的labels和C切片,重新定义色标再绘制即可。
以下是实例:
首先读取数据:

geo=Dataset('F:\wrfout\geo_em.d01.nc')
landuse=getvar(geo,'LU_INDEX')

获取当前模拟域内的所有土地类型值,并简单构造对应索引:

landclass=np.unique(to_np(landuse))
landidx=np.unique(to_np(landuse))-1
landidx=landidx.tolist()

重新构造色标标签:

 cm,labels = LU_MODIS21()  
    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
    c=C.take(landidx,0)
    labels=[labels[int(landidx[i])] for i in range(len(landidx))]
    cm=ListedColormap(c)
   

开始绘图,这里使用了本人自行封装的arcticplot函数,具体可参见本人之前的博客:cartopy绘制北极

fig,f1_ax1,x,y,masked_v=arcticplot(z_masked_overlap,lat2d,lon2d,to_np(landuse))
    font = {'family' : 'serif',
             'color'  : 'darkred',
             'weight' : 'normal',
             'size'   : 16,
             }
   
c7=f1_ax1.pcolormesh(x, y, to_np(landuse), cmap=cm)
cbar=fig.colorbar(c7,shrink=1)#设置图标
cbar.set_ticks(landclass)
cbar.ax.set_yticklabels(labels)
 #cbar.ax.set_xticklabels(labels) #If using a horizontal colorbar orientation
cbar.ax.invert_yaxis()
cbar.ax.tick_params(labelsize=15)
plt.show()

绘图如下:
在这里插入图片描述

Logo

开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!

更多推荐