WRF后处理：python cartopy绘制土地利用/土地分类图//python绘制WRF下垫面类型（以北极为例）

下垫面类型对于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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87

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

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

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

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

重新构造色标标签：

 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)
   
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

开始绘图，这里使用了本人自行封装的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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

绘图如下：