本節提要：仿制中央氣象臺氣象服務圖片、關于cartopy里的投影與轉換、cartopy中extent與boundary。

一、仿制中央氣象臺圖片

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

今天就接著上次談論的圖例和色條，來談談py下仿制央臺圖片流程。

我仿制的圖片如下：

A、仿制的高溫圖

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

首先涉及到資料的問題，地質災害不在常規預報里，但是氣象局也必須發這項預警。這里只能用事先做好的實驗數據，預報在清江兩岸有一定的山洪泥石流風險。使用的仍然是不均勻的站點資料，所以先要將站點資料網格化，變為格點資料后再用等值線填色的方式畫出危險區域。最后，通過前面幾節提到的添加圖例的方法，完善圖形。關鍵代碼如下：

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]#危險等級colordict=['white','orange','red']#顏色列表danger_map=mcolors.ListedColormap(colordict)#產生顏色映射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('地質災害風險(sample data)',fontsize=12)danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="red")low_danger = mpatches.Rectangle((0, 0), 1, 1, facecolor="orange")labels = ['泥石流高風險','內澇積水風險']ax.legend([danger, low_danger], labels,fontsize=7,frameon=False,title='圖例',loc='lower left', bbox_to_anchor=(-0.01, 0), fancybox=False)

B、降水量圖的仿制

使用的是最開始實驗自定義colorbar時的那張圖的數據，但是當時用的是色條來表示降水量，這次我們用圖例的方式表示降水量，前面的步驟和A中的類似。最后添加圖例時候，我用的是逐個添加，實際上可以用循環的方式添加。代碼如下：

filename=r'C:\Users\lenovo\Desktop\ex1.xlsx'df=pd.read_excel(filename)lon=df['lon']#經度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]#雨量等級colordict=['#A6F28F','#3DBA3D','#61BBFF','#0000FF','#FA00FA','#800040']#顏色列表rain_map=mcolors.ListedColormap(colordict)#產生顏色映射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這個語句就是添加圖例上的色條，rectangle就是矩形的意思，延伸的，你也可以使用其他形狀添加圖例(包括自定義形狀)。如何提供更豐富的標記，只能煩請各位在官網文檔上查找了。

二、Cartopy里的投影與轉換

我在兩個月前經常碰到這個問題，有兩個關于投影的—crs、transform。前一個設定投影方式，后一個涉及投影與數據轉換。

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

'AxesSubplot' object has no attribute 'add_feature'

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

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

如果主地圖投影為默認投影方式(PlateCarree)，那么，你在后續的繪制中是不需要進行投影轉換的。

但是，在主地圖不是默認投影方式時，需要進行轉換。比如：

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')

這種情況下，將無法繪制出數據：

這是因為，數據默認的格式為PlateCarree下的投影，需要transform的介入：

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

transform=PlateCarree()的意思是，這個數據原來是PlateCarree()投影下的，轉化為主圖的格式后(即LambertConformal)，再繪制等值線。

主圖ax的投影方式

數據的投影方式

轉換命令

LambertConformal

PlateCarree

transform=PlateCarree

LambertConformal

Mercator

transform=Mercator

PlateCarree

Mercator

transform=Mercator

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

cartopy官網實例有一個

Custom Boundary Shape，通過該例子可以裁剪邊框形狀?，F在嘗試裁剪東亞地區的蘭勃脫投影。

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)]#五個點，但是首尾是一樣的，以連接為封閉的四邊形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])#裁剪地圖范圍

這時，可能因為投影轉換的原因，出現了幾個明顯的折痕，于是，進一步提高path上點的數量：

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)

此時，因為vertices中的邊界點的數量增多，邊界連接更加圓滑。請注意，在這種方式下裁剪了邊界，一定要使用set_extent來裁剪地圖大小，不然就會出現地圖蜷縮在圖里這種情況：

在cartopy=0.17中，不能使用draw_labels=True來為除PlateCarree、Mercator之外的投影添加經緯標簽，不過據說在0.18版本中已經優化，讀者可以試試。通過前面提到的數據轉換，繪制了一張氣溫圖：

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

往期回顧：

歡迎關注云臺書使獲取更多

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。