cartopy繪圖指南

2022-07-20 19:22 作者:八八不愛叭叭 0人讀過 | 我要投稿

前言

嗨，你好，我是來自點點GIS的南南

我與Cartopy的認識起源于“氣象水文科研貓”的這個推文，那時候的我覺得，用代碼畫地圖好酷，arcgis就感覺low爆了。但是一直沒有時間學(xué)習(xí)。前段時間放暑假，磕磕絆絆裝完包以后就不想動彈了，折騰環(huán)境折騰的我頭皮發(fā)麻。gdal和cartopy真的是我學(xué)python以來裝的最麻煩的包。

很多時候我想學(xué)習(xí)，卻感覺無從下手，大佬們的博文一大堆，我看了以后只有一句“臥槽，牛逼”，然后沒然后了。前幾天某人欺負我，把我整抑郁了，大晚上睡不著覺，就想著填個坑，雖然我技術(shù)菜，不能“炫技”，但寫個入門教程難不倒我吧。然后就有了這篇文章。

這篇文章參考了諸多大佬的博文，如氣象學(xué)家，云臺書使，氣象學(xué)人，好奇心Log，等等公眾號大佬。我作為初學(xué)者，寫的文肯定錯誤不少，歡迎大家給我留言指正，要是有老板叫我去打工就更好了，嘿嘿嘿

Emil：2531788189@qq.com

Cartopy介紹

Cartopy 是一個開源免費的第三方 Python 擴展包，由英國氣象辦公室的科學(xué)家們開發(fā)，支持 Python 2.7 和 Python 3，致力于使用最簡單直觀的方式生成地圖，并提供對 matplotlib 友好的協(xié)作接口。

在常用的python繪圖庫中，basemap，geopandas，pyecharts等，其中basemap已經(jīng)停止維護了，在前文中我已經(jīng)講到過，pyecharts是用于數(shù)據(jù)可視化等專業(yè)圖表的繪制，之前我的文章也介紹過，pyecharts在地學(xué)可視化中功能實在過于簡陋;geopandas是基于pandas的，一般用于商業(yè)數(shù)據(jù)分析，所以我更推薦大家學(xué)習(xí)Cartopy，雖然Cartopy氣象專業(yè)用得很多哈哈哈

Cartopy安裝

https://mp.weixin.qq.com/s/GW6odDBGLPVRZTKx0YZLVg

Cartopy繪圖入門

在Cartopy的使用過程中，我們常常需要搭配其他庫一起使用，常用的有如Matplotlib ，numpy，pandas，gma等

Cartopy投影

Cartopy 是利用 Matplotlib 來畫圖的，因此首先要導(dǎo)入 pyplot 模塊。在 Cartopy 中，每種投影都是一個類，被存放在 cartopy.crs 模塊中，crs 即坐標參考系統(tǒng)（Coordinate Reference Systems）之意。所以接著要導(dǎo)入這個模塊。

投影用法示例：

proj?=?ccrs.PlateCarree()???#proj?=?ccrs.投影()

這里簡單說明一下常用的三種投影，Cartopy 大概有三四十種投影，如國想詳細了解可以參考下面的文章

https://scitools.org.uk/cartopy/docs/v0.15/crs/projections.html默認投影PlateCarree墨卡托投影Mercator蘭勃脫投影Lambert

默認投影適合單獨省份或者地級市的繪制，蘭勃脫投影適合中緯度大范圍繪制，比如繪制全中國大公雞、東亞形勢、西北太平洋等；墨卡托投影適合低緯度赤道附近的繪制，一般研究臺風(fēng)、緯向環(huán)流等

繪制地圖示例

將 cartopy 集成到 matplotlib 的主要類是 GeoAxes，它是普通 matplotlib Axes 的子類。GeoAxes 類為特定于繪制地圖的軸添加了額外的功能。創(chuàng)建一個 GeoAxes 對象的辦法是，在創(chuàng)建 axes（或 subplot）時，通過參數(shù) projection 指定一個 ccrs 中的投影。這里便利用這一方法生成了一個等距圓柱投影下的 ax。

最后調(diào)用 ax 的方法 coastlines 畫出海岸線，默認以本初子午線為中心，比例尺為 1:110m（m 表示 million）。

因此用 Cartopy 畫地圖的基本流程并不復(fù)雜：

