本節(jié)提要：仿制中央氣象臺(tái)氣象服務(wù)圖片、關(guān)于cartopy里的投影與轉(zhuǎn)換、cartopy中extent與boundary。

一、仿制中央氣象臺(tái)圖片

從鄙人高三填報(bào)了南信的志愿開始，就一直持續(xù)的關(guān)注中央氣象臺(tái)，也算是一個(gè)老看客了。期間還經(jīng)歷了比較大型的改版(我感覺改版后變丑了O(∩_∩)O哈哈~)，我覺得央臺(tái)的氣象服務(wù)圖片做得比較好，基本上以填色圖和散點(diǎn)圖為主，輔以少量雷達(dá)圖，絕大多數(shù)圖片都能很好的向各層次的讀者展現(xiàn)天氣預(yù)報(bào)與實(shí)況。

今天就接著上次談?wù)摰膱D例和色條，來(lái)談?wù)刾y下仿制央臺(tái)圖片流程。

我仿制的圖片如下：

A、仿制的高溫圖

恩施想要這樣的高溫還是比較困難的(恩施算得上矮高原了，鄂西涼都)，所以修改成了地質(zhì)災(zāi)害的預(yù)報(bào)。

首先涉及到資料的問題，地質(zhì)災(zāi)害不在常規(guī)預(yù)報(bào)里，但是氣象局也必須發(fā)這項(xiàng)預(yù)警。這里只能用事先做好的實(shí)驗(yàn)數(shù)據(jù)，預(yù)報(bào)在清江兩岸有一定的山洪泥石流風(fēng)險(xiǎn)。使用的仍然是不均勻的站點(diǎn)資料，所以先要將站點(diǎn)資料網(wǎng)格化，變?yōu)楦顸c(diǎn)資料后再用等值線填色的方式畫出危險(xiǎn)區(qū)域。最后，通過前面幾節(jié)提到的添加圖例的方法，完善圖形。關(guān)鍵代碼如下：

olon=np.linspace(108,111,90)olat=np.linspace(29,32,90)olon,olat=np.meshgrid(olon,olat)func=Rbf(lon,lat,danger,function='linear')danger_new=func(olon,olat)colorlevel=[0,60,80,100]#危險(xiǎn)等級(jí)colordict=['white','orange','red']#顏色列表danger_map=mcolors.ListedColormap(colordict)#產(chǎn)生顏色映射norm=mcolors.BoundaryNorm(colorlevel,danger_map.N)#生成索引ax.contourf(olon,olat,danger_new,levels=colorlevel,cmap=danger_map,norm=norm)#填色圖clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000)#白化ax.set_title('地質(zhì)災(zāi)害風(fēng)險(xiǎn)(sample data)',fontsize=12)danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")low_danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="orange")labels = ['泥石流高風(fēng)險(xiǎn)','內(nèi)澇積水風(fēng)險(xiǎn)']ax.legend([danger, low_danger], labels,fontsize=7,frameon=False,title='圖例',loc='lower left', bbox_to_anchor=(-0.01, 0), fancybox=False)

B、降水量圖的仿制

使用的是最開始實(shí)驗(yàn)自定義colorbar時(shí)的那張圖的數(shù)據(jù)，但是當(dāng)時(shí)用的是色條來(lái)表示降水量，這次我們用圖例的方式表示降水量，前面的步驟和A中的類似。最后添加圖例時(shí)候，我用的是逐個(gè)添加，實(shí)際上可以用循環(huán)的方式添加。代碼如下：

filename=r'C:\Users\lenovo\Desktop\ex1.xlsx'df=pd.read_excel(filename)lon=df['lon']#經(jīng)度lat=df['lat']#緯度olon=np.linspace(108,111,90)olat=np.linspace(29,32,90)olon,olat=np.meshgrid(olon,olat)rain=df['rain_new']func=Rbf(lon,lat,rain,function='linear')rain_new=func(olon,olat)colorlevel=[0.1,10.0,25.0,50.0,100.0,250.0,500.0]#雨量等級(jí)colordict=['#A6F28F','#3DBA3D','#61BBFF','#0000FF','#FA00FA','#800040']#顏色列表rain_map=mcolors.ListedColormap(colordict)#產(chǎn)生顏色映射norm=mcolors.BoundaryNorm(colorlevel,rain_map.N)#生成索引cs= ax.contourf(olon,olat,rain_new,levels=colorlevel,cmap=rain_map,norm=norm)clip=maskout.shp2clip(cs, ax,r'E:\dijishi\cn_province.shp' ,420000)#白化############添加圖例#######################larger1?=?mpatches.Rectangle((0,?0),?1,?1,?facecolor="#A6F28F")larger2 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#3DBA3D")larger3 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#61BBFF")larger4 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#0000FF")larger5 = mpatches.Rectangle((0, 0), 1, 1, facecolor="#FA00FA")labels = ['0~10','10~25','25~50','50~100','100~250']ax.legend([larger1,larger2,larger3,larger4,larger5], labels,fontsize=12,frameon=False,title='圖例(mm)',loc='lower left', bbox_to_anchor=(-0.01, -0.03), fancybox=False)

其中，mpatches.Rectangle這個(gè)語(yǔ)句就是添加圖例上的色條，rectangle就是矩形的意思，延伸的，你也可以使用其他形狀添加圖例(包括自定義形狀)。如何提供更豐富的標(biāo)記，只能煩請(qǐng)各位在官網(wǎng)文檔上查找了。

