WRF后处理/Python处理nc数据与可视化/极坐标网格绘制(Cartopy、netcdf4)——以北极雪水当量数据为例
时间:2022-09-28 18:30:00
试了下用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数据信息
得到的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)]
阅读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#重新将年雪水当量填充,便于输出
好了,我们计算出了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()
以上便完成了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 #轴线边框粗细
开始绘制底图:
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边界,为圆形边界,否则为正方形的极地投影
此时,你绘制的图是这样的:
这是地图的底图,我们还没有将数据画上去。
年均雪水当量绘制
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')#保存
最后,出的图:
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])
四、总结
数据处理方面
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元器件数据手册
、IC替代型号,打造电子元器件IC百科大全!