創(chuàng)建畫布。
通過指定 projection 參數(shù)，創(chuàng)建 GeoAxes 對象。
調(diào)用 GeoAxes 的方法畫圖。
GeoAxes 用法擴展（部分常用）

set_global：讓地圖的顯示范圍擴展至投影的最大范圍。例如，對 PlateCarree 投影的 ax 使用后，地圖會變成全球的。
set_extent：給出元組 (x0, x1, y0, y1) 以限制地圖的顯示范圍。
set_xticks：設(shè)置 x 軸的刻度。
set_yticks：設(shè)置 y 軸的刻度。
gridlines：給地圖添加網(wǎng)格線。
coastlines：在地圖上繪制海岸線。
stock_img：給地圖添加低分辨率的地形圖背景。
add_feature：給地圖添加特征（例如陸地或海洋的填充、河流等）。

import?matplotlib.pyplot?as?plt###引入庫包 import?cartopy.crs?as?ccrs proj?=?ccrs.PlateCarree() fig?=?plt.figure(figsize=(4,?4),?dpi=200)??#?創(chuàng)建畫布 ax?=?plt.axes(projection=ccrs.PlateCarree())#?創(chuàng)建子圖 ax.coastlines()#?調(diào)用ax的方法畫海岸線 ''' fig語法： figure(num=None,?figsize=None,?dpi=None,?facecolor=None,?edgecolor=None,?frameon=True) num:圖像編號或名稱，數(shù)字為編號?，字符串為名稱 figsize:指定figure的寬和高，單位為英寸； dpi參數(shù)指定繪圖對象的分辨率，即每英寸多少個像素，缺省值為80??????1英寸等于2.5cm,A4紙是?21*30cm的紙張? facecolor:背景顏色 edgecolor:邊框顏色 frameon:是否顯示邊框 '''

如需修改中央經(jīng)線，可以在創(chuàng)建子圖的時候指定一些參數(shù)

ax?=?plt.axes(projection=ccrs.PlateCarree(central_longitude=180))#全球圖像的中央位于太平洋的?180?度經(jīng)線處

也可以在繪制線條的時候指定線寬

ax.coastlines(lw=0.3)#線條寬度0.3

在地圖上添加特征

上面的地圖是一個十分粗糙的地圖，僅有海岸線，嘗試一下添加更多地理信息

添加預(yù)定義要素

首先需要導(dǎo)入一個cartopy.feature 常量，為了簡化一些非常常見的情況，如大陸海洋國界等cartopy都已經(jīng)進行了預(yù)定義，但需要注意的是直接導(dǎo)入的中國國界線并不是標準的，在使用時需要注意

這些信息帶有三種分辨率，分別為110m，50m和10m。

import?matplotlib.pyplot?as?plt###引入庫包 import?cartopy.crs?as?ccrs import?cartopy.feature?as?cfeature#預(yù)定義常量 proj?=?ccrs.PlateCarree() fig?=?plt.figure(figsize=(4,?4),?dpi=200)??#?創(chuàng)建畫布 ax?=?plt.axes(projection=ccrs.PlateCarree(central_longitude=180))#全球圖像的中央位于太平洋的?180?度經(jīng)線處 ax.add_feature(cfeature.LAND)#添加陸地 ax.add_feature(cfeature.COASTLINE,lw?=?0.3)#添加海岸線 ax.add_feature(cfeature.RIVERS,lw?=?0.25)#添加河流 #ax.add_feature(cfeat.RIVERS.with_scale('50m'),lw?=?0.25)??#?加載分辨率為50的河流 ax.add_feature(cfeature.LAKES)#添加湖泊 ax.add_feature(cfeature.BORDERS,?linestyle?=?'-',lw?=?0.25)#不推薦，因為該默認參數(shù)會使得我國部分領(lǐng)土丟失 ax.add_feature(cfeature.OCEAN)#添加海洋 ax.coastlines(lw=0.3)

在cartopy中，所有的線條都可以改變他的粗細。只用改變其參數(shù)lw = *即可

ax.add_feature(cfeature.RIVERS,lw?=?2)#改變河流粗細

同樣的，所有的線，面也可以改變它的顏色

ax.add_feature(cfeature.LAKES,color='r')#指定湖泊顏色為紅色

添加經(jīng)緯度標簽

這里還是要用到前文說過的GeoAxes 用法

set_xticks：設(shè)置 x 軸的刻度。
set_yticks：設(shè)置 y 軸的刻度。