二、Cartopy里的投影與轉(zhuǎn)換

我在兩個(gè)月前經(jīng)常碰到這個(gè)問題，有兩個(gè)關(guān)于投影的—crs、transform。前一個(gè)設(shè)定投影方式，后一個(gè)涉及投影與數(shù)據(jù)轉(zhuǎn)換。

首先說crs，這個(gè)是GeoAxes的基礎(chǔ)，只有在projection設(shè)置了投影之后，才能添加地圖。如果沒有一步轉(zhuǎn)化，有一些專屬于GeoAxes的命令就無(wú)法生效：

'AxesSubplot' object has no attribute 'add_feature'

比較典型的如add_feature、add_geometries添加地圖和地理特征的命令不能使用。

crs還控制著數(shù)據(jù)繪圖與邊界的裁剪，比如set_extent(crs=ccrs.PlateCarree())?，就使裁剪的方式按照PlateCarree()的方式進(jìn)行邊界的裁剪，一個(gè)經(jīng)典的案例即蘭勃脫下的使用extent進(jìn)行裁剪，那么將會(huì)在地圖中出現(xiàn)范圍之外的地圖：

如果主地圖投影為默認(rèn)投影方式(PlateCarree)，那么，你在后續(xù)的繪制中是不需要進(jìn)行投影轉(zhuǎn)換的。

但是，在主地圖不是默認(rèn)投影方式時(shí)，需要進(jìn)行轉(zhuǎn)換。比如：

ax=fig.add_subplot(1,1,1,projection=ccrs.LambertConformal(central_longitude=110,?central_latitude=35))ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r')

這種情況下，將無(wú)法繪制出數(shù)據(jù)：

這是因?yàn)?#xff0c;數(shù)據(jù)默認(rèn)的格式為PlateCarree下的投影，需要transform的介入：

ax.contourf(lon, lat, temp,levels=np.arange(-40,40),cmap='RdBu_r',transform=ccrs.PlateCarree())

transform=PlateCarree()的意思是，這個(gè)數(shù)據(jù)原來(lái)是PlateCarree()投影下的，轉(zhuǎn)化為主圖的格式后(即LambertConformal)，再繪制等值線。

主圖ax的投影方式

數(shù)據(jù)的投影方式

轉(zhuǎn)換命令

LambertConformal

PlateCarree

transform=PlateCarree

LambertConformal

Mercator

transform=Mercator

PlateCarree

Mercator

transform=Mercator

三、通過set_boundary和set_extent繪制扇形圖

cartopy官網(wǎng)實(shí)例有一個(gè)

Custom Boundary Shape，通過該例子可以裁剪邊框形狀。現(xiàn)在嘗試裁剪東亞地區(qū)的蘭勃脫投影。

import numpy as npfrom matplotlib.path import Pathimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeatvertices = [(65, 10), (65, 60), (145, 60), (145, 10), (65, 10)]#五個(gè)點(diǎn)，但是首尾是一樣的，以連接為封閉的四邊形boundary = Path(vertices)#邊界形狀fig=plt.figure(figsize=(5,5),dpi=500)ax=plt.axes([0,0,1,1],projection=ccrs.LambertConformal(central_longitude=110, central_latitude=35))ax.add_feature(cfeat.OCEAN)#添加海洋ax.add_feature(cfeat.LAND)#添加陸地ax.set_boundary(boundary, transform=ccrs.PlateCarree())#裁剪邊界ax.set_extent([65,145,0,60])#裁剪地圖范圍

這時(shí)，可能因?yàn)橥队稗D(zhuǎn)換的原因，出現(xiàn)了幾個(gè)明顯的折痕，于是，進(jìn)一步提高path上點(diǎn)的數(shù)量：

latmin = 10latmax = 60lonmin = 70lonmax = 140lats = np.linspace(latmax, latmin, latmax - latmin + 1)lons = np.linspace(lonmin, lonmax, lonmax - lonmin + 1)vertices = [(lon, latmin) for lon in range(lonmin, lonmax + 1, 1)] + \[(lon, latmax) for lon in range(lonmax, lonmin - 1, -1)]boundary = mpath.Path(vertices)ax.set_boundary(boundary, transform=ccrs.PlateCarree())ax.set_extent(extent)

此時(shí)，因?yàn)関ertices中的邊界點(diǎn)的數(shù)量增多，邊界連接更加圓滑。請(qǐng)注意，在這種方式下裁剪了邊界，一定要使用set_extent來(lái)裁剪地圖大小，不然就會(huì)出現(xiàn)地圖蜷縮在圖里這種情況：

在cartopy=0.17中，不能使用draw_labels=True來(lái)為除PlateCarree、Mercator之外的投影添加經(jīng)緯標(biāo)簽，不過據(jù)說在0.18版本中已經(jīng)優(yōu)化，讀者可以試試。通過前面提到的數(shù)據(jù)轉(zhuǎn)換，繪制了一張氣溫圖：

ax.contourf(data_lon,data_lat,temp,levels=np.arange(-40,40),cmap='Spectral_r',transform=ccrs.PlateCarree())

往期回顧：

歡迎關(guān)注云臺(tái)書使獲取更多

總結(jié)

以上是生活随笔為你收集整理的怎么用python画天气图_Python气象绘图教程（十五）—Cartopy_5的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

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