1. 程式人生 > >python繪製地圖的利器Cartopy使用說明

python繪製地圖的利器Cartopy使用說明

python繪製地圖一般使用Basemap繪圖包,但該包配置相對較繁瑣,自定義性不強,這裡介紹一個繪製地圖的利器Cartopy,個人認為該工具方便、快捷,附上一些自己寫的程式。
準備工作,工欲善其事,必先利其器
(1)先下載主角:Cartopy
a)下載地址:
linux平臺直接去官網下載:http://scitools.org.uk/cartopy/download.html
windows平臺下的Cartopy下載地址:
http://www.lfd.uci.edu/~gohlke/pythonlibs/#cartopy
還需要一些必須的軟體包:Shapely、pyshp也都從上面的網址下載
官網自帶的示例網址:
http://scitools.org.uk/cartopy/docs/latest/gallery.html


b)Cartopy下載安裝完畢後,開啟python,輸入以下程式碼:
  1. import matplotlib.pyplot as plt
  2. import cartopy.crs as ccrs
  3. plt.figure(figsize=(6, 3))
  4. ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
  5. ax.coastlines(resolution='110m')
  6. ax.gridlines()
複製程式碼 執行後,如果順利匯入了cartopy包,不管有沒有出現地圖,都不要緊,第二步走起。
(2)再下載配角:地圖資料
下載地址:
http://www.naturalearthdata.com/downloads/

裡面有三種解析度的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資料夾,具體該資料夾在哪,搜尋電腦檔案找找看
),我安裝的是winpython2.7.9.3,physical目錄就位於上面這個目錄中,所以我把所有shape地圖檔案拷貝到了該physical目錄下。

地圖檔案列表.jpg (40.09 KB, 下載次數: 3)

下載附件  儲存到相簿

地圖書檔案列表

2015-3-25 19:40 上傳



準備工作完成後,進入正題:
(3)繪製地圖
前面兩步雖有些繁瑣,但會一勞永逸,下面是示例程式,scale引數用於調整使用哪種解析度的地圖,全球建議用1:110的,小區域可以用1:50的或1:10的。
a)繪製全球地圖程式:
  1. #===================================================
  2. #使用cartopy繪製地圖
  3. #需要從http://www.naturalearthdata.com/downloads/下載shape檔案
  4. #下載後,解壓縮,檔名統一去掉"ne_"開頭,拷貝至D:\Program Files\
  5. #WinPython-32bit-2.7.9.3\settings\.local\share\cartopy\shapefiles\natural_earth\physical\
  6. #路徑下面,coastline檔案對應ax.coastlines命令,land檔案對應land命令
  7. #===================================================
  8. scale='110m'
  9. import matplotlib.pyplot as plt
  10. import cartopy.crs as ccrs
  11. import cartopy.feature as cfeature
  12. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
  13. fig=plt.figure(figsize=(8, 10))
  14. ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
  15. ax.set_global()
  16. #===================================================
  17. #需要填充陸地顏色時使用
  18. #ax.add_feature(cfeature.LAND, facecolor='0.75') #預設為110m,其它解析度需用下面命令
  19. land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  20.                                                               facecolor=cfeature.COLORS['land'])
  21. ax.add_feature(land, facecolor='0.75')
  22. #===================================================
  23. #改變ax.add_feature和ax.coastlines的先後使用順序可實現邊界線的顯示或完全填充覆蓋
  24. ax.coastlines(scale)
  25. #===================================================
  26. #標註座標軸
  27. ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
  28. ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
  29. #zero_direction_label用來設定經度的0度加不加E和W
  30. lon_formatter = LongitudeFormatter(zero_direction_label=False)
  31. lat_formatter = LatitudeFormatter()
  32. ax.xaxis.set_major_formatter(lon_formatter)
  33. ax.yaxis.set_major_formatter(lat_formatter)
  34. #新增網格線
  35. gl = ax.gridlines()
複製程式碼

地圖1.png (49.05 KB, 下載次數: 5)

下載附件  儲存到相簿

2015-3-25 19:54 上傳



