Cartopy繪圖系列 | 繪制全球溫度場(chǎng)+風(fēng)矢量場(chǎng)

世上總會(huì)有“猝不及防的再見和毫無留戀的散場(chǎng)”， 但也會(huì)有 “突如其來的遇見和始料未及的歡喜”， 無論何時(shí)，記得不要失了好心情！

空間數(shù)據(jù)的可視化展示是空間數(shù)據(jù)可視化制圖的常見需求。目前常見的地圖繪圖工具和軟件有:
  • (1) 支持 多種操作系統(tǒng)的命令行工具 GMT（Generic Mapping Tools）
  • (2) ArcGIS、GrDAS等。
  • (3) NCAR Command Language): 常用于氣象數(shù)據(jù)讀取和可視化分析的高級(jí)語言（目前也已經(jīng)遷移到PyNGL上）
  • (4)Python 繪圖工具 matplotlib 的擴(kuò)展包 Basemap（作者在前面的文章中有簡(jiǎn)單介紹）

除此之外，小編想給大家推薦 Python 的一種制圖工具包 Cartopy,（內(nèi)容比Basemap更加豐富和實(shí)用）

cartopy.png

Cartopy簡(jiǎn)介與安裝

Cartopy 是一個(gè)開源免費(fèi)的第三方 Python 擴(kuò)展包，由英國(guó)氣象辦公室的科學(xué)家們開發(fā)，支持 Python 2.7 和 Python 3，致力于使用最簡(jiǎn)單直觀的方式生成地圖，并提供對(duì) matplotlib 的接口，兩者可以完美的協(xié)作。
1、Cartopy常用于用于地理空間數(shù)據(jù)處理，以便生成地圖和其他地理空間數(shù)據(jù)分析。Cartopy利用了強(qiáng)大的PROJ.4、NumPy和Shapely庫，并在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)欽定的繼承人，無疑是最佳選擇（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 如果沒有報(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é)的常見繪圖案例

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é)果展示

test.jpg

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é)果展示

test.jpg

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é)果展示

test.jpg

感謝支持，本人能力有限，如文章存在錯(cuò)誤和高見歡迎留言！

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容