Cartopy畫地圖第八天（冷空氣南下，NCL色標(biāo)使用）

相信很多朋友都看過這張漂亮的圖吧

這是中國氣象愛好者公眾號(hào)畫的一張反映冷空氣的圖片，當(dāng)時(shí)就被驚艷到了，本帖就來復(fù)刻一下這張圖。
本次的知識(shí)點(diǎn)大概總結(jié)一下就是：
1.畫風(fēng)的箭頭流場(chǎng)
2.單畫某一根等值線，比如588，0℃
3.引入NCL色標(biāo)
4.制作gif圖片
先上個(gè)效果圖

一、分析圖片

中國氣象愛好者的圖片拿到手以后開始分析
1.500的位勢(shì)高度畫等值線
2.850的溫度填色
3.0℃等值線
4.588畫等值線
5.850風(fēng)場(chǎng)的箭頭流場(chǎng)
6.注意投影方式
其實(shí)也不難，一個(gè)一個(gè)畫就好了。

二、函數(shù)準(zhǔn)備

（一）、使用NLC色標(biāo)的函數(shù)

import re import matplotlib.cm from matplotlib import colors import matplotlib.pyplot as plt import numpy as np class Colormap(colors.ListedColormap):def __init__(self, c, name='from_list', n=None):self.colors = cself.name = nameself.N = nsuper(Colormap, self).__init__(self.colors, name=self.name, N=self.N)def __getitem__(self, item):return Colormap(self.colors[item], name='sliced_' + self.name)def show(self):a = np.outer(np.ones(10), np.arange(0, 1, 0.001))plt.figure(figsize=(2.5, 0.5))plt.subplots_adjust(top=0.95, bottom=0.05, left=0.01, right=0.99)plt.subplot(111)plt.axis('off')plt.imshow(a, aspect='auto', cmap=self, origin='lower')plt.text(0.5, 0.5, self._name,verticalalignment='center', horizontalalignment='center',fontsize=12, transform=plt.gca().transAxes)plt.show() def coltbl(cmap_file):pattern = re.compile(r'(\d\.?\d*)\s+(\d\.?\d*)\s+(\d\.?\d*).*')with open(cmap_file) as cmap:cmap_buff = cmap.read()cmap_buff = re.compile('ncolors.*\n').sub('', cmap_buff)if re.search(r'\s*\d\.\d*', cmap_buff):return np.asarray(pattern.findall(cmap_buff), 'f4')else:return np.asarray(pattern.findall(cmap_buff), 'u1') / 255. def get_my_cmaps(cname):try:if cname in matplotlib.cm._cmap_registry:return matplotlib.cm.get_cmap(cname)except:passcmap_file = os.path.join( cname+ ".rgb")cmap = Colormap(coltbl(cmap_file), name=cname)matplotlib.cm.register_cmap(name=cname, cmap=cmap)return cmap

一個(gè)類，倆函數(shù)，大家可以研究一下怎么實(shí)現(xiàn)引入NCL色標(biāo)的，這里因?yàn)槠年P(guān)系，就不介紹了，有興趣的朋友可以留言，下次單獨(dú)出一期來介紹引入NCL色標(biāo)的原理。

（二）、數(shù)據(jù)讀取的函數(shù)

本次數(shù)據(jù)使用的是Micaps導(dǎo)出的二進(jìn)制歐細(xì)資料，讀取函數(shù)如下

from read_mdfs import MDFS_Grid from scipy.interpolate import interpolate def get_lat_lon(filepath):#獲取經(jīng)緯度數(shù)據(jù)a = MDFS_Grid(filepath)lon = a.data['Lon']lat = a.data['Lat'][::-1]#翻了緯度lon_scipy = np.arange(lon.min(), lon.max(), 0.05)#插值lat_scipy = np.arange(lat.min(), lat.max(), 0.05)return lon,lat,lon_scipy,lat_scipy def read_data_hgt_or_tmp(filepath,lon,lat,lon_scipy,lat_scipy):#讀位勢(shì)高度和溫度a = MDFS_Grid(filepath)data = a.data['Grid'][::8,::8]#為了數(shù)據(jù)更平滑，降低了精度spline = interpolate.RectBivariateSpline(lat[::8], lon[::8],data,)data = spline(lat_scipy,lon_scipy)#插值return data[::-1,:]#翻了緯度 def read_data(filepath):#讀取一般數(shù)據(jù)，這里是風(fēng)a = MDFS_Grid(filepath)data = a.data['Grid']return data[::-1,:]#翻了緯度

（三）、地圖數(shù)據(jù)讀取

后面畫圖為了做gif圖片，所以要畫好幾次地圖，但是讀取地圖數(shù)據(jù)可以一次完成，想辦法減少運(yùn)行時(shí)間是個(gè)好習(xí)慣哦。