b)繪製全球地圖(函式形式)
  1. #===================================================
  2. #函式形式,呼叫cartopy,繪製全球地圖
  3. #===================================================
  4. import matplotlib.pyplot as plt
  5. import cartopy.crs as ccrs
  6. import cartopy.feature as cfeature
  7. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
  8. def make_map(scale):
  9.     fig=plt.figure(figsize=(8, 10))
  10.     ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
  11.     ax.set_global()
  12.     land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  13.                                                               facecolor=cfeature.COLORS['land'])
  14.     ax.add_feature(land, facecolor='0.75')
  15.     ax.coastlines(scale)
  16.     #標註座標軸
  17.     ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
  18.     ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
  19.     #zero_direction_label用來設定經度的0度加不加E和W
  20.     lon_formatter = LongitudeFormatter(zero_direction_label=False)
  21.     lat_formatter = LatitudeFormatter()
  22.     ax.xaxis.set_major_formatter(lon_formatter)
  23.     ax.yaxis.set_major_formatter(lat_formatter)
  24.     #新增網格線
  25.     #gl = ax.gridlines()
  26.     ax.grid()
  27.     return fig,ax
  28. fig,ax=make_map(scale='110m')
複製程式碼 此程式結果與上面一樣。
c)繪製區域地圖(函式形式)
  1. #===================================================
  2. #函式形式,呼叫cartopy,繪製區域地圖
  3. #===================================================
  4. import numpy as np
  5. import matplotlib.pyplot as plt
  6. import cartopy.crs as ccrs
  7. import cartopy.feature as cfeature
  8. from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
  9. def make_map(scale,box,xstep,ystep):
  10.     fig=plt.figure(figsize=(8, 10))
  11.     ax=plt.axes(projection=ccrs.PlateCarree())
  12.     #set_extent需要配置相應的crs,否則出來的地圖範圍不準確
  13.     ax.set_extent(box,crs=ccrs.PlateCarree())
  14.     land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  15.                                                               facecolor=cfeature.COLORS['land'])
  16.     ax.add_feature(land, facecolor='0.75')
  17.     ax.coastlines(scale)
  18.     #===================================================
  19.     #影象地址D:\Program Files\WinPython-32bit-2.7.9.3\python-2.7.9\Lib\site-packages\
  20.     #cartopy\data\raster\natural_earth\50-natural-earth-1-downsampled.png
  21.     #如果有其它高精度影象檔案,改名替換即可
  22.     ax.stock_img()
  23.     #===================================================
  24.     #標註座標軸
  25.     ax.set_xticks(np.arange(box[0],box[1]+xstep,xstep), crs=ccrs.PlateCarree())
  26.     ax.set_yticks(np.arange(box[2],box[3]+ystep,ystep), crs=ccrs.PlateCarree())
  27.     #zero_direction_label用來設定經度的0度加不加E和W
  28.     lon_formatter = LongitudeFormatter(zero_direction_label=False)
  29.     lat_formatter = LatitudeFormatter()
  30.     ax.xaxis.set_major_formatter(lon_formatter)
  31.     ax.yaxis.set_major_formatter(lat_formatter)
  32.     #新增網格線
  33.     ax.grid()
  34.     return fig,ax
  35. box=[100,150,0,50]
  36. fig,ax=make_map(scale='50m',box=box,xstep=10,ystep=10)
複製程式碼

地圖2.png (149.94 KB, 下載次數: 2)

下載附件  儲存到相簿

2015-3-25 19:57 上傳



