創作內容

5 GP

(筆記) python的cartopy使用、清除已畫的資料方法

作者：已經改掉暱稱的米奇│2019-03-09 02:53:08│巴幣：10│人氣：3444

記錄一下

之前畫地理地圖相關的東西時都是用basemap

但是basemap已經停止更新了，維護工作也只會持續到2020

而官方之後的心力都會放在cartopy上

雖然就算停止維護了還是可以正常使用basemap啦

但是保險起見，再加上既然官方都那樣建議了

看來還是先乖乖學一下cartopy比較好

basemap的流程大概是這樣：

利用Basemap(...)建立basemap物件

然後利用這個物件的各種method作圖

而cortopy則不一樣

將cortopy import進來後

只需要在平常使用的matplotlib函式中加個'projection'轉換成cartopy的格式

之後就可以跟平常一樣的方式畫圖

(有點類似用matploblib畫3D圖那樣)

所以原本用basemap的話會像這樣：

m = Basemap(....)

m.contourf(...)

m.quiver(....)

....

cartopy則會是這樣

fig, ax = plt.subplots(subplot_kw={'projection': .....})

ax.conourf(.....)

ax.quiver(......)

就不需要再透過其他物件來畫了

---

拿以前的例子來試試

這三張沒什麼特別的

就跟之前那篇basemap的差不多

只是import變成cartopy

這張就是開始用cartopy去畫了

可以看到在subplots()裡面，利用subplot_kw，傳入一個dictionary

其中包含了{'projection' : ccrs.PlateCarree()}

而ccrs.PlateCarree當然就是要畫的地圖投影方式了

另外我這邊是用subplots來做

不過只要任何可以建立axes、而且包含project選項的函數都可以

像是這樣也是可以的：

ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

轉換成cartopy的座標後

再利用ax.set_extent將範圍限制在我們想要的範圍內

也就是cartopy會先建立一個全世界的地圖，再依照想要的範圍去限制

不像basemap是一開始就會需要給定範圍才能建立地圖

不過有趣的是雖然cartopy一開始會先建出全世界，但是執行速度似乎不會因此比較慢

剩下就蠻直觀的

ax.coastlines()就是畫出海岸線，設定解析度是用'110m'、'50m'、'10m'這樣

ax.gridlines()是畫出經緯度線，draw_labels就是設定要不要順便畫出xy的tickslabel

不過畫出來會整個圖四個邊都有label

我自己希望只要左邊和下面有就好，所以特別設為False

(不過預設也是False就是了)

換個例子看看，這次是溫度

再看一個例子，這次是風場

用箭頭(風向風速)疊contourf(風速)

--------------------------------------------------------------------------------

接下來是看要怎麼樣清除圖上已經畫的資料

因為常常會遇到一種情況，就是需要畫非常多不同但是類似的圖

比方說同樣畫降雨，但是要把0~60個時間的降雨畫出來，所以共61張圖這樣

那最簡單的作法就是直接寫個迴圈

每一張圖就重畫一次

但是這樣其實蠻浪費時間的

畢竟每跑一次就要重新畫地圖、海岸線、格線之類的

偏偏這些東西其實每張圖都是一樣的，但是卻要一直重畫

所以如果想要把已經畫的資料清掉

可以利用ax.collections.clear()去清掉

ax.collections是一個list，包含了當下圖上已畫出來的東西

像是這樣

以這個前面畫風速的為例子

這個colletions包含了10個已經畫出來的東西

前7個都是contourf相關的

第8個是quiver

後面兩個則是x和y方向的格線

所以如果利用ax.collections.clear()

這整個list會被清空

所以像contourf、quiver之類的也會從圖上消失

但是像地圖、海岸線等等都還會留著

像是這樣

清掉之後就像這樣

原本的contourf都不見了

所以如果要畫各個時間的圖

就不需要一直全部重畫了

只需要最開始先畫完地圖、海岸線等等

之後每畫完一次資料就清掉，一直重複到結束就可以了

最後比較一下

如果用原本的方式，畫完一張就關掉再全部重畫一次 (作法1)

或是用ax.collections.clear()的方式清掉資料 (作法2)

這兩種方法的速度會差多少

可以看到一直重畫的做法1，需要31秒來畫61張圖

相較之下，利用ax.collections.clear()來處理，則只需要18.9秒

這個差距在圖的內容更多、圖片張數更多，或是更精細時會差更多

-------

整體來說
cartopy的使用上有比basemap更簡潔一點

而且執行效率也明顯比較高

不過cartopy看起來還是有缺一些東西

但是之後basemap的開發改成去做cartopy

應該還是會越來越完整的

喜歡 5 收藏 2 引用留言推上首頁檢舉

引用網址：https://home.gamer.com.tw/TrackBack.php?sn=4318912
All rights reserved. 版權所有，保留一切權利

相關創作

北藝動畫準備大秘辛-與教授的面對面訪談

【漫畫】シロサワ-別班女生是忠於慾望的傻瓜(2/6)

NVIDIA GeForce Game Ready Driver 552.44 WHQL

友人糯米贈2024生日賀圖、魔術烏骨雞

【七曜】

留言共 9 篇留言

青蛙-ちずる婆爆：
霸...你的程式能力根本屌打一堆人包跨我嗚嗚電機系不會寫程式是不是要ㄘ屎ㄌ

03-18 00:12