with open('CN-border-La.dat') as src:context = src.read()blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

（四）、其他準(zhǔn)備

fig = plt.figure(figsize=(47, 30), dpi=30)#畫布的設(shè)置，也是一次就夠了，畫完一張圖清空一次就好了，減少運(yùn)行內(nèi)存。 frames = []#用來存放圖片的列表，為的是制作gif TIMES = ['000','003','006','009','012']#做循環(huán)用的時(shí)次

三、循環(huán)畫圖

（一）、讀取數(shù)據(jù)

for TIME in TIMES:#讀取經(jīng)緯度、500位勢(shì)高度file_path = r"ECMWF_HR\ECMWF_HR_HGT_500_21122420." + TIMElon ,lat,lon_scipy,lat_scipy = get_lat_lon(file_path)hgt_500 = read_data_hgt_or_tmp(file_path,lon ,lat,lon_scipy,lat_scipy)#讀取850溫度file_path = r"ECMWF_HR\ECMWF_HR_TMP_850_21122420." + TIMEtem_850 = read_data_hgt_or_tmp(file_path,lon ,lat,lon_scipy,lat_scipy)#讀取850U風(fēng)file_path = r"ECMWF_HR\ECMWF_HR_UGRD_850_21122420." + TIMEu_850 = read_data(file_path)#讀取850V風(fēng)file_path = r"ECMWF_HR\ECMWF_HR_VGRD_850_21122420." + TIMEv_850 = read_data(file_path)

（二）、畫圖和地圖設(shè)置

#加子圖，做一些地圖設(shè)置ax = fig.add_axes([0.02, 0.05, 0.9, 0.9], projection=ccrs.LambertConformal(central_latitude=90, central_longitude=105))extent = [80,130,20,55]ax.set_extent(extent,crs=ccrs.Geodetic())gl = ax.gridlines( draw_labels=True, linewidth=2, color='k', alpha=0.5, linestyle='--')gl.xformatter = LONGITUDE_FORMATTER ##坐標(biāo)刻度轉(zhuǎn)換為經(jīng)緯度樣式gl.yformatter = LATITUDE_FORMATTERgl.xlabel_style = {'size': 30}gl.ylabel_style = {'size': 30}resolution_map = '50m'ax.add_feature(cfeature.OCEAN.with_scale(resolution_map))ax.add_feature(cfeature.LAND.with_scale(resolution_map))ax.add_feature(cfeature.RIVERS.with_scale(resolution_map))ax.add_feature(cfeature.LAKES.with_scale(resolution_map))for line in borders:ax.plot(line[0::2], line[1::2], '-', color='k',transform=ccrs.Geodetic())

這些也不介紹了，在之前的帖子有介紹，可以往前學(xué)習(xí)。

（三）、500位勢(shì)高度和588、0℃線繪制

#畫500位勢(shì)高度ct = ax.contour(lon_scipy, lat_scipy, hgt_500, 20, colors='k',linewidths=3, transform=ccrs.PlateCarree())ax.clabel(ct, inline=True, fontsize=30, fmt='%d')#單畫一遍588，紅色clev = [588]ct = ax.contour(lon_scipy, lat_scipy, hgt_500, clev, colors='r',linewidths=3, transform=ccrs.PlateCarree())ax.clabel(ct, inline=True, fontsize=30, fmt='%d')#畫0℃線，白色clev = [0.0]ct = ax.contour(lon_scipy, lat_scipy, tem_850, clev, colors='w',linewidths=5, transform=ccrs.PlateCarree())ax.clabel(ct, inline=True, fontsize=30, fmt='%d')

（四）、850溫度繪制

#畫850溫度clevs = range(-45,46,3)#引入NCL色標(biāo)my_cmap = get_my_cmaps("MPL_hsv")#翻轉(zhuǎn)色標(biāo)my_cmap = my_cmap.reversed()#填色cf = ax.contourf(lon_scipy, lat_scipy, tem_850, clevs, transform=ccrs.PlateCarree(), cmap=my_cmap)#設(shè)置色標(biāo)position = fig.add_axes([0.92, 0.05, 0.02, 0.9])cb = fig.colorbar(cf, cax=position, orientation='vertical',ticks=range(-45,46,3))cb.set_label('Temperature ℃', fontdict={'size': 30})cb.ax.tick_params(which='major', direction='in', length=6, labelsize=30)