例：x軸為(-180,-120,-60,0,60,120,180),y軸為(-90,-60, -30, 0, 30, 60,90)

x_extent?=?[?0,?60,?120,?180,?240,?300,?360]#經(jīng)緯度范圍,#直接使用(-180,-120,-60,0,60,120,180)會異常，需要寫成(0,?60,?120,?180,?240,?300,?360)的形式 y_extent?=?[?-90,-60,?-30,?0,?30,?60,90] fig?=?plt.figure(figsize=(4,?4),?dpi=200)?? ax?=?plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) ax.add_feature(cfeature.LAND)#添加陸地 ax.add_feature(cfeature.COASTLINE,lw?=?0.3)#添加海岸線 ax.add_feature(cfeature.RIVERS,lw?=?0.25)#添加河流 ax.add_feature(cfeature.LAKES)#指定湖泊顏色為紅色#添加湖泊 #ax.add_feature(cfeature.BORDERS,?linestyle?=?'-',lw?=?0.25)#不推薦，因為該默認參數(shù)會使得我國部分領(lǐng)土丟失 ax.add_feature(cfeature.OCEAN)#添加海洋 ax.set_xticks(x_extent,?crs=ccrs.PlateCarree())#添加經(jīng)緯度 ax.set_yticks(y_extent,?crs=ccrs.PlateCarree())

將經(jīng)緯度標簽轉(zhuǎn)換為具有單位的形式

from?cartopy.mpl.ticker?import?LongitudeFormatter,?LatitudeFormatter?#導(dǎo)入Cartopy專門提供的經(jīng)緯度的Formatter x_extent?=?[?0,?60,?120,?180,?240,?300,?360]#經(jīng)緯度范圍,#直接使用(-180,-120,-60,0,60,120,180)會異常，需要寫成(0,?60,?120,?180,?240,?300,?360)的形式 y_extent?=?[?-90,-60,?-30,?0,?30,?60,90] fig?=?plt.figure(figsize=(4,?4),?dpi=200)?? ax?=?plt.axes(projection=ccrs.PlateCarree(central_longitude=180)) ax.add_feature(cfeature.LAND)#添加陸地 ax.add_feature(cfeature.COASTLINE,lw?=?0.3)#添加海岸線 ax.add_feature(cfeature.RIVERS,lw?=?0.25)#添加河流 ax.add_feature(cfeature.LAKES)#指定湖泊顏色為紅色#添加湖泊 #ax.add_feature(cfeature.BORDERS,?linestyle?=?'-',lw?=?0.25)#不推薦，因為該默認參數(shù)會使得我國部分領(lǐng)土丟失 ax.add_feature(cfeature.OCEAN)#添加海洋 ax.set_xticks(x_extent,?crs=ccrs.PlateCarree())#添加經(jīng)緯度 ax.set_yticks(y_extent,?crs=ccrs.PlateCarree()) #?利用Formatter格式化刻度標簽 lon_formatter?=?LongitudeFormatter(zero_direction_label=False) lat_formatter?=?LatitudeFormatter() ax.xaxis.set_major_formatter(lon_formatter) ax.yaxis.set_major_formatter(lat_formatter)

經(jīng)緯網(wǎng)

在前文中GeoAxes 用法提到過添加經(jīng)緯網(wǎng)的函數(shù)gridlines，只用在上述代碼末尾中添加ax.gridlines()，即

ax.gridlines()#添加網(wǎng)格

這實在是太丑了，給它換個線樣式

ax.gridlines(linestyle='--')#添加網(wǎng)格，線樣式為'--'

配上網(wǎng)格總感覺刻度朝外怪怪的

ax.tick_params(color?=?'blue',direction='in')#更改刻度指向為朝內(nèi)，顏色設(shè)置為藍色

更多GeoAxes 用法可以參考以下文章

https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html

局部地圖

很多時候我們繪圖不會用到世界地圖，所以可以嘗試更改你的范圍。cartopy提供更改的方式是通過定義經(jīng)緯度范圍來進行更改的。在上述代碼末尾添加下方樣例代碼試試

set_extent：給出元組 (x0, x1, y0, y1) 以限制地圖的顯示范圍。

ax.set_extent([0,180,0,35],crs?=?ccrs.PlateCarree())?#選取經(jīng)度為0°E-180°E，緯度為0°N-90°N的區(qū)域