已經改掉暱稱的米奇：
沒ㄅ這很多人都會我只是整理了一下... 而且今天去企博還被人資鄙視唉03-18 00:24

Alan：
你好我想請問一下如果我想判斷降雨天數畫出網格圖要如何選取每一個網格點，我圖是用pcolor畫的，謝謝。

08-19 00:30

已經改掉暱稱的米奇：
有點不太懂你的意思耶有詳細一點的說明嗎
08-19 12:15

Alan：
抱歉我講的不是很清楚，我講詳細一點
你這篇的降雨量圖是用contourf畫的，我用pcolor把降雨量的空間分布變成網格點，呈現出來的樣子就會變成像馬賽克一樣每一格都有對應不同降雨量的顏色。
那我想要知道每一格的降雨量並去判斷是否有降雨，但我現在是卡在無法知道每一格的降雨量，所以想問問看有沒有辦法解決，謝謝。

08-19 14:55

已經改掉暱稱的米奇：
不過你要畫pcolor的話，會自己先有對應的raw data不是嗎？直接從raw data判斷每格有沒有降雨就可以了ㄅ？還是有簡單的code可以看一下你是怎麼畫圖的嗎，不然我還是不太懂卡住的點會在哪[e3]08-19 15:51

Alan：
這邊好像無法放圖片，我簡單打一下好了，我使用的資料是5km解析的日降雨資料
rainday = [] #回傳位置
for itime in range(0, 365): #時間為一年365天
for ilon in range(0, 60): #經度資料有60個
for ilat in range(0, 81): #緯度資料有81個
m = data_rain[itime, ilon, ilat] #data_rain是我的降雨資料
time = 0
if m > 0.1: #判斷是否降雨
time = time + 1
rainday.append(time) #回傳是否有降雨(1:降雨，0:無降雨)
我想知道空間降雨日數然後像降雨圖一樣畫出來，但前面光是判斷就先卡住了，不知道這樣寫對不對，最近一個月才學使用python，所以寫法有問題還請見諒

08-19 16:27

已經改掉暱稱的米奇：
你的說明跟code好像對不太起來，不過你要算的應該是每個經緯度座標上，各自有多少天降雨量超過門檻，最後把各經緯度的降雨天數當作畫圖的顏色對ㄅ？是的話我大概寫了一下你參考看看 https://ideone.com/RFiksB
08-20 02:22

已經改掉暱稱的米奇：
(前面都只是我隨便建資料，看22行以後就好)
然後我假設你有裝numpy而且你的data_rain是numpy array
08-20 02:24

Alan：
成功了，我一直以為要用迴圈下去寫，原來不用這麼麻煩，真的非常感謝幫忙!!!

08-20 17:25

已經改掉暱稱的米奇：
不會 python裡能不用迴圈就別用尤其是針對numpy array更是如此08-20 20:43

Alan：
不好意思又再麻煩你，想請問一下如果想要加上縣市邊界要怎麼寫?謝謝。

08-29 19:40

已經改掉暱稱的米奇：
可以查一下shapefile，下載帶有台灣縣市邊界的shapefile，然後cartopy可以讀取並畫出來，可以先試看看08-30 12:52

已經改掉暱稱的米奇：
我都還沒測過，不過shapefile可以試試這份 https://data.gov.tw/dataset/7442，選SHP那個下載，cartopy怎麼畫可以參考這篇 https://techoverflow.net/2021/04/25/how-to-plot-shapefile-data-in-cartopy/08-30 13:05

已經改掉暱稱的米奇：
然後我印象中coastline跟縣市邊界會不重疊，所以如果要畫縣市邊界的話可以考慮不要畫coastline08-30 13:06

Alan：
好的，我試試看，感謝你。

08-31 14:00

彥：
你好，我在學習metpy 時無意間找到您的文章
並試著講您的code照著打
但是我在執行到
cn = ax.contourf(lon[0,:,:], lat[0,:,:], rain[itime,:,:],
cmap=precip_colormap, alpha=0.65,
norm=norm, levels=cLevel, linestyles=None)
這一段時出現以下錯誤，我找不出原因，請問是發生什麼問題了?
IndexError: index 18 is out of bounds for axis 0 with size 1

11-10 14:42

已經改掉暱稱的米奇：
這個IndexError就是指有某個陣列(假設叫A)，A的維度0的大小只有1，但它卻想要取A維度0的第18個位子，所以就跳這個錯誤
我猜可能是前面準備各變數的時候有問題，這可能要看一下你裡面各個參數的shape才知道
你可以在執行這行之前把各變數的大小確定一下
我猜應該是cLevel大小不對，或是準備precip_colormap跟norm的地方有錯，要檢查一下11-11 12:53

彥：
感謝，解答
回家試試看

11-11 12:55

已經改掉暱稱的米奇：
OK，如果有看到問題再跟我說，搞不好是我這邊有bug哈11-11 12:59

我要留言提醒：您尚未登入，請先登入再留言

5 喜歡★micky85lu 可決定是否刪除您的留言，請勿發表違反站規文字。

前一篇：一維波動方程模擬... 後一篇：(筆記) python ...

追蹤私訊切換新版閱覽

作品資料夾

未分類 (164)

happy545 給晚上好~：
希望大家會喜歡我寫的故事，歡迎多加分享喔~~XD看更多我要大聲說昨天19:54