WRF后处理:python cartopy绘制土地利用/土地分类图//python绘制WRF下垫面类型(以北极为例)
下垫面类型对于WRF的地表过程十分重要,而在我们研究WRF的地表过程之前,需要对输入的土地利用类型进行一些绘制,以便后续的修改。LU_MODIS21函数定义了绘制的土地类型、标签与对应色标,我们使用时,直接引用LU_MODIS21()即可返回。其实很简单,只要将你拥有的土地类型数据提取出来,将原本函数中的labels和C切片,重新定义色标再绘制即可。在这里,我将以北极地区为例,绘制北极地区的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
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()
绘图如下:
开放原子开发者工作坊旨在鼓励更多人参与开源活动,与志同道合的开发者们相互交流开发经验、分享开发心得、获取前沿技术趋势。工作坊有多种形式的开发者活动,如meetup、训练营等,主打技术交流,干货满满,真诚地邀请各位开发者共同参与!
更多推荐
所有评论(0)