添加標題

ax.set_title('Cartopy')?#添加標題Cartopy

Cartopy繪圖進階

在前文中提到過，Cartopy的中國地圖邊界是有問題的，那么在日常使用中，我們該如何避免這些問題呢?

Cartopy中國地圖繪制方法

方法一

用 Cartopy 的 shapereader 讀取后即可繪制。 Cartopy 的 shapereader用法示例可參考以下文章

https://scitools.org.uk/cartopy/docs/v0.15/tutorials/using_the_shapereader.html

關(guān)于正確的行政邊界獲取可以參考以下博文

https://mp.weixin.qq.com/s/aUCjTRrV6Cz-k7mtXOEiGw

使用cartopy讀取shapefile繪制共有兩種方法,分別是add_feature和add_geometries，但無論你使用哪種方法，你都需要先讀取文件才能夠添加。Cartopy提供了一個基于pyshp的接口以實現(xiàn)對地理文件的簡單讀取和操作：

首先導(dǎo)入相關(guān)包

import?numpy?as?np import?cartopy.crs?as?ccrs import?matplotlib.pyplot?as?plt import?cartopy.feature?as?cfeat import?cartopy.io.shapereader?as?shapereader

你找他總得給他來一個地址吧

shp_path?=?"D:/data2/china/china/province_9south.shp"#切記路徑一定要英文，血的教訓(xùn)?。?！

繪圖三部曲，畫布，投影，子圖

fig??=?plt.figure(figsize=(6,6))?#創(chuàng)建畫布 proj?=?ccrs.PlateCarree()#創(chuàng)建投影，選擇cartopy的默認投影 ax???=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})??#子圖 extent?=?[70,140,2,60]#限定繪圖范圍

add_feature加載自定義地圖

china_map?=?cfeat.ShapelyFeature(shapereader.Reader(china_map).geometries(),?proj,?edgecolor='k',?facecolor='none') ax.add_feature(china_map,?linewidth=1) ax.set_extent(extent,?crs=proj)

add_geometries加載自定義地圖

ax.add_geometries(shapereader.Reader(china_map).geometries(),?crs=proj,?facecolor='none',edgecolor='k',linewidth=1) ax.set_extent(extent,?crs=proj)

完整代碼如下

import?numpy?as?np import?cartopy.crs?as?ccrs import?matplotlib.pyplot?as?plt import?cartopy.feature?as?cfeat import?cartopy.io.shapereader?as?shapereader shp_path?=?"D:/data2/china/china/province_9south.shp"#切記路徑一定要英文，血的教訓(xùn)?。?！ fig??=?plt.figure(figsize=(6,6))?#創(chuàng)建畫布 proj?=?ccrs.PlateCarree()#創(chuàng)建投影，選擇cartopy的默認投影 ax???=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})??#子圖 extent?=?[70,140,2,60]#限定繪圖范圍 china_map?=?cfeat.ShapelyFeature(shapereader.Reader(shp_path).geometries(),?proj,?edgecolor='k',?facecolor='none') ax.add_feature(china_map,?linewidth=1)#添加讀取的數(shù)據(jù)，設(shè)置線寬 ax.set_extent(extent,?crs=proj) #ax.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,?#facecolor='none',edgecolor='k',linewidth=1) #ax.set_extent(extent,?crs=proj)

方法二

使用GMT 中文社區(qū)上提供的省界的經(jīng)緯度數(shù)據(jù)進行繪圖。GMT 中文社區(qū)為用戶提供了一套名為CN-border的基本準確的數(shù)字化中國國界，、省界數(shù)據(jù)。在此向 GMT 中文社區(qū)的維護和貢獻者表示感謝！

但需要注意的是，即便正確繪制了國界、省界，所繪地圖如果要在境內(nèi)公開展示，依然需要審圖。個人沒有提交審圖申請的資格，需要付費，委托給有資質(zhì)的企事業(yè)單位代為提交審圖申請。下為GMT 中文社區(qū)網(wǎng)址，大家可自行下載

https://docs.gmt-china.org/latest/dataset-CN/CN-border/

CN-border 數(shù)據(jù)提供了三個數(shù)據(jù)文件：

CN-border-La.gmt：中國國界、省界、十段線以及南海諸島數(shù)據(jù)
CN-border-L1.gmt：中國國界、十段線以及南海諸島數(shù)據(jù)，不含省界數(shù)據(jù)
ten-dash-line.gmt：僅十段線數(shù)據(jù)

