炸鸡人博客 基本上无害
Posts with the tag cartopy:

Cartopy 系列:画中国地图的工具箱 frykit

contourf

前言

最早笔者用 Python 画中国地图时,会准备 bou2_4p.shp 文件,然后封装一个读取 shapefile 并添加到 GeoAxes 上的函数,别的项目要用时就把数据和函数复制粘贴过去。Cartopy 系列:从入门到放弃 里就是这么做的。

后来工作中用到了 Clarmy 开发的 cnmaps 包,只用两行就能快速绘制地图,非常方便。同时萌生了自己实现一个功能类似的包的想法,遂开发出了 frykit

Cartopy 系列:裁剪填色图出界问题

前言

裁剪或者说白化,就是让填色图只显示在多边形里面,不显示在多边形外面,例如只显示 GeoAxes.contourf 在中国境内的结果。实现方法为:

from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
from cartopy.io.shapereader import Reader

reader = Reader(filepath)
geom = next(reader.geometries())
reader.close()

cf = ax.contourf(X, Y, Z, transform=crs)
geom = ax.projection.project_geometry(geom, crs)
path = Path.make_compound_path(*geos_to_path(geom))
for col in cf.collections:
    col.set_clip_path(path, ax.transData)
  • crs 坐标系上的多边形对象变换到 data 坐标系上。
  • 利用 geos_to_pathmake_compound_path 将变换后的多边形转为 Path 对象。
  • QuadContourSet.collections 里的每个成员调用 set_clip_path 方法,并且指定 data 坐标系。

fig1

Cartopy 系列:利用多边形生成掩膜数组

二维平面上一系列点的坐标由 xy 数组描述,同时准备一个形状相同的 mask 数组。若第 i 个点落入了平面上一个多边形的内部,则令 mask[i] = True;若在多边形外,则令 mask[i] = False。由此得到的 mask 数组即掩膜(mask)数组,它可以作为布尔索引分出多边形内外的点

x_in, y_in = x[mask], y[mask]
x_out, y_out = x[mask], y[mask]

它可以作为掩膜,掩盖多边形范围外的值——即把外面的值设为 NaN,以便进行后续的计算

z[~mask] = np.nan
z_mean = np.nanmean(z)

下图展示了两个应用:左小图的多边形是一个中心带洞的正方形,给定一系列散点的坐标,计算出掩膜后可以把多边形内的散点画成红色,多边形外的散点画成蓝色;右小图的多边形是中国全域,给定 (50, 50) 形状的经纬度网格,计算出掩膜后用橙色画出掩膜为 True 的部分,这张掩膜之后可以用来处理网格上的其它变量。

diagram

本文的目的是介绍如何用 Python 制作掩膜数组,并尽量优化其运行时间。从 shapefile 中读取中国国界并转化为 Shapely 中的多边形对象等操作,已经在博文 Cartopy 系列:探索 shapefile 中详细介绍过了,本文是对其的一个补充。

Cartopy 系列:探索 shapefile

前言

Cartopy 可以通过 feature 模块向地图添加国界 BORDER 和省界 STATES,因其底层采用的 Natural Earth 地图数据并不符合我国的政治主张,所以我们经常需要自备 shapefile 文件来画中国省界,以下面的代码为例

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

extents = [70, 140, 0, 60]
crs = ccrs.PlateCarree()
fig = plt.figure()
ax = fig.add_subplot(111, projection=crs)
ax.set_extent(extents, crs)

filepath = './data/bou2_4/bou2_4p.shp'
reader = shpreader.Reader(filepath)
geoms = reader.geometries()
ax.add_geometries(geoms, crs, lw=0.5, fc='none')
reader.close()

plt.show()

图就不放了,这段代码足以应付大部分需要画省界的情况。然而我在无脑粘贴代码的过程中逐渐产生了疑惑:为什么 shapefile 会由三个文件组成?省界是以何种形式存储在文件中?Cartopy 和 Matplotlib 又是怎样将省界画出来的?调查一番源码后总结出了这段代码底层实现的流程:

  • 利用 PyShp 包读取 shapefile 文件中的每个形状。
  • 利用 Shapely 包将形状转换为几何对象。
  • 利用 Cartopy 包将几何对象投影到地图所在的坐标系上。
  • 用投影后的坐标构造 Matplotlib 的 Path 对象,最后画在地图上。

本文的目的即是从头到尾解说一下这段流程,希望加深对 shapefile 格式,Matplotlib 和 Cartopy 包的理解。令人意外的是,随着探索的不断深入,我发现自己自然而然地学会了如何实现省份填色、省份合并,地图白化等,以前看起来十分困难的操作。本文也会一并介绍这些应用。

Matplotlib 系列:图解 quiver

前言

Matplotlib 中用箭头表示风场或电磁场等矢量场时需要用到 quiver 方法,据字典,quiver 一词的意思是颤动、颤抖或箭袋,貌似也就最后一个意思跟箭头搭得上边。相比于其它画图方法,quiver 的参数又多又容易混淆,所以本文将以图解的方式逐一介绍。这些参数按功能可分为三种:控制箭头位置和数值的、控制箭头长度和角度的,以及控制箭头尺寸和形状的。下面会按照这个分组顺序来解说。本文代码基于 Matplotlib 3.3.4。

Cartopy 系列:对入门教程的补充

前言

几年前曾写过 Cartopy 系列:从入门到放弃,但现在来看还是遗漏了不少细节,比如初学者可能会遇到以下问题

  • 经度是用 [-180°, 180°] 还是 [0°, 360°] 范围?
  • 为什么有时候设置的刻度显示不全?
  • 怎么截取跨越地图边界的区域,画图又怎么跨越边界?

本文将会用坐标变换的思想来解答以上问题,希望能给读者一些实用的启示。本来应该把这些内容写到入门教程里的,但可能会太长,所以现在单独成篇。文中的讨论主要针对最常用的 Plate Carrée 投影,其它投影需要读者自己测试。代码基于 Cartopy 0.18.0,虽然现在已经更新到 0.20.0 了,但基本思想是一致的。

Cartopy 系列:为 Lambert 投影地图添加刻度

前言

Cartopy 中的 Plate Carrée 投影使用方便,但在展示中国地图时会使中国的形状显得很瘪,与之相比,Lambert 投影的效果会更加美观,下图显示了两种投影的差异

comparison

所以本文将会介绍如何在 Cartopy 中实现 Lambert 投影,并为地图添上合适的刻度。文中 Cartopy 的版本是 0.18.0。

Cartopy 系列:从入门到放弃

简介

常用的地图可视化的编程工具有 MATLAB、IDL、GrADS、GMT、NCL 等。我之前一直使用的是脚本语言 NCL,易用性不错,画地图的效果也很好。然而 2019 年初,NCAR 宣布 NCL 将停止更新,并会在日后转为 Python 的绘图包。于是我开始考虑转投 Python,同时觉得在 Python 环境下如果还是用 PyNGL 那一套语法的话,未免有些换汤不换药。因此我选择用 Python 环境下专有的 Cartopy 包来画地图。

cartopy_log

此前 Python 最常用的地图包是 Basemap,然而它将于 2020 年被弃用,官方推荐使用 Cartopy 包作为替代。Cartopy 是英国气象局开发的地图绘图包,实现了 Basemap 的大部分功能,还可以通过 Matplotlib 的 API 实现丰富的自定义效果。

本文将会从一个 NCL 转 Python 的入门者的角度,介绍如何安装 Cartopy,如何绘制地图,并实现一些常用的效果。代码基于 0.18.0 版本的 Cartopy。