利用 Python 地理制图与 NetCDF 数据处理
本文介绍了如何使用 Python 的 cartopy 库来绘制 CGCS2000 坐标系下的中国地图，以及如何使用 netCDF4 库来读取 NetCDF 数据并将其绘制在地图上。
预备知识
NetCDF 简介
NetCDF（Network Common Data Format，网络通用数据格式）是由美国大学大气研究中心（University Corporation for Atmospheric Research，UCAR）的 Unidata 项目支持的一种面向阵列的数据格式，历史可以追溯到 1988 年。UCAR 也是 netCDF 软件、标准开发和更新维护的主力军。NetCDF 格式标准是一个开放的标准，并被开放地理空间联盟（Open Geospatial Consortium）收录。
NetCDF 旨在提供一种简单的、可移植的、自描述的方法来存储科学数据。NetCDF 由三个部分组成：数据模型、API 和文件格式。NetCDF 数据模型是一种面向阵列的数据模型，它可以表示任意维度的数据。NetCDF API 是一组用于创建、操作和访问 NetCDF 数据的函数。NetCDF 文件格式是一种用于存储 NetCDF 数据的文件格式，它是一种二进制文件格式，它可以存储任意维度的数据。
CGCS2000 坐标系
现行《中华人民共和国测绘法》第十条规定，“国家建立全国统一的大地坐标系统、平面坐标系统、高程系统、地心坐标系统和重力测量系统，确定国家大地测量等级和精度以及国家基本比例尺地图的系列和基本精度”。1
建国以来，我国于上世纪 50 年代和 80 年代分别建立了 1954 年北京坐标系（俗称北京 54）和 1980 年西安坐标系（俗称西安 80），并基于此两种坐标系测制了各种比例尺地形图。但是限于当时的技术条件，我国大地坐标系基本上是依赖于传统技术手段实现的。北京 54 坐标系采用的是 Красовский 椭球体（一般称作 SK-42 椭球体）。该椭球在计算和定位的过程中，没有采用中国的数据，该系统在我国范围内符合得不好。后来到了 70 年代，经过整体平差，采用 1975 年 IUGG（International Union of Geodesy and Geophysics，国际大地测量学与地球物理学联合会）第十六届大会推荐的参考椭球参数，我国建立了 1980 西安坐标系。但是由于近年来的发展，西安 80 坐标系存在以下不足：2
-
二维坐标系统。1980 西安坐标系是经典大地测量成果的归算及其应用，它的表现形式为平面的二维坐标。用现行坐标系只能提供点位平面坐标，而且表示两点之间的距离精确度也比用现代手段测得的低10倍左右。高精度、三维与低精度、二维之间的矛盾是无法协调的。比如将卫星导航技术获得的高精度的点的三维坐标表示在现有地图上，不仅会造成点位信息的损失（三维空间信息只表示为二维平面位置），同时也将造成精度上的损失。
-
参考椭球参数。随着科学技术的发展，国际上对参考椭球的参数已进行了多次更新和改善。1980 西安坐标系所采用的 IAG 1975 椭球，其长半轴要比现在国际公认的 WGS84 椭球长半轴的值大 3 米左右，而这可能引起地表长度误差达 10 倍左右。
在此基础上，为了修正西安 80 坐标系带来的偏差，国家出台了 CGCS 2000 坐标系。
Plate carrée 投影法 和 WGS84 地理坐标系
等角投影（也称为等距圆柱投影或 la carte parallélogramatique 投影），包括 plate carrée 投影的特殊情况，是一种简单的地图投影，由 Μαρῖνος ὁ Τύριος（泰尔的马里努斯，古希腊地理学家）发明。据 Πτολεμαῖος（托勒密，古希腊数学家）称 Μαρῖνος 在公元 100 年左右发明了这种投影。
该投影将经线映射为间距恒定的垂直直线（用于恒定间距的经线间隔），将纬线圈映射为间距恒定的水平直线（用于恒定间距的平行线）。这就导致该投影既不是等面积的，也不是保形的。由于这种投影带来的扭曲，它在导航或地籍测绘中实际上没有什么用处。但是由于点在地图上的相应地理位置之间的关系特别简单，因此它在地理信息系统中有着重要的作用。
世界大地测量系统（World Geodetic System, WGS）是一个用于制图学、大地测量学和卫星导航（包括GPS）的标准。目前的版本是 WGS84，定义了一个以地球为中心、地球固定的坐标系统和大地测量基准，还描述了相关的地球引力模型（Earth Gravitational Model, EGM）和世界磁场模型（World Magnetic Model, WMM）。该标准由美国国家地理空间情报局出版和维护。它是目前世界上最常用的地球椭球体标准。
有关制图的国家标准
国家标准 GB/T 33181-2016《国家基本比例尺地图 1:250 000 1:500 000 1:1 000 000 地形图》规定，
- 1:250 000、1:500 000 地图采用 Gauss-Krüger 投影，从东经 0° 起算，按经差 6° 分带。各分带以中央子午线为直角坐标纵轴，向西平移 500 km，以赤道为直角坐标系横轴。
- 1:1 000 000 地形图各图幅单独采用双标准纬线正轴等角圆锥投影，按纬差 4° 分带。每幅图具有两条标准纬线，其纬度为：（ 和 分别代表图幅的南北边纬度）。
Python 部分
利用 cartopy 绘制 CGCS2000 坐标系下的中国地图
cartopy 是一个 Python 的地图绘制库，它提供了一系列的地图投影和坐标系，可以用来绘制地图。cartopy 依赖于 matplotlib，因此在使用 cartopy 之前，需要先安装 matplotlib。cartopy 的安装可以通过 pip 安装，也可以通过 conda 安装。在本文中，我们使用 pip 安装 cartopy。
# Windows
$ pip install cartopy
# MacOS, Linux
$ pip3 install cartopy在本文中，我们使用两个现成的 shapefile（一个是中国省界、一个是九段线）和一个自带的 shapefile（cartopy 在 Amazon 上存储好的）来绘制正轴等角圆锥投影下的中国地图。
首先我们来导入 shapefile 矢量地图。值得注意的是，我这里使用的 shapefile 是 WGS84 坐标系下的，因此在绘制地图时，需要指定地图的坐标系和数据的坐标系。
# demo.py
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
# 省界
province_line = 'province/province.shp'
province_feature = ShapelyFeature(
Reader(province_line).geometries(),
crs=ccrs.PlateCarree(), # 地图原始坐标系
facecolor='none' # 不填充颜色
)
# 九段线
jiuduanxian = 'jiuduanxian/jiuduanxian.shp'
jiuduanxian_feature = ShapelyFeature(
Reader(jiuduanxian).geometries(),
crs=ccrs.PlateCarree(), # 地图原始坐标系
facecolor='none' # 不填充颜色
)接下来，要根据 CGCS2000 坐标系的参数定义我们要使用的投影坐标系。
根据 epsg.io 的介绍，CGCS2000 使用的椭球体是 GRS80 椭球体，因此我们需要先定义一个 GRS80 椭球体。
globe = ccrs.Globe(ellipse='GRS80')接下来是定义一个正轴等角圆锥投影的坐标系。对于完整的中国地图，正轴等角圆锥投影（Lambert 等角圆锥投影，保留角度）的参数一般是：
- 中央经线：105°E
- 中央纬线：36°N
- 标准纬线：25°N、47°N
我们使用 cartopy.crs 中的 LambertConformal 类来定义一个正轴等角圆锥投影的坐标系：
lcc = ccrs.LambertConformal(
central_longitude=105,
central_latitude=36,
standard_parallels=(25, 47),
globe=globe
)剩下就非常简单了：
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=lcc)
ax.add_feature(province_feature, edgecolor='black', linewidth=0.5) # 绘制省界
ax.add_feature(jiuduanxian_feature, edgecolor='black', linewidth=0.5) # 绘制九段线
ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.5) # 添加海岸线
ax.set_extent([73, 135, 15, 55], crs=ccrs.PlateCarree()) # 设置地图范围
ax.gridlines(draw_labels=True, dms=False, x_inline=False, y_inline=False) # 绘制经纬线
plt.savefig('china.png', dpi=180)结果如下：
Python 读取 NetCDF 数据
接下来就是如何读取 NetCDF 数据，并把它绘制在地图上。在本文中，我们使用 netCDF4 库来读取 NetCDF 数据，安装：
$ pip install netCDF4这里我们手头有一个 EDGARv6.0 的 2018 年全球甲烷排放数据3，精度为 0.1°x0.1°，文件名为默认的 v6.0_CH4_2018_TOTALS.0.1x0.1.nc。但是，在写代码之前我们得弄明白这个数据的结构才行。我们可以使用 ncdump 命令来查看 NetCDF 文件的结构：
$ ncdump -h v6.0_CH4_2018_TOTALS.0.1x0.1.nc
netcdf v6.0_CH4_2018_TOTALS.0.1x0.1 {
dimensions:
lat = 1800 ;
lon = 3600 ;
variables:
float lat(lat) ;
lat:standard_name = "latitude" ;
lat:long_name = "latitude" ;
lat:units = "degrees_north" ;
lat:comment = "center_of_cell" ;
float lon(lon) ;
lon:standard_name = "longitude" ;
lon:long_name = "longitude" ;
lon:units = "degrees_east" ;
lon:comment = "center_of_cell" ;
float emi_ch4(lat, lon) ;
emi_ch4:standard_name = "tendency_of_atmosphere_mass_content_of_methane_due_to_emission" ;
emi_ch4:long_name = "Emissions of CH4 - " ;
emi_ch4:units = "kg m-2 s-1" ;
emi_ch4:cell_method = "time: mean (interval: 1 year, 365 days)" ;
emi_ch4:total_emi_ch4 = " 3.75268e+011 kg/year" ;
emi_ch4:comment = " (see http://edgar.jrc.ec.europa.eu/methodology.php#12sou for the definitions of the single sources)" ;
// global attributes:
...
}global attributes 被我省略了，因为没什么用。我们可以看到，这个 NetCDF 文件有两个维度，分别是 lat 和 lon，分别对应纬度和经度。这个 NetCDF 文件中有一个变量，名为 emi_ch4，它的维度是 lat 和 lon，也就是说，这个变量是一个二维数组，每个元素都是一个浮点数。这个变量的单位是 kg m-2 s-1，也就是说，这个变量的每个元素都是一个甲烷排放量，单位是每平方米每秒。这个变量的 cell_method 属性是 time: mean (interval: 1 year, 365 days)，也就是说，这个变量是一个年平均值。这个变量的 total_emi_ch4 属性是 3.75268e+011 kg/year，也就是说，这个变量的总排放量是 375.268 Tg。
了解到了这些后，我们来进行作图：
import netCDF4
import cartopy.crs as ccrs
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import matplotlib.pyplot as plt
import numpy as np
data = netCDF4.Dataset('v6.0_CH4_2018_TOTALS.0.1x0.1.nc').variables['emi_ch4'][:]
lons = netCDF4.Dataset('v6.0_CH4_2018_TOTALS.0.1x0.1.nc').variables['lon'][:]
lats = netCDF4.Dataset('v6.0_CH4_2018_TOTALS.0.1x0.1.nc').variables['lat'][:]
point_lons = np.tile(lons, (lats.size, 1))
point_lats = np.tile(lats, (lons.size, 1)).T
globe = ccrs.Globe(ellipse='GRS80')
lcc = ccrs.LambertConformal(
central_longitude=105,
central_latitude=36,
standard_parallels=(25, 47),
globe=globe
)
fname = 'province/province.shp'
shape_feature = ShapelyFeature(
Reader(fname).geometries(),
crs=ccrs.PlateCarree(),
facecolor='none'
)
plt.figure(figsize=(8, 7))
ax = plt.axes(projection=lcc)
ax.add_feature(shape_feature, edgecolor='black', linewidth=0.25)
colormap = plt.cm.get_cmap('BuPu')
vmin = 0; vmax = 1e-9
sm = plt.cm.ScalarMappable(cmap=colormap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
ax.pcolormesh(lcc.transform_points(ccrs.PlateCarree(), point_lons, point_lats)[:, :, 0], lcc.transform_points(ccrs.PlateCarree(), point_lons, point_lats)[:, :, 1], data, cmap=colormap, vmin=vmin, vmax=vmax)
ax.coastlines(resolution='50m', color='black', linewidth=0.25)
ax.set_extent([79, 127, 15, 55], crs=ccrs.PlateCarree())
cbar = plt.colorbar(sm, shrink=0.78, pad=0.05, extend='max')
cbar.set_label('$\\mathrm{CH_4\,emission\,[kg/(m^2·s)]}$', fontsize=10)
plt.savefig('EDGARv6.png', dpi=180, bbox_inches='tight')结果如下：
FAQ
一般来说，cartopy 的安装需要设备上已经安装了 GDAL 和 PROJ 库。如果你的设备上没有安装这两个库，那么你需要先安装这两个库，然后再安装 cartopy。对于 Windows 系统，非常推荐使用 conda 来安装 GDAL 和 PROJ 库，因为 conda 会自动处理好这两个库的依赖关系。对于 MacOS 和 Linux 系统，用包管理器安装即可。以 MacOS 为例，可以使用 brew 来安装 GDAL 和 PROJ 库：
brew install gdal && brew install proj欢迎联系我获取本文涉及的数据！
参考资料
参考资料 & 附注
-
全国人大常委会. 《中华人民共和国测绘法》：中华人民共和国主席令第 67 号[A]. ↩
-
中华人民共和国中央人民政府. 我国自 2008 年 7 月 1 日起将启用 2000 国家大地坐标系[E/OL] (2008-06-27) [2022-12-31]. http://www.gov.cn/gzdt/2008-06/27/content_1029576.htm ↩
-
Joint Research Centre Data Catalogue, European Commission. EDGAR v6.0 Greenhouse Gas Emissions[E/OL] (2021-05-17) [2023-01-02]. https://data.jrc.ec.europa.eu/dataset/97a67d67-c62e-4826-b873-9d972c4f670b ↩