import?numpy?as?np import?matplotlib.pyplot?as?plt import?cartopy.crs?as?ccrs import?cartopy.feature?as?cfeature #加載邊界數(shù)據(jù)，CN-border-La.dat下載 #https://github.com/gmt-china/china-geospatial-data/releases with?open('C:/Users/啊啊啊啊/Desktop/china-geospatial-data-GB2312/CN-border-La.gmt')?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] #?畫布 fig?=?plt.figure(figsize=[10,?8]) #?設(shè)置投影并繪制主圖形 ax?=?plt.axes(projection=ccrs.LambertConformal(central_latitude=90, ???????????????????????????????????????????????central_longitude=105)) #?繪制線條 for?line?in?borders: ????ax.plot(line[0::2],?line[1::2],?'-',?lw=1,?color='k', ????????????transform=ccrs.Geodetic()) #?繪制邊界線 for?line?in?borders: ????sub_ax.plot(line[0::2],?line[1::2],?'-',?lw=1,?color='k', ????????????????transform=ccrs.Geodetic()) #?設(shè)置圖形范圍 sub_ax.set_extent([105,?125,?0,?25]) #?顯示圖形 plt.show()

方法三

cnmaps繪圖，cnmaps是Cartopy的一個擴展包，內(nèi)置了中國地圖數(shù)據(jù)源，該數(shù)據(jù)源自高德地圖API，大家可以自行去該項目進行學(xué)習(xí)

https://github.com/Clarmy/cnmaps

重點區(qū)域突出顯示

該示例是在前文方法一的基礎(chǔ)上進行繪制

Cartopy 默認對所有地物使用相同的外觀，如果需要突出顯示某些地物，就必須進行篩選。這里從province_9south.shp這份全國省級行政區(qū)劃數(shù)據(jù)中選中山西（通過屬性表SHENG_ID字段），然后使用 ax.add_geometries() 方法將它們加入到原有地圖元素中。

for?state,?geos?in?zip(shapereader.Reader(shp_path).records(),?shapereader.Reader(shp_path).geometries())?: ????if?state.attributes['SHENG_ID']?in?[14?]:#如果'14'在'SHENG_ID'這個字段中 ????????ax.add_geometries([geos],?crs=proj,facecolor='None',?edgecolor='blue')??#?重點區(qū)域上色

關(guān)于屬性表，你可以通過gdal來打印出該shpfile文件所具有的屬性表，也可以通過arcgis，qgis等GIS軟件打開查詢

僅繪制山西省

注意在前文中我加粗的那個詞“加入”這意味著上文中重點區(qū)域突出顯示的山西省圖層是經(jīng)過篩選而添加到地圖中的，而不是對山西省這個范圍進行符號化，如果我們要只顯示山西省的話，就將上述代碼中顯示全國范圍的代碼注釋掉即可，參見下方示例代碼

fig?=?plt.figure(figsize=(6,6))?? ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? #ax.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,facecolor='none',edgecolor='k',linewidth=1) ax.set_extent(extent,?crs=proj) for?state,?geos?in?zip(shapereader.Reader(shp_path).records(),?shapereader.Reader(shp_path).geometries())?: ????if?state.attributes['SHENG_ID']?in?[14]: ????????ax.add_geometries([geos],?crs=proj,facecolor='None',?edgecolor='red')??#?重點區(qū)域上色

添加南海小地圖

在之前的的學(xué)習(xí)中我們知道，cartopy繪制的地圖稱為子圖，在繪制中國地圖時候，有時候由于地圖大小的限制，我們無法展示部分地區(qū)如南海，常規(guī)的方法是繪制兩幅地圖，比如一張為全國地圖，一張為局部地圖，也就是常說的南海小地圖。

常見的subplot和subplot2grid函數(shù)一般來說繪制的地圖大小是一樣的，不容易展示比例大小，所以我們選擇add_axes()命令來繪制兩個大小不一樣的子圖。需要注意的是在繪制時我們需要定義兩個extent參數(shù)，即分別為總的地圖和南海小地圖分別定義畫布大小繪圖范圍