這里使用了 get_my_cmaps(）函數(shù)，這是自己寫的，在文章開頭有說明，要使用他得到NCL的色標(biāo)MPL_hsv，同時(shí)需要在文件夾中放入MPL_hsv.rgb文件，有需要的朋友可以留下郵箱，作者給你們發(fā)。

（五）、850風(fēng)場(chǎng)箭頭繪制

#畫風(fēng)流場(chǎng)step = 5#設(shè)置步長，防止箭頭過密ax.quiver(lon[::step],lat[::step],u_850[::step,::step],v_850[::step,::step],color='w',scale=600, width=0.002, pivot='mid', transform=ccrs.PlateCarree())

（六）、其他設(shè)置，制作GIF

#顯示預(yù)報(bào)時(shí)次ax.text(0.85, -0.03, '21122420+' + TIME,transform=ax.transAxes, fontsize=50)#設(shè)置圖片標(biāo)題ax.set_title('ECMWF 850mb Temperature(shaded),Wind(vector)&500mb Geopotential Height(contour)',color='k',fontsize= 60)#保存單張圖片plt.savefig(TIME+'.png', dpi=30)#清理畫布plt.clf()#將繪制gif需要的靜態(tài)圖片名放入列表frames.append(imageio.imread(TIME+'.png')) #跳出循環(huán)，制作GIF并保存 imageio.mimsave('aa.gif', frames, 'GIF', duration=1.0)

三、完整代碼

from read_mdfs import MDFS_Grid import re import matplotlib.cm from matplotlib import colors import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import interpolateimport cartopy.crs as ccrs from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER import cartopy.feature as cfeature import os import imageioclass Colormap(colors.ListedColormap):def __init__(self, c, name='from_list', n=None):self.colors = cself.name = nameself.N = nsuper(Colormap, self).__init__(self.colors, name=self.name, N=self.N)def __getitem__(self, item):return Colormap(self.colors[item], name='sliced_' + self.name)def show(self):a = np.outer(np.ones(10), np.arange(0, 1, 0.001))plt.figure(figsize=(2.5, 0.5))plt.subplots_adjust(top=0.95, bottom=0.05, left=0.01, right=0.99)plt.subplot(111)plt.axis('off')plt.imshow(a, aspect='auto', cmap=self, origin='lower')plt.text(0.5, 0.5, self._name,verticalalignment='center', horizontalalignment='center',fontsize=12, transform=plt.gca().transAxes)plt.show() def coltbl(cmap_file):pattern = re.compile(r'(\d\.?\d*)\s+(\d\.?\d*)\s+(\d\.?\d*).*')with open(cmap_file) as cmap:cmap_buff = cmap.read()cmap_buff = re.compile('ncolors.*\n').sub('', cmap_buff)if re.search(r'\s*\d\.\d*', cmap_buff):return np.asarray(pattern.findall(cmap_buff), 'f4')else:return np.asarray(pattern.findall(cmap_buff), 'u1') / 255. def get_my_cmaps(cname):try:if cname in matplotlib.cm._cmap_registry:return matplotlib.cm.get_cmap(cname)except:passcmap_file = os.path.join( cname+ ".rgb")cmap = Colormap(coltbl(cmap_file), name=cname)matplotlib.cm.register_cmap(name=cname, cmap=cmap)return cmapdef get_lat_lon(filepath):#獲取經(jīng)緯度數(shù)據(jù)a = MDFS_Grid(filepath)lon = a.data['Lon']lat = a.data['Lat'][::-1]#翻了緯度lon_scipy = np.arange(lon.min(), lon.max(), 0.05)#插值lat_scipy = np.arange(lat.min(), lat.max(), 0.05)return lon,lat,lon_scipy,lat_scipy def read_data_hgt_or_tmp(filepath,lon,lat,lon_scipy,lat_scipy):#讀位勢(shì)高度和溫度a = MDFS_Grid(filepath)data = a.data['Grid'][::8,::8]#為了數(shù)據(jù)更平滑，降低了精度spline = interpolate.RectBivariateSpline(lat[::8], lon[::8],data,)data = spline(lat_scipy,lon_scipy)#插值return data[::-1,:]#翻了緯度 def read_data(filepath):#讀取一般數(shù)據(jù)，這里是風(fēng)a = MDFS_Grid(filepath)data = a.data['Grid']return data[::-1,:]#翻了緯度with open('CN-border-La.dat') as src:context = src.read()blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks] #畫布的設(shè)置，也是一次就夠了，畫完一張圖清空一次就好了，減少運(yùn)行內(nèi)存。 fig = plt.figure(figsize=(47, 30), dpi=30) #用來存放圖片的列表，為的是制作gif frames = [] #做循環(huán)用的時(shí)次 TIMES = ['000','003','006','009','012'] #循環(huán)開始 for TIME in TIMES:#讀取經(jīng)緯度、500位勢(shì)高度file_path = r"ECMWF_HR\ECMWF_HR_HGT_500_21122420." + TIMElon ,lat,lon_scipy,lat_scipy = get_lat_lon(file_path)hgt_500 = read_data_hgt_or_tmp(file_path,lon ,lat,lon_scipy,lat_scipy)#讀取850溫度file_path = r"ECMWF_HR\ECMWF_HR_TMP_850_21122420." + TIMEtem_850 = read_data_hgt_or_tmp(file_path,lon ,lat,lon_scipy,lat_scipy)#讀取850U風(fēng)file_path = r"ECMWF_HR\ECMWF_HR_UGRD_850_21122420." + TIMEu_850 = read_data(file_path)#讀取850V風(fēng)file_path = r"ECMWF_HR\ECMWF_HR_VGRD_850_21122420." + TIMEv_850 = read_data(file_path)#加子圖，做一些地圖設(shè)置ax = fig.add_axes([0.02, 0.05, 0.9, 0.9], projection=ccrs.LambertConformal(central_latitude=90, central_longitude=105))extent = [80,130,20,55]ax.set_extent(extent,crs=ccrs.Geodetic())gl = ax.gridlines( draw_labels=True, linewidth=2, color='k', alpha=0.5, linestyle='--')gl.xformatter = LONGITUDE_FORMATTER ##坐標(biāo)刻度轉(zhuǎn)換為經(jīng)緯度樣式gl.yformatter = LATITUDE_FORMATTERgl.xlabel_style = {'size': 30}gl.ylabel_style = {'size': 30}resolution_map = '50m'ax.add_feature(cfeature.OCEAN.with_scale(resolution_map))ax.add_feature(cfeature.LAND.with_scale(resolution_map))ax.add_feature(cfeature.RIVERS.with_scale(resolution_map))ax.add_feature(cfeature.LAKES.with_scale(resolution_map))for line in borders:ax.plot(line[0::2], line[1::2], '-', color='k',transform=ccrs.Geodetic())#畫500位勢(shì)高度ct = ax.contour(lon_scipy, lat_scipy, hgt_500, 20, colors='k',linewidths=3, transform=ccrs.PlateCarree())ax.clabel(ct, inline=True, fontsize=30, fmt='%d')#單畫一遍588，紅色clev = [588]ct = ax.contour(lon_scipy, lat_scipy, hgt_500, clev, colors='r',linewidths=3, transform=ccrs.PlateCarree())ax.clabel(ct, inline=True, fontsize=30, fmt='%d')#畫0℃線，白色clev = [0.0]ct = ax.contour(lon_scipy, lat_scipy, tem_850, clev, colors='w',linewidths=5, transform=ccrs.PlateCarree())ax.clabel(ct, inline=True, fontsize=30, fmt='%d')#畫850溫度clevs = range(-45,46,3)#引入NCL色標(biāo)my_cmap = get_my_cmaps("MPL_hsv")#翻轉(zhuǎn)色標(biāo)my_cmap = my_cmap.reversed()#填色cf = ax.contourf(lon_scipy, lat_scipy, tem_850, clevs, transform=ccrs.PlateCarree(), cmap=my_cmap)#設(shè)置色標(biāo)position = fig.add_axes([0.92, 0.05, 0.02, 0.9])cb = fig.colorbar(cf, cax=position, orientation='vertical',ticks=range(-45,46,3))cb.set_label('Temperature ℃', fontdict={'size': 30})cb.ax.tick_params(which='major', direction='in', length=6, labelsize=30)#畫風(fēng)流場(chǎng)step = 5#設(shè)置步長，防止箭頭過密ax.quiver(lon[::step],lat[::step],u_850[::step,::step],v_850[::step,::step],color='w',scale=600, width=0.002, pivot='mid', transform=ccrs.PlateCarree())#顯示預(yù)報(bào)時(shí)次ax.text(0.85, -0.03, '21122420+' + TIME,transform=ax.transAxes, fontsize=50)#設(shè)置圖片標(biāo)題ax.set_title('ECMWF 850mb Temperature(shaded),Wind(vector)&500mb Geopotential Height(contour)',color='k',fontsize= 60)#保存單張圖片plt.savefig(TIME+'.png', dpi=30)#清理畫布plt.clf()#將繪制gif需要的靜態(tài)圖片名放入列表frames.append(imageio.imread(TIME+'.png')) #跳出循環(huán)，制作GIF并保存 imageio.mimsave('aa.gif', frames, 'GIF', duration=1.0)

這樣，圖片就復(fù)刻成功了

需要數(shù)據(jù)文件和代碼的請(qǐng)留言郵箱，作者給發(fā)

總結(jié)

以上是生活随笔為你收集整理的Cartopy画地图第八天（冷空气南下，NCL色标使用）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。