世上總會(huì)有“猝不及防的再見(jiàn)和毫無(wú)留戀的散場(chǎng)”肃弟， 但也會(huì)有 “突如其來(lái)的遇見(jiàn)和始料未及的歡喜”氨鹏， 無(wú)論何時(shí)争拐，記得不要失了好心情粮坞！
空間數(shù)據(jù)的可視化展示是空間數(shù)據(jù)可視化制圖的常見(jiàn)需求。目前常見(jiàn)的地圖繪圖工具和軟件有:
- (1) 支持 多種操作系統(tǒng)的命令行工具 GMT（Generic Mapping Tools）
- (2) ArcGIS赤拒、GrDAS等秫筏。
- (3) NCAR Command Language): 常用于氣象數(shù)據(jù)讀取和可視化分析的高級(jí)語(yǔ)言（目前也已經(jīng)遷移到PyNGL上）
- (4)Python 繪圖工具 matplotlib 的擴(kuò)展包 Basemap（作者在前面的文章中有簡(jiǎn)單介紹）
除此之外，小編想給大家推薦 Python 的一種制圖工具包 Cartopy,（內(nèi)容比Basemap更加豐富和實(shí)用）
Cartopy簡(jiǎn)介與安裝
Cartopy 是一個(gè)開(kāi)源免費(fèi)的第三方 Python 擴(kuò)展包挎挖，由英國(guó)氣象辦公室的科學(xué)家們開(kāi)發(fā)这敬，支持 Python 2.7 和 Python 3，致力于使用最簡(jiǎn)單直觀的方式生成地圖蕉朵，并提供對(duì) matplotlib 的接口崔涂，兩者可以完美的協(xié)作。
1始衅、Cartopy常用于用于地理空間數(shù)據(jù)處理堪伍，以便生成地圖和其他地理空間數(shù)據(jù)分析。Cartopy利用了強(qiáng)大的PROJ.4觅闽、NumPy和Shapely庫(kù)，并在Matplotlib之上提供了一個(gè)編程接口涮俄，用于創(chuàng)建出版高質(zhì)量的地圖蛉拙。
2、Cartopy的關(guān)鍵特性是它面向?qū)ο蟮耐队岸x彻亲，以及在這些投影之間轉(zhuǎn)換點(diǎn)孕锄、線吮廉、向量、多邊形和圖像的能力畸肆。
為什么要選用Cartopy ?
- Cartopy 的工作流非常簡(jiǎn)單：設(shè)置地圖投影宦芦，添加地圖要素，最后顯示地圖轴脐。
- Basemap繪圖包配置相對(duì)較繁瑣调卑，自定義性不強(qiáng)而且Basemap 將在 2020 年版本的更新。因此大咱，如果你在尋找一個(gè)新的Python制圖工具包恬涧，Cartopy是Basemap官網(wǎng)欽定的繼承人，無(wú)疑是最佳選擇（https://matplotlib.org/basemap/users/intro.html#cartopy-new-management-and-eol-announcement）碴巾。
Cartopy安裝：
官網(wǎng)教程：https://scitools.org.uk/cartopy/docs/latest/installing.html#installing
1）Anaconda環(huán)境
如果你正在使用 Python 的科學(xué)計(jì)算發(fā)行版 Anaconda溯捆，安裝 Cartopy 非常容易。
命令行輸入: conda install -c conda-forge cartopy
2）Windows環(huán)境
命令行輸入：pip install cartopy
3）可以切換國(guó)內(nèi)像源厦瓢，安裝速度更快
命令行輸入：pip install -i https://pypi.tuna.tsinghua.edu.cn/simple cartopy
測(cè)試是否安裝成功：
啟動(dòng)python命令行工具提揍，輸入 import cartopy 如果沒(méi)有報(bào)錯(cuò)，則安裝成功煮仇！
Cartopy 制圖簡(jiǎn)介
-
Cartopy官方網(wǎng)站中列舉了許多繪圖案例劳跃，并且有完整的代碼演示。鏈接https://scitools.org.uk/cartopy/docs/latest/gallery/index.html
image.png Cartopy地圖繪制的總體流程就是：(1) 選擇合適的投影欺抗； (2) 添加地圖要素（自定義shp\海岸線\邊界線等）售碳；(3) 疊加空間數(shù)據(jù)；(4) 設(shè)置地圖要素（經(jīng)緯度標(biāo)簽绞呈、比例尺贸人、文本等）
下面介紹幾種作者自己總結(jié)的常見(jiàn)繪圖案例
1、全球地圖顯示
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : cartopy_script.py
# Author : zengsk in NanJing
# Created: 2020/6/15 0:54
"""
Details: Cartopy package 繪圖示例
"""
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
def main():
fig = plt.figure(figsize=(8, 8))
# Label axes of a Plate Carree projection with a central longitude of 180:
ax1 = fig.add_subplot(2, 1, 1,
projection=ccrs.PlateCarree(central_longitude=180))
ax1.set_global()
ax1.coastlines() # 添加海岸線
ax1.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
ax1.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
ax2 = fig.add_subplot(2, 1, 2,
projection=ccrs.Robinson(central_longitude=0))
# make the map global rather than have it zoom in to
# the extents of any plotted data
ax2.set_global()
ax2.stock_img()
ax2.coastlines()
ax2.add_feature(cfeature.BORDERS, linestyle='-') # 添加國(guó)家邊界
ax2.plot(-0.08, 51.53, 'o', transform=ccrs.PlateCarree())
plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
plt.show()
if __name__ == '__main__':
main()
結(jié)果展示
2佃声、顯示自定義 shapefile文件**
往往我們需要繪制自定義范圍的研究區(qū)域艺智，需要繪制指定shapefile文件的邊界！
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : cartopy_test2.py
"""
Created by s.k zeng in Chengdu on 2020/07/07
"""
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as mticker
import matplotlib.pyplot as plt
def main():
# 設(shè)置投影
proj = ccrs.PlateCarree() # 圓柱投影圾亏， 默認(rèn)WGS1984
extent = [70, 140, 0, 60] # 顯示范圍
shpfn = r'E:\GisData\SHP\China_shp\bou2_4m\bou2_4l.shp'
tpshp = r'E:\GisData\SHP\Tibetan\TP_boundary\TP_Clip.shp'
glshp = r'E:\GisData\SHP\Global\country.shp'
reader = shpreader.Reader(shpfn)
states_provinces = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')
reader = shpreader.Reader(tpshp)
tpfeat = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')
reader = shpreader.Reader(glshp)
glfeat = cfeature.ShapelyFeature(reader.geometries(), proj, facecolor='none')
fig = plt.figure(figsize=(8, 8), frameon=True)
ax1 = fig.add_subplot(1, 1, 1, projection=proj)
# Label axes of a Mercator projection without degree symbols in the labels
# and formatting labels to include 1 decimal place:
ax1.set_extent(extent, proj)
ax1.add_feature(glfeat, linewidth=1.5, edgecolor='grey')
ax1.add_feature(states_provinces, linewidth=1.5, edgecolor='blue')
ax1.add_feature(tpfeat, linewidth=2.0, edgecolor='red')
# 設(shè)置經(jīng)緯網(wǎng)以及標(biāo)簽
ax1.gridlines(proj, draw_labels=False, linewidth=1.2, color='k', alpha=0.5, linestyle='--')
ax1.set_xticks(np.arange(extent[0], extent[1] + 10, 10), crs=proj)
ax1.set_yticks(np.arange(extent[2], extent[3] + 10, 10), crs=proj)
ax1.xaxis.set_major_formatter(mticker.LongitudeFormatter())
ax1.yaxis.set_major_formatter(mticker.LatitudeFormatter())
plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
plt.show()
if __name__ == '__main__':
main()
結(jié)果展示
3十拣、綜合案例：Cartopy繪制全球溫度場(chǎng)+風(fēng)矢量場(chǎng)數(shù)據(jù)**
NCEP/NCAR Reanalysis 1 再分析資料的地面10m風(fēng)場(chǎng)數(shù)據(jù)（u v）和地面溫度數(shù)據(jù)繪圖示例。（以2020年01月01日為例）
數(shù)據(jù)可在官網(wǎng)下載：https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.html
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name : uvwind_plot.py
"""
Details: NCEP/NCAR Reanalysis 1: Surface 10m風(fēng)速數(shù)據(jù)繪制風(fēng)矢量場(chǎng)圖
Created by s.k zeng in Chengdu on 2020/07/07
"""
import netCDF4 as nc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.mpl.ticker as mticker
def read_uv():
'''
:return: lon, lat, uwind, vwind
'''
uwind = r"E:\Scripts\Test\uvdata\uwnd.10m.gauss.2020.nc"
vwind = r"E:\Scripts\Test\uvdata\vwnd.10m.gauss.2020.nc"
ufid = nc.Dataset(uwind)
vfid = nc.Dataset(vwind)
print(ufid.variables)
lon = ufid.variables['lon'][:]
lat = ufid.variables['lat'][:]
u = ufid.variables['uwnd'][0, :, :] # unit: m/s
v = vfid.variables['vwnd'][0, :, :]
return lon, lat, u, v
def read_airtemper():
temperFn = r"E:\Scripts\Test\uvdata\air.sig995.2020.nc"
tfid = nc.Dataset(temperFn)
print(tfid.variables)
lon = tfid.variables['lon'][:]
lat = tfid.variables['lat'][:]
air = tfid.variables['air'][0, :, :] # air temperature unit: degK
return lon, lat, air
def main():
lon, lat, u, v = read_uv()
lon2, lat2, air = read_airtemper()
proj = ccrs.PlateCarree(central_longitude=180)
fig = plt.figure(figsize=(9, 6), dpi=300)
ax = fig.add_subplot(1, 1, 1, projection=proj)
cs = ax.contourf(lon2, lat2, air, transform=proj, cmap='RdBu_r') # RdBu_r nipy_spectral
ax.quiver(lon, lat, u, v, transform=ccrs.PlateCarree(), regrid_shape=35, width=0.002)
ax.coastlines(color='dimgray')
ax.set_global()
cbar = fig.colorbar(cs, orientation='vertical', pad=0.02, aspect=20, shrink=0.6)
cbar.set_label('degK')
xticks = [-180, -120, -60, 0, 60, 120, 180]
ax.set_xticks(xticks, crs=proj)
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], crs=proj)
lon_formatter = mticker.LongitudeFormatter(zero_direction_label=True)
lat_formatter = mticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
plt.savefig(r'C:\Users\zengsk\Desktop\test.jpg', dpi=300)
plt.show()
if __name__ == '__main__':
main()
結(jié)果展示