extent?=?[105,?133,?15,?45]#限定繪圖范圍 fig?=?plt.figure(figsize=(6,8))??#注意畫布大小 ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? ax.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,facecolor='none',edgecolor='k',linewidth=1) ax.set_extent(extent,?crs=proj) left,?bottom,?width,?height?=?0.67,?0.16,?0.23,?0.27#南海小地圖大小 ax2?=?fig.add_axes([left,?bottom,?width,?height],?projection=proj) ax2.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,?facecolor='none',edgecolor='k',linewidth=1) ax2.set_extent([105,?125,?0,?25])#注意南海小地圖的范圍

添加地形圖背景

一個簡單的嘗試

在看完了上述內(nèi)容后還是覺得有些不足，讓我們結(jié)合來Cartopy繪圖入門中所講述的內(nèi)容來給他添加一個背景吧

import?numpy?as?np import?cartopy.crs?as?ccrs import?matplotlib.pyplot?as?plt import?cartopy.feature?as?cfeature import?cartopy.io.shapereader?as?shapereader shp_path?=?"D:/data2/china/china/province_9south.shp"#切記路徑一定要英文，血的教訓(xùn)?。?！ fig?=?plt.figure(figsize=(6,8))??#注意畫布大小 proj?=?ccrs.PlateCarree()#創(chuàng)建投影，選擇cartopy的默認投影 ax???=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})??#子圖 extent?=?[105,?133,?15,?45]#限定繪圖范圍 ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? ax.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,facecolor='none',edgecolor='k',linewidth=1) ax.set_extent(extent,?crs=proj) #添加其他要素 ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) ax.add_feature(cfeature.RIVERS.with_scale('50m')) ax.add_feature(cfeature.LAKES.with_scale('50m')) left,?bottom,?width,?height?=?0.67,?0.16,?0.23,?0.27#南海小地圖大小 ax2?=?fig.add_axes([left,?bottom,?width,?height],?projection=proj) ax2.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,?facecolor='none',edgecolor='k',linewidth=1) ax2.set_extent([105,?125,?0,?25])#注意南海小地圖的范圍

這時候我們會疑惑，為什么南海小地圖還是這樣。在上文中提到過，南海小地圖是另一幅地圖，他與原來的中國地圖你可以理解為是兩幅地圖。所以也需要在南海小地圖中添加要素，也就是ax2中

import?numpy?as?np import?cartopy.crs?as?ccrs import?matplotlib.pyplot?as?plt import?cartopy.feature?as?cfeature import?cartopy.io.shapereader?as?shapereader shp_path?=?"D:/data2/china/china/province_9south.shp"#切記路徑一定要英文，血的教訓(xùn)?。?！ fig?=?plt.figure(figsize=(6,8))??#注意畫布大小 proj?=?ccrs.PlateCarree()#創(chuàng)建投影，選擇cartopy的默認投影 ax???=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})??#子圖 extent?=?[105,?133,?15,?45]#限定繪圖范圍 ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? ax.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,facecolor='none',edgecolor='k',linewidth=1) ax.set_extent(extent,?crs=proj) #添加其他要素_中國地圖 ax.add_feature(cfeature.OCEAN.with_scale('50m')) ax.add_feature(cfeature.LAND.with_scale('50m')) ax.add_feature(cfeature.RIVERS.with_scale('50m')) ax.add_feature(cfeature.LAKES.with_scale('50m')) left,?bottom,?width,?height?=?0.67,?0.16,?0.23,?0.27#南海小地圖大小 ax2?=?fig.add_axes([left,?bottom,?width,?height],?projection=proj) ax2.add_geometries(shapereader.Reader(shp_path).geometries(),?crs=proj,?facecolor='none',edgecolor='k',linewidth=1) ax2.set_extent([105,?125,?0,?25])#注意南海小地圖的范圍 #添加其他要素_南海小地圖 ax2.add_feature(cfeature.OCEAN.with_scale('50m')) ax2.add_feature(cfeature.LAND.with_scale('50m')) ax2.add_feature(cfeature.RIVERS.with_scale('50m')) ax2.add_feature(cfeature.LAKES.with_scale('50m'))

添加圖片作為背景

方法一

使用cartopy的add_image方法，這里以添加osm在線地圖的圖片為例子

imagery?=?OSM() ax.add_image(imagery,?4)#四為影像級別

完整代碼如下