(4)繪製極地投影地圖
繪製該地圖需要使用最新版的cartopy-0.12版本,windows下暫無此版本,可以安裝好0.11後,刪除cartopy安裝目錄下的mpl資料夾下的全部檔案,然後拷貝0.12目錄cartopy-0.12.0rc1\lib\cartopy\mpl資料夾下的全部檔案進行替換,即可使用新版的功能。此程式長了點,因為Cartopy對極座標投影沒有自動標註經緯度功能,需要自己設定,調了好久標註位置,請大家切用且珍惜哈。
  1. import matplotlib.path as mpath
  2. import matplotlib.pyplot as plt
  3. import numpy as np
  4. import matplotlib.ticker as mticker
  5. import cartopy.crs as ccrs
  6. import cartopy.feature as cfeature
  7. fig=plt.figure(figsize=(6, 6))
  8. ax=plt.axes(projection=ccrs.NorthPolarStereo())
  9. box=[-180, 180, 55, 90]
  10. xstep,ystep=30,15
  11.     # Limit the map to -60 degrees latitude and below.
  12. ax.set_extent(box, crs=ccrs.PlateCarree())
  13. scale='50m'
  14. land = cfeature.NaturalEarthFeature('physical', 'land', scale,edgecolor='face',
  15.                                                               facecolor=cfeature.COLORS['land'])
  16. ocean = cfeature.NaturalEarthFeature('physical', 'ocean', scale,edgecolor='face',
  17.                                                              facecolor=cfeature.COLORS['water'])
  18. ax.add_feature(land,facecolor='0.75')
  19. ax.add_feature(ocean,facecolor='blue')
  20. ax.coastlines(scale,linewidth=0.9)
  21. #標註座標軸
  22. line=ax.gridlines(draw_labels=False)
  23. line.ylocator=mticker.FixedLocator(np.arange(40,90,20))#手動設定x軸刻度
  24. line.xlocator=mticker.FixedLocator(np.arange(-180,210,30))#手動設定x軸刻度
  25.     # Compute a circle in axes coordinates, which we can use as a boundary
  26.     # for the map. We can pan/zoom as much as we like - the boundary will be
  27.     # permanently circular.
  28. theta = np.linspace(0, 2*np.pi, 100)
  29. center, radius = [0.5, 0.5], 0.5
  30. verts = np.vstack([np.sin(theta), np.cos(theta)]).T
  31. circle = mpath.Path(verts * radius + center)
  32. ax.set_boundary(circle, transform=ax.transAxes)
  33. #建立要標註的labels字串
  34. ticks=np.arange(0,210,30)
  35. etick=['0']+['%d$^\circ$E'%tick for tick in ticks if (tick !=0) & (tick!=180)]+['180']
  36. wtick=['%d$^\circ$W'%tick for tick in ticks if (tick !=0) & (tick!=180)]
  37. labels=etick+wtick
  38. #建立與labels對應的經緯度標註位置
  39. #xticks=[i for i in np.arange(0,210,30)]+[i for i in np.arange(-32,-180,-30)]
  40. xticks=[-0.8,28,58,89.1,120,151,182.9,-36,-63,-89,-114,-140]
  41. yticks=[53]+[53]+[54]+[55]*2+[54.5]+[54]+[50]+[49]*3+[50.6]
  42. #標註經緯度   
  43. #ax.text(0.01,0.23,'60$^\circ$W',transform=ax.transAxes,rotation=25)
  44. #ax.text(-63,50,'60$^\circ$W',transform=ccrs.Geodetic(),rotation=25)
  45. for xtick,ytick,label in zip(xticks,yticks,labels):
  46.     ax.text(xtick,ytick,label,transform=ccrs.Geodetic())
  47. x=[180, 180, 0, 0]
  48. y=[50, 90, 90, 50]
  49. ax.plot([-180,0],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.4)
  50. ax.plot([-90,90],[80,80],':',transform=ccrs.Geodetic(),color='k',linewidth=0.5)
  51. #ax.plot([90,0],[50,50],'-.',transform=ccrs.Geodetic(),color='r',linewidth=6)
  52. ax.text(11.9333,78.9166,r'$\bigstar,transform=ccrs.Geodetic(),size=15,color='r')
  53. fig.savefig(u'c:\\北極.png',dpi=300)
複製程式碼

北極地圖.png (64.49 KB, 下載次數: 3)

下載附件  儲存到相簿

北極投影地圖

2015-3-25 20:24 上傳



介紹完畢,寫一個帖子還真挺耗時間,不多說了,更多的功能,請去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即呼叫到了自己的地圖,不用擔心邊界不準確了。
不知道為什麼最後多出來一張圖片,麻煩版主幫忙刪了吧。