python繪製地圖的利器Cartopy使用說明
阿新 • • 發佈:2019-01-22
準備工作,工欲善其事,必先利其器
(1)先下載主角:Cartopy
a)下載地址:
linux平臺直接去官網下載:http://scitools.org.uk/cartopy/download.html
windows平臺下的Cartopy下載地址:
http://www.lfd.uci.edu/~gohlke/pythonlibs/#cartopy
還需要一些必須的軟體包:Shapely、pyshp也都從上面的網址下載
官網自帶的示例網址:
b)Cartopy下載安裝完畢後,開啟python,輸入以下程式碼:
- import matplotlib.pyplot as plt
- import cartopy.crs as ccrs
- plt.figure(figsize=(6, 3))
- ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
- ax.coastlines(resolution='110m')
- ax.gridlines()
(2)再下載配角:地圖資料
下載地址:
裡面有三種解析度的shape地圖資料可選,方便起見,分別下載三種解析度中的physical資料中的Coastline和Land資料,每種資料下載後都是一個壓縮包,如下載了1:10解析度的physical中的coastline資料壓縮包:ne_10m_coastline.zip,解壓縮後有6個檔案,其中“ne_10m_coastline.README”和“ne_10m_coastline.VERSION”可直接刪除,剩餘4個,進行改名操作,副檔名前面的檔名,如“ne_10m_coastline”,修改為“10m_coastline”,即去掉“ne_”,4個檔案分別這樣更改。再下載1:50和1:110的檔案分別進行此操作。所有地圖檔案下載、解壓、更名完畢後,拷貝到一個資料夾下。我的資料夾列表如下圖,把這些檔案全選(注意進入資料夾目錄,全選檔案,不帶資料夾),複製貼上到D:\Program Files\WinPython-32bit-2.7.9.3\settings\.local\share\cartopy\shapefiles\natural_earth\physical 目錄下(該目錄根據自己所裝的python而定,執行(1)中的程式後,程式會自動建立physical資料夾,具體該資料夾在哪,搜尋電腦檔案找找看
地圖檔案列表.jpg
地圖書檔案列表
準備工作完成後,進入正題:
(3)繪製地圖
前面兩步雖有些繁瑣,但會一勞永逸,下面是示例程式,scale引數用於調整使用哪種解析度的地圖,全球建議用1:110的,小區域可以用1:50的或1:10的。
a)繪製全球地圖程式:
- #===================================================
- #使用cartopy繪製地圖
- #需要從http://www.naturalearthdata.com/downloads/下載shape檔案
- #下載後,解壓縮,檔名統一去掉"ne_"開頭,拷貝至D:\Program Files\
- #WinPython-32bit-2.7.9.3\settings\.local\share\cartopy\shapefiles\natural_earth\physical\
- #路徑下面,coastline檔案對應ax.coastlines命令,land檔案對應land命令
- #===================================================
- scale='110m'
- import matplotlib.pyplot as plt
- import cartopy.crs as ccrs
- import cartopy.feature as cfeature
- from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
- fig=plt.figure(figsize=(8, 10))
- ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
- ax.set_global()
- #===================================================
- #需要填充陸地顏色時使用
- #ax.add_feature(cfeature.LAND, facecolor='0.75') #預設為110m,其它解析度需用下面命令
- land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
- facecolor=cfeature.COLORS['land'])
- ax.add_feature(land, facecolor='0.75')
- #===================================================
- #改變ax.add_feature和ax.coastlines的先後使用順序可實現邊界線的顯示或完全填充覆蓋
- ax.coastlines(scale)
- #===================================================
- #標註座標軸
- ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
- ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
- #zero_direction_label用來設定經度的0度加不加E和W
- lon_formatter = LongitudeFormatter(zero_direction_label=False)
- lat_formatter = LatitudeFormatter()
- ax.xaxis.set_major_formatter(lon_formatter)
- ax.yaxis.set_major_formatter(lat_formatter)
- #新增網格線
- gl = ax.gridlines()
地圖1.png
b)繪製全球地圖(函式形式)
- #===================================================
- #函式形式,呼叫cartopy,繪製全球地圖
- #===================================================
- import matplotlib.pyplot as plt
- import cartopy.crs as ccrs
- import cartopy.feature as cfeature
- from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
- def make_map(scale):
- fig=plt.figure(figsize=(8, 10))
- ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
- ax.set_global()
- land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
- facecolor=cfeature.COLORS['land'])
- ax.add_feature(land, facecolor='0.75')
- ax.coastlines(scale)
- #標註座標軸
- ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
- ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
- #zero_direction_label用來設定經度的0度加不加E和W
- lon_formatter = LongitudeFormatter(zero_direction_label=False)
- lat_formatter = LatitudeFormatter()
- ax.xaxis.set_major_formatter(lon_formatter)
- ax.yaxis.set_major_formatter(lat_formatter)
- #新增網格線
- #gl = ax.gridlines()
- ax.grid()
- return fig,ax
- fig,ax=make_map(scale='110m')
c)繪製區域地圖(函式形式)
- #===================================================
- #函式形式,呼叫cartopy,繪製區域地圖
- #===================================================
- import numpy as np
- import matplotlib.pyplot as plt
- import cartopy.crs as ccrs
- import cartopy.feature as cfeature
- from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
- def make_map(scale,box,xstep,ystep):
- fig=plt.figure(figsize=(8, 10))
- ax=plt.axes(projection=ccrs.PlateCarree())
- #set_extent需要配置相應的crs,否則出來的地圖範圍不準確
- ax.set_extent(box,crs=ccrs.PlateCarree())
- land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
- facecolor=cfeature.COLORS['land'])
- ax.add_feature(land, facecolor='0.75')
- ax.coastlines(scale)
- #===================================================
- #影象地址D:\Program Files\WinPython-32bit-2.7.9.3\python-2.7.9\Lib\site-packages\
- #cartopy\data\raster\natural_earth\50-natural-earth-1-downsampled.png
- #如果有其它高精度影象檔案,改名替換即可
- ax.stock_img()
- #===================================================
- #標註座標軸
- ax.set_xticks(np.arange(box[0],box[1]+xstep,xstep), crs=ccrs.PlateCarree())
- ax.set_yticks(np.arange(box[2],box[3]+ystep,ystep), crs=ccrs.PlateCarree())
- #zero_direction_label用來設定經度的0度加不加E和W
- lon_formatter = LongitudeFormatter(zero_direction_label=False)
- lat_formatter = LatitudeFormatter()
- ax.xaxis.set_major_formatter(lon_formatter)
- ax.yaxis.set_major_formatter(lat_formatter)
- #新增網格線
- ax.grid()
- return fig,ax
- box=[100,150,0,50]
- fig,ax=make_map(scale='50m',box=box,xstep=10,ystep=10)
地圖2.png
(4)繪製極地投影地圖
繪製該地圖需要使用最新版的cartopy-0.12版本,windows下暫無此版本,可以安裝好0.11後,刪除cartopy安裝目錄下的mpl資料夾下的全部檔案,然後拷貝0.12目錄cartopy-0.12.0rc1\lib\cartopy\mpl資料夾下的全部檔案進行替換,即可使用新版的功能。此程式長了點,因為Cartopy對極座標投影沒有自動標註經緯度功能,需要自己設定,調了好久標註位置,請大家切用且珍惜哈。
- import matplotlib.path as mpath
- import matplotlib.pyplot as plt
- import numpy as np
- import matplotlib.ticker as mticker
- import cartopy.crs as ccrs
- import cartopy.feature as cfeature
- fig=plt.figure(figsize=(6, 6))
- ax=plt.axes(projection=ccrs.NorthPolarStereo())
- box=[-180, 180, 55, 90]
- xstep,ystep=30,15
- # Limit the map to -60 degrees latitude and below.
- ax.set_extent(box, crs=ccrs.PlateCarree())
- scale='50m'
- land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
- facecolor=cfeature.COLORS['land'])
- ocean = cfeature.NaturalEarthFeature('physical', 'ocean', scale,edgecolor='face',
- facecolor=cfeature.COLORS['water'])
- ax.add_feature(land,facecolor='0.75')
- ax.add_feature(ocean,facecolor='blue')
- ax.coastlines(scale,linewidth=0.9)
- #標註座標軸
- line=ax.gridlines(draw_labels=False)
- line.ylocator=mticker.FixedLocator(np.arange(40,90,20))#手動設定x軸刻度
- line.xlocator=mticker.FixedLocator(np.arange(-180,210,30))#手動設定x軸刻度
- # Compute a circle in axes coordinates, which we can use as a boundary
- # for the map. We can pan/zoom as much as we like - the boundary will be
- # permanently circular.
- theta = np.linspace(0, 2*np.pi, 100)
- center, radius = [0.5, 0.5], 0.5
- verts = np.vstack([np.sin(theta), np.cos(theta)]).T
- circle = mpath.Path(verts * radius + center)
- ax.set_boundary(circle, transform=ax.transAxes)
- #建立要標註的labels字串
- ticks=np.arange(0,210,30)
- etick=['0']+['%d$^\circ$E'%tick for tick in ticks if (tick !=0) & (tick!=180)]+['180']
- wtick=['%d$^\circ$W'%tick for tick in ticks if (tick !=0) & (tick!=180)]
- labels=etick+wtick
- #建立與labels對應的經緯度標註位置
- #xticks=[i for i in np.arange(0,210,30)]+[i for i in np.arange(-32,-180,-30)]
- xticks=[-0.8,28,58,89.1,120,151,182.9,-36,-63,-89,-114,-140]
- yticks=[53]+[53]+[54]+[55]*2+[54.5]+[54]+[50]+[49]*3+[50.6]
- #標註經緯度
- #ax.text(0.01,0.23,'60$^\circ$W',transform=ax.transAxes,rotation=25)
- #ax.text(-63,50,'60$^\circ$W',transform=ccrs.Geodetic(),rotation=25)
- for xtick,ytick,label in zip(xticks,yticks,labels):
- ax.text(xtick,ytick,label,transform=ccrs.Geodetic())
- x=[180, 180, 0, 0]
- y=[50, 90, 90, 50]
- ax.plot([-180,0],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.4)
- ax.plot([-90,90],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.5)
- #ax.plot([90,0],[50,50],'-.',transform=ccrs.Geodetic(),color='r',linewidth=6)
- ax.text(11.9333,78.9166,r'$\bigstar,transform=ccrs.Geodetic(),size=15,color='r')
- fig.savefig(u'c:\\北極.png',dpi=300)
北極地圖.png
北極投影地圖
介紹完畢,寫一個帖子還真挺耗時間,不多說了,更多的功能,請去cartopy官網http://scitools.org.uk/cartopy/docs/latest/gallery.html,看更多的例子。
另外一個notebook網址也很多示例可以借鑑:http://nbviewer.ipython.org/gist/pelson
注:如果想使用自己的地圖,比如我有一個帶我國法定邊界的地圖shape檔案,名為:World.shp、World.shx、World.dbf、World.prj,重新命名為10m_coastline或10m_land檔案(這裡要根據該shape檔案是line型還是polygon型),替換原目錄下的shape檔案,畫圖時使用scale為10m即呼叫到了自己的地圖,不用擔心邊界不準確了。
不知道為什麼最後多出來一張圖片,麻煩版主幫忙刪了吧。