import?cartopy.crs?as?ccrs import?cartopy.feature?as?cfeat from?cartopy.io.shapereader?import?Reader import?matplotlib.pyplot?as?plt from?matplotlib.image?import?imread from?cartopy.io.img_tiles?import?OSM def?create_map(): ????#?--設(shè)置shp路徑，數(shù)據(jù)集已公開 ????shp_path?=?'D:/data2/china/china/province_9south.shp' ????#?--設(shè)置tif路徑，數(shù)據(jù)集已公開 ????tif_path?=?'D:/data2/china/china/NE1_50M_SR_W.tif' ????#?--創(chuàng)建畫圖空間 ????proj?=?ccrs.PlateCarree()??#?創(chuàng)建坐標系 ????fig?=?plt.figure(figsize=(6,?8))??#?創(chuàng)建頁面 ????ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? ????#?--設(shè)置地圖屬性 ????provinces?=?cfeat.ShapelyFeature(Reader(shp_path).geometries(),proj,?edgecolor='k',facecolor='none') ????#?--加載省界線 ????ax.add_feature(provinces,?lw=0.6,?zorder=2) ????ax.set_extent([105,?133,?15,?45])??#可根據(jù)需求自行定義 ???#添加osm圖片 ????imagery?=?OSM() ????ax.add_image(imagery,?4)#四為影像級別 ax?=?create_map() plt.show()

方法二

添加圖片所使用的第二種方法是 ax.imshow，是基于matplotlib的，這里以加載地形圖圖片為例

ax.imshow(imread(tif_path),extent=[-180,?180,?-90,?90])注意范圍

image-20220720154440290

繪制的完整代碼如下：

import?cartopy.crs?as?ccrs import?cartopy.feature?as?cfeat from?cartopy.io.shapereader?import?Reader import?matplotlib.pyplot?as?plt from?matplotlib.image?import?imread def?create_map(): ????#?--設(shè)置shp路徑，數(shù)據(jù)集已公開 ????shp_path?=?'D:/data2/china/china/province_9south.shp' ????#?--設(shè)置tif路徑，數(shù)據(jù)集已公開 ????tif_path?=?'D:/data2/china/china/NE1_50M_SR_W.tif' ????#?--創(chuàng)建畫圖空間 ????proj?=?ccrs.PlateCarree()??#?創(chuàng)建坐標系 ????fig?=?plt.figure(figsize=(6,?8))??#?創(chuàng)建頁面 ????ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? ????#?--設(shè)置地圖屬性 ????provinces?=?cfeat.ShapelyFeature(Reader(shp_path).geometries(),proj,?edgecolor='k',facecolor='none') ????#?--加載省界線 ????ax.add_feature(provinces,?lw=0.6,?zorder=2) ????ax.set_extent([105,?133,?15,?45])??#可根據(jù)需求自行定義 ????#?ax.stock_img()??##可直接使用低分辨率背景 ????#?--加載高分辨率地形 ????ax.imshow(imread(tif_path),extent=[-180,?180,?-90,?90]) ax?=?create_map() plt.show()

你可以試試結(jié)合上述所講的內(nèi)容，做出一副這樣的地圖

image-20220720154648087

示例代碼如下：

import?cartopy.crs?as?ccrs import?cartopy.feature?as?cfeat from?cartopy.io.shapereader?import?Reader import?matplotlib.pyplot?as?plt from?matplotlib.image?import?imread def?create_map(): ????#?--設(shè)置shp路徑，數(shù)據(jù)集已公開 ????shp_path?=?'D:/data2/china/china/province_9south.shp' ????#?--設(shè)置tif路徑，數(shù)據(jù)集已公開 ????tif_path?=?'D:/data2/china/china/NE1_50M_SR_W.tif' ????#?--創(chuàng)建畫圖空間 ????proj?=?ccrs.PlateCarree()??#?創(chuàng)建坐標系 ????fig?=?plt.figure(figsize=(6,?8))??#?創(chuàng)建頁面 ????ax?=?fig.subplots(1,?1,?subplot_kw={'projection':?proj})?? ????#?--設(shè)置地圖屬性 ????provinces?=?cfeat.ShapelyFeature(Reader(shp_path).geometries(),proj,?edgecolor='k',facecolor='none') ????#?--加載省界線 ????ax.add_feature(provinces,?lw=0.6,?zorder=2) ????ax.set_extent([105,?133,?15,?45])??#可根據(jù)需求自行定義 ????#?ax.stock_img()??##可直接使用低分辨率背景 ????#?--加載高分辨率地形 ????ax.imshow(imread(tif_path),origin='upper',?transform=proj,extent=[-180,?180,?-90,?90]) ????#?--設(shè)置網(wǎng)格點屬性 ????gl?=?ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=True,linewidth=1.2,color='k',alpha=0.5,linestyle='--') ????#?--設(shè)置南海子圖 ????left,?bottom,?width,?height?=?0.67,?0.16,?0.23,?0.27 ????ax2?=?fig.add_axes([left,?bottom,?width,?height],?projection=proj) ????ax2.add_feature(provinces,?linewidth=0.6,?zorder=2) ????ax2.set_extent([105,?125,?0,?25]) ????#?ax2.stock_img()??##可直接使用低分辨率背景 ????#?--加載高分辨率地形 ????ax2.imshow(imread(tif_path),origin='upper',transform=proj,extent=[-180,?180,?-90,?90]) ????return?ax ax?=?create_map() plt.show()

