WRF后处理/Python处理nc数据与可视化/极坐标网格绘制(Cartopy、netcdf4）——以北极雪水当量数据为例

试了下用python处理并绘制北极雪水当量数据（来源：北极雪水当量格网数据集，，以往数据处理与图像绘制我习惯于使用matlab或R，绘制使用ArcGIS。不过python毕竟是万金油语言，试一试如何处理，将来以备不时之需也好，就当编程复建了。
网上python处理与绘制代码很多，我看了眼，感觉都比较杂乱简洁，且缺少极地投影绘制，对于可能遇到的问题以及整个流程的介绍目前我没找到让我满意的资料，决定现写一篇笔记留着自用。
nc数据算是气象里最常用的数据类型，无论是再分析资料还是模式结果都以这类数据储存，其处理与可视化过程也是基本操作了，希望这篇博文能给初学者一些帮助，毕竟绘图与数据处理是最基本的内容，所涉及的代码也很简单，我这里简单汇总一下。
本文内容包括：
1、netcdf4库对nc数据进行简单读写与处理
2、matplotlib和Cartopy对数据进行可视化
3、matplotlib库绘图细节
4、补充内容：使用numpy库进行经纬度换算
可根据自己需要至不同部分查看。

nc数据读写与处理

首先对对下载好的nc数据进行读取，并查看nc信息，以便后续处理

import numpy as np
import netCDF4 as nc
d=nc.Dataset('E:/arctic/XDA19070302_154/Arctic_SWE_2019.nc')#读取nc数据，2019年泛北极地区雪水当量数据
all_vars = d.variables.keys()   #获取所有变量名称

all_vars_info = d.variables.items()  #获取所有变量信息
print(type(all_vars_info))   #输出为： odict_items 。这里将其转化为 list列表
all_vars_info = list(all_vars_info)
print(all_vars_info)  #显示nc数据信息

1
2
3
4
5
6
7
8
9
10

得到的nc数据信息为：

<class 'dict_items'>
[('time', <class 'netCDF4._netCDF4.Variable'>
float64 time(time)
    standard_name: time
    long_name: time
    units: days since 2019-01-01
    calendar: standard
    axis: T
unlimited dimensions: time
current shape = (365,)
filling off), ('x', <class 'netCDF4._netCDF4.Variable'>
float32 x(x)
    standard_name: projection_x_coordinate
    long_name: x coordinate of projection
    units: m
    axis: X
unlimited dimensions: 
current shape = (978,)
filling off), ('y', <class 'netCDF4._netCDF4.Variable'>
float32 y(y)
    standard_name: projection_y_coordinate
    long_name: y coordinate of projection
    units: m
    axis: Y
unlimited dimensions: 
current shape = (978,)
filling off), ('lambert_azimuthal_equal_area', <class 'netCDF4._netCDF4.Variable'>
int32 lambert_azimuthal_equal_area()
    grid_mapping_name: lambert_azimuthal_equal_area
    false_easting: 0.0
    false_northing: 0.0
    latitude_of_projection_origin: 90.0
    longitude_of_projection_origin: 0.0
    long_name: CRS definition
    longitude_of_prime_meridian: 0.0
    semi_major_axis: 6378137.0
    inverse_flattening: 298.257223563
    spatial_ref: PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",90],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH]]
    GeoTransform: -4889336.08287 10000 0 4890662.63721 0 -10000 
unlimited dimensions: 
current shape = ()
filling off), ('swe', <class 'netCDF4._netCDF4.Variable'>
float64 swe(time, y, x)
    grid_mapping: lambert_azimuthal_equal_area
    _FillValue: -9999.0
    missing_value: -9999.0
unlimited dimensions: time
current shape = (365, 978, 978)
filling off)]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

阅读nc信息时，我们应当注意以下信息：变量、变量维度、地理信息。
在本个nc数据中，我们共有4个变量:time（时间，每日）、x（投影x坐标，经度）y(投影y坐标）、swe(雪水当量，它是时间、x和y共同描述的函数），shape则为变量的维度。
下面，我们开始提取各变量，并进行简单计算，我们需要做的是：根据2019每日雪水当量数据，计算2019年年均雪水当量数据。

time=np.array(d.variables['time'])#时间
d_lon=np.array(d.variables['x'])#经度
d_lat=np.array(d.variables['y'])#维度
swe=np.array(d.variables['swe'])#x雪水当量
FillValue_E=d.variables['swe']._FillValue#雪水当量中存在填充与缺失数据，将其找出
print(FillValue_E)
swe=np.ma.masked_values(swe,FillValue_E)#将填充部分掩盖，方便计算
swe_year=np.transpose(np.sum(swe,axis=0))#将SWE按照时间维度相加，axis=0表示按列求和，求和后将数据转置，否则绘图时经纬度将错乱
swe_year=np.ma.filled(swe_year,FillValue_E#重新将年雪水当量填充，便于输出
1
2
3
4
5
6
7
8
9

好了，我们计算出了2019年的泛北极（45°N-90°N）地区的年雪水当量，接下来我们将其保存，便于之后使用。

f_w = nc.Dataset('E:/arctic/XDA19070302_154/swe_year.nc','w',format = 'NETCDF4')   #创建一个格式为.nc的，名字为 ‘swe_year.nc’的文件 
f_w.createDimension('x',978)   #设置变量x维度
f_w.createDimension('y',978)  #设置变量y维度
f_w.createVariable('x',np.float32,('x'))  #创造变量x，为长度为978，数据类型为float32的数组
f_w.createVariable('y',np.float32,('y'))#同上
f_w.variables['x'][:] = d_lon#指定x值
f_w.variables['y'][:] = d_lat#指定y值
f_w.createVariable( 'swe_year', np.float64, ('x','y'))#创造变量swe_year，由变量x和有描述

f_w.variables['swe_year'][:] = swe_year#赋值
f_w.close()
1
2
3
4
5
6
7
8
9
10
11

以上便完成了nc数据的基本操作。

二、数据可视化(Cartopy与matplotlib）

cartopy与Matplotlib关系
关于Cartopy的基本绘图流程，网络上有许多资料，在这里我推荐一篇：Cartopy入门到放弃
简单来讲，Cartipy绘制地图的底图，而Matplotlib则负责将数据画在地图上。
在绘制图形时，你需要知道：画布大小、地图投影方式、地图显示范围、数据投影方式（重要）
如果投影设置错误，或者经纬度范围设置出现问题，会导致你的图只有地图的底图，也会影响绘图美观性（血泪教训……）
地图底图绘制
先加载个包，进行一些基本设置

import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.path as mpath
import cmaps
mpl.rcParams["font.family"] = 'Arial'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 16   #字体大小
mpl.rcParams["axes.linewidth"] = 1   #轴线边框粗细
1
2
3
4
5
6
7
8
9
10
11
12
13
14

开始绘制底图：

proj =ccrs.NorthPolarStereo(central_longitude=0)#设置地图投影
#在圆柱投影中proj = ccrs.PlateCarree(central_longitude=xx)
leftlon, rightlon, lowerlat, upperlat = (-180,180,45,90)#经纬度范围

img_extent = [leftlon, rightlon, lowerlat, upperlat]
fig1 = plt.figure(figsize=(8,8))#设置画布大小
f1_ax1 = fig1.add_axes([0.2, 0.3, 0.5, 0.5],projection = ccrs.NorthPolarStereo())#绘制地图位置
#注意此处添加了projection = ccrs.NorthPolarStereo()，指明该axes为北半球极地投影
#f1_ax1.gridlines()
f1_ax1.set_extent(img_extent, ccrs.PlateCarree())
f1_ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
#######以下为网格线的参数######
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
##############################
f1_ax1.set_boundary(circle, transform=f1_ax1.transAxes)
#设置axes边界，为圆形边界，否则为正方形的极地投影
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

此时，你绘制的图是这样的：

这是地图的底图，我们还没有将数据画上去。
年均雪水当量绘制

swe=nc.Dataset('E:/arctic/XDA19070302_154/swe_year.nc')
lon=np.array(swe.variables['x'])
lat=np.array(swe.variables['y'])
swe_year=np.array(swe.variables['swe_year'])
#导入数据
d=np.ma.masked_values(swe_year,-9999.0)
d=d/365#年均雪水当量
data_proj =ccrs.LambertAzimuthalEqualArea(central_latitude=90, central_longitude=0)#设置数据的投影信息，注意查看原始nc文件夹中的投影信息，此处为LambertAzimuthalEqualArea，中心经纬度为90°N，0°
c7=f1_ax1.pcolormesh(lon,lat,d,cmap=cmaps.BlAqGrYeOrRe,transform=data_proj,vmin=0,vmax=350)#绘制年均雪水当量
position=fig1.add_axes([0.2, 0.25, 0.55, 0.025])#图标位置

font = {'family' : 'serif',
        'color'  : 'darkred',
        'weight' : 'normal',
        'size'   : 16,
        }
cb=fig1.colorbar(c7,cax=position,orientation='horizontal',format='%.1f',extend='both')#设置图标
cb.set_label('SWE(mm)',fontdict=font) #添加图标标签
#fig1.savefig('E:/arctic/2019SWE.jpg')#保存
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

最后，出的图：

2019年泛北极地区年均雪水当量
对比一下数据文档给出的：

集成的泛北极地区年均雪水当量（2002-2017年）
大致的分布还是正确的，具体的绘图细节可自行设置。
色标
python自带的色标并不美观，气象家园有大神将NCL的色标移至了python中，库为cmaps，我在绘图中使用的便是这个库的色标，详见cmaps
坐标添加
Cartopy本身在绘制极地投影时一个bug，由于自带的边界是方形的，在绘制时，我们要给它绘制上圆形边界，此时经纬度变无法添加，只能自己根据经纬度大致位置，当作文本添加至图中，添加坐标可见：
How to Add More Lines of Longitude/latitude in NorthPolarStereo Projection
在这里我没有添加，偷个懒……

三、 经纬度换算(XY坐标换至经纬度坐标）

这个nc数据的经纬度用XY坐标进行描述的，虽然对于图形绘制没有影响，但考虑到经纬度坐标对于经纬度的设置更加直观方便，这里贴上我将此数据坐标换算的代码：

import numpy as np
import math
'''
def XYtoGPS(x, y, ref_lat, ref_lon):
        CONSTANTS_RADIUS_OF_EARTH = [6371000 for i in range(978)]   # meters (m)
        CONSTANTS_RADIUS_OF_EARTH=np.array(CONSTANTS_RADIUS_OF_EARTH,dtype=np.float32)
        x_rad = np.divide(x,CONSTANTS_RADIUS_OF_EARTH)
        y_rad = np.divide(x,CONSTANTS_RADIUS_OF_EARTH)
        c = np.sqrt(x_rad * x_rad + y_rad * y_rad)

        ref_lat_rad = np.radians(ref_lat)
        ref_lon_rad = np.radians(ref_lon)

        ref_sin_lat = np.sin(ref_lat_rad)
        ref_cos_lat = np.cos(ref_lat_rad)

        if abs(c) > 0:
            sin_c = np.sin(c)
            cos_c = np.cos(c)

            lat_rad = np.asin(cos_c * ref_sin_lat + (x_rad * sin_c * ref_cos_lat) / c)
            lon_rad = (ref_lon_rad + np.atan2(y_rad * sin_c, c * ref_cos_lat * cos_c - x_rad * ref_sin_lat * sin_c))

            lat = np.degrees(lat_rad)
            lon = np.degrees(lon_rad)

        else:
            lat = np.degrees(ref_lat)
            lon = np.degrees(ref_lon)

        return lat, lon
'''
ref_lat=[90 for i in range (978)]
ref_lon=[0 for i in range (978)]
CONSTANTS_RADIUS_OF_EARTH = [6371000 for i in range(978)]   # meters (m)
CONSTANTS_RADIUS_OF_EARTH=np.array(CONSTANTS_RADIUS_OF_EARTH,dtype=np.float32)
x=lon
y=lat
x_rad =np.divide(x,CONSTANTS_RADIUS_OF_EARTH)
y_rad = np.divide(x,CONSTANTS_RADIUS_OF_EARTH)
c = np.sqrt(x_rad * x_rad + y_rad * y_rad)
ref_lat_rad = np.radians(ref_lat)
ref_lon_rad = np.radians(ref_lon)
lat1=[0 for i in range (978)]
lon1=[0 for i in range (978)]
ref_sin_lat = np.sin(ref_lat_rad)
ref_cos_lat = np.cos(ref_lat_rad)
A=np.arange(0,977,1)
for i in A:
    if abs(c[i])>0:
        sin_c = math.sin(c[i])
        cos_c = math.cos(c[i])
        lat_rad = math.asin(cos_c * ref_sin_lat[i] + (x_rad[i] * sin_c * ref_cos_lat[i]) / c[i])
        lon_rad = (ref_lon_rad[i] + math.atan2(y_rad[i] * sin_c, c[i] * ref_cos_lat[i] * cos_c - x_rad[i] * ref_sin_lat[i] * sin_c))
        lat1[i] = math.degrees(lat_rad)
        lon1[i]= math.degrees(lon_rad)
    else:
        lat1[i] = math.degrees(ref_lat[i])
        lon1[i]= math.degrees(ref_lon[i])
        
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

四、总结

数据处理方面
python对于nc数据的处理，我个人觉得无功无过，相对于Matlab和R而言，没有多方便，也没有多繁琐，数据类型也类似，实际上，对于任何的数据处理，这些语言都是大同小异，只要掌握了常见的数据类型：列表、数组、字典等，任何语言处理数据都什么类似，当然，不同的语言会有一些特色数据类型，比如matlab的cell，R的数据框，不过实际操作起来差别不大，熟悉哪个用哪个吧。
绘图方面
我的评价是：不如ArcGIS(×）
图标设置和画图大小、相对位置设置非常繁琐，配色也不大好看，绘制非等间距配色很麻烦，远不如ArcGIS直观简洁。
不过ArcGIS是专业的地理信息软件，也许不存在可比性。
如果需要ArcGIS绘图，可以用gdal库将nc转为tif，再在GIS中绘图，但是导出为tif文件matlav和R也能做……
综上，python处理nc与进行WRF后处理的优势是：万金油语言，可以同时实现下载数据、处理数据、绘图，比较方便。
但是，无论是数据处理还是绘图都没有什么优势，如果像我一样习惯用matlab和R处理数据，用ArcGIS出图的话，用python大可不必，反而会因为不熟悉debug很久。
不过闲的没事跑一跑也没坏处，我当作简单的编程复建还是有点用的（×）
如果之后有空可能会尝试用GDAL处理NC，再用GIS出图来对比一下。
最后，附上完整代码：

import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.path as mpath
import cmaps
d=nc.Dataset('E:/arctic/XDA19070302_154/Arctic_SWE_2019.nc')
all_vars = d.variables.keys()   #获取所有变量名称

all_vars_info = d.variables.items()  #获取所有变量信息
print(type(all_vars_info))   #输出为： odict_items 。这里将其转化为 list列表
all_vars_info = list(all_vars_info)
print(all_vars_info)  
d_lon=np.array(d.variables['x'])
d_lat=np.array(d.variables['y'])
time=np.array(d.variables['time'])
d_lon=np.array(d.variables['x'])
d_lat=np.array(d.variables['y'])
swe=np.array(d.variables['swe'])
print("E's _FillValue are:")
FillValue_E=d.variables['swe']._FillValue
print(FillValue_E)
swe=np.ma.masked_values(swe,FillValue_E)
swe_year=np.transpose(np.sum(swe,axis=0))
swe_year=np.ma.filled(swe_year,FillValue_E)

f_w = nc.Dataset('E:/arctic/XDA19070302_154/swe_year.nc','w',format = 'NETCDF4')   #创建一个格式为.nc的，名字为 ‘hecheng.nc’的文件 
f_w.createDimension('x',978)   
f_w.createDimension('y',978)  
f_w.createVariable('x',np.float32,('x'))  
f_w.createVariable('y',np.float32,('y'))
f_w.variables['x'][:] = d_lon
f_w.variables['y'][:] = d_lat
f_w.createVariable( 'swe_year', np.float64, ('x','y'))

f_w.variables['swe_year'][:] = swe_year
f_w.close()

mpl.rcParams["font.family"] = 'Arial'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 16   #字体大小
mpl.rcParams["axes.linewidth"] = 1   #轴线边框粗细（默认的太粗了）

swe=nc.Dataset('E:/arctic/XDA19070302_154/swe_year.nc')
lon=np.array(swe.variables['x'])
lat=np.array(swe.variables['y'])
swe_year=np.array(swe.variables['swe_year'])

proj =ccrs.NorthPolarStereo(central_longitude=0)
#在圆柱投影中proj = ccrs.PlateCarree(central_longitude=xx)
leftlon, rightlon, lowerlat, upperlat = (-180,180,45,90)
#仅画60°E-90°E部分
img_extent = [leftlon, rightlon, lowerlat, upperlat]
fig1 = plt.figure(figsize=(8,8))
f1_ax1 = fig1.add_axes([0.2, 0.3, 0.5, 0.5],projection = ccrs.NorthPolarStereo())
#注意此处添加了projection = ccrs.NorthPolarStereo()，指明该axes为北半球极地投影
#f1_ax1.gridlines()
f1_ax1.set_extent(img_extent, ccrs.PlateCarree())
f1_ax1.add_feature(cfeature.COASTLINE.with_scale('110m'))
#######以下为网格线的参数######
theta = np.linspace(0, 2*np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)
##############################
f1_ax1.set_boundary(circle, transform=f1_ax1.transAxes)
#设置axes边界，为圆形边界，否则为正方形的极地投影

d=np.ma.masked_values(swe_year,-9999.0)
d=d/365
#fig1.text(x, y, r'0$^\circ$',fontsize=14, horizontalalignment='center',verticalalignment='center')

data_proj =ccrs.LambertAzimuthalEqualArea(central_latitude=90, central_longitude=0)
c7=f1_ax1.pcolormesh(lon,lat,d,cmap=cmaps.BlAqGrYeOrRe,transform=data_proj,vmin=0,vmax=350)

#f1_ax1.contourf(lon,lat,swe_year, zorder=0, levels =np.arange(-0.6,0.7,0.1) , extend = 'both',transform=data_proj, cmap=plt.cm.RdBu_r)
position=fig1.add_axes([0.2, 0.25, 0.55, 0.025])

font = {'family' : 'serif',
        'color'  : 'darkred',
        'weight' : 'normal',
        'size'   : 16,
        }
cb=fig1.colorbar(c7,cax=position,orientation='horizontal',format='%.1f',extend='both')
cb.set_label('SWE(mm)',fontdict=font) 
#fig1.savefig('E:/arctic/2019SWE.jpg')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92