Cartopy繪圖系列 | 繪制全球溫度場+風矢量場

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

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

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

cartopy.png

Cartopy簡介與安裝

Cartopy 是一個開源免費的第三方 Python 擴展包，由英國氣象辦公室的科學家們開發(fā)，支持 Python 2.7 和 Python 3，致力于使用最簡單直觀的方式生成地圖，并提供對 matplotlib 的接口，兩者可以完美的協(xié)作。
1、Cartopy常用于用于地理空間數(shù)據(jù)處理，以便生成地圖和其他地理空間數(shù)據(jù)分析。Cartopy利用了強大的PROJ.4、NumPy和Shapely庫，并在Matplotlib之上提供了一個編程接口，用于創(chuàng)建出版高質量的地圖。
2、Cartopy的關鍵特性是它面向對象的投影定義，以及在這些投影之間轉換點、線、向量、多邊形和圖像的能力。

為什么要選用Cartopy ?

Cartopy安裝：

官網教程：https://scitools.org.uk/cartopy/docs/latest/installing.html#installing
1）Anaconda環(huán)境
如果你正在使用 Python 的科學計算發(fā)行版 Anaconda，安裝 Cartopy 非常容易。
命令行輸入: conda install -c conda-forge cartopy
2）Windows環(huán)境
命令行輸入：pip install cartopy
3）可以切換國內像源，安裝速度更快
命令行輸入：pip install -i https://pypi.tuna.tsinghua.edu.cn/simple cartopy

測試是否安裝成功：
啟動python命令行工具，輸入 import cartopy 如果沒有報錯，則安裝成功！

Cartopy 制圖簡介

  • Cartopy官方網站中列舉了許多繪圖案例，并且有完整的代碼演示。鏈接https://scitools.org.uk/cartopy/docs/latest/gallery/index.html

    image.png

  • Cartopy地圖繪制的總體流程就是：(1) 選擇合適的投影； (2) 添加地圖要素（自定義shp\海岸線\邊界線等）；(3) 疊加空間數(shù)據(jù)；(4) 設置地圖要素（經緯度標簽、比例尺、文本等）

下面介紹幾種作者自己總結的常見繪圖案例

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='-') # 添加國家邊界
    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()

結果展示

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():

    # 設置投影
    proj = ccrs.PlateCarree() # 圓柱投影， 默認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')

    # 設置經緯網以及標簽
    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()

結果展示

test.jpg

3、綜合案例：Cartopy繪制全球溫度場+風矢量場數(shù)據(jù)**

NCEP/NCAR Reanalysis 1 再分析資料的地面10m風場數(shù)據(jù)（u v）和地面溫度數(shù)據(jù)繪圖示例。（以2020年01月01日為例）
數(shù)據(jù)可在官網下載：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風速數(shù)據(jù)繪制風矢量場圖
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()

結果展示

test.jpg

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

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

相關閱讀更多精彩內容

友情鏈接更多精彩內容