注：本例教程參考至和鯨大佬【大自然搬運工】，針對大佬的代碼做了一些簡化方便大家理解；大佬繪制經(jīng)緯網(wǎng)好像使用的是matplotlib庫，與我之前所講述的cartopy庫繪制不一樣。原文章如下：

https://www.heywhale.com/mw/project/5f3c95a3af3980002cbf560b

學(xué)習(xí)心得

本文作為基礎(chǔ)入門教程就寫到這里吧，這幾天學(xué)習(xí)強度有點大了。在本文的學(xué)習(xí)過程中，我在和鯨社區(qū)找到了大量優(yōu)質(zhì)的學(xué)習(xí)博文，十分建議大家可以去看一看，同時以也才cartopy的文檔里得到了非常多的幫助，里面還有很多地圖的繪制方式，如果有機會，我希望我可以學(xué)習(xí)一下。

和鯨社區(qū)

https://www.heywhale.com/home

cartopyAPI參考

https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.mpl.geoaxes.GeoAxes.html#cartopy.mpl.geoaxes.GeoAxes.add_feature

cartopy更多地圖

https://scitools.org.uk/cartopy/docs/latest/gallery/index.html

數(shù)據(jù)下載

本文PDF下載學(xué)習(xí)

https://wwc.lanzouw.com/i71EE08450xg

行政區(qū)劃_省界.shp

https://wwc.lanzouw.com/iQ4TT0843gpc

GMT 中文社區(qū)提供的數(shù)據(jù)

https://wwc.lanzouw.com/iYLrc0843gvi

地形圖圖片

https://wwc.lanzouw.com/ipPrM0843g8f

參考文檔

https://zhajiman.github.io/post/cartopy_introduction/ https://mp.weixin.qq.com/s/FKNBOD2Yo4sawWTO76zZgA https://mp.weixin.qq.com/s/m4Ao0--o1umSICEWO9UlKw https://mp.weixin.qq.com/s/oe2ncVdqfdjaPSDZ2-YLQg https://cloud.tencent.com/developer/article/1790590 https://mp.weixin.qq.com/s/wbyrkmVoaWgz6MdBI_6BZg https://mp.weixin.qq.com/s/5RiuM2AQxNIvgNcLf5fCDw https://mp.weixin.qq.com/s/FGCPIoiC5OUkGGf_IreWZQ https://www.cnblogs.com/youxiaogang/p/14247184.html https://www.heywhale.com/mw/project/629eefe52d8168866e54adc9 https://www.heywhale.com/mw/project/5f1a4df794d484002d2db14a https://www.heywhale.com/mw/project/5f3c95a3af3980002cbf560b https://blog.csdn.net/m0_37362454/article/details/81511427 https://www.osgeo.cn/pygis/cartopy-feature.html https://vimsky.com/examples/detail/python-module-cartopy.mpl.ticker.html https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html https://mp.weixin.qq.com/s/jpZOpnFvMwi4ZTafflXIzw https://gnss.help/2018/04/24/cartopy-gallery/index.html https://blog.csdn.net/weixin_42372313/article/details/113572786 https://scitools.org.uk/cartopy/docs/latest/index.html https://scitools.org.uk/cartopy/docs/v0.13/matplotlib/geoaxes.html https://waterhackweek.github.io/visualization/04-geopandas-and-cartopy/ https://scitools.org.uk/cartopy/docs/v0.14/matplotlib/feature_interface.html https://mp.weixin.qq.com/s/ASV-VI6gfbea90vtJbdpIw

標簽：