资讯详情

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

试了下用python处理和绘制北极雪水当量数据(来源:北极雪水当量格网数据集,我习惯于使用过去的数据处理和图像绘制。matlab或R,绘制使用ArcGIS。不过python毕竟是万金油语言,试试怎么处理,以后以备不时之需,就当编程复建了。 网上python有很多处理和绘制代码。当我看到它们时,我感到凌乱和简单,缺乏极地投影和绘制。目前,我还没有找到令我满意的信息来介绍可能遇到的问题和整个过程,所以我决定写一个笔记来自己使用。 nc数据是气象中最常用的数据类型,无论是再分析数据还是模式结果都存储在这些数据中,其处理和可视化过程也是基本操作,希望博客能给初学者一些帮助,毕竟,绘图和数据处理是最基本的内容,代码也很简单,我在这里简要总结。 本文包括: 1、 2、 3、 4、 可根据自己的需要到不同的部分查看。

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的基本绘图流程,网络上有许多资料,在这里我推荐一篇: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')#保存

最后,出的图: 对比一下数据文档给出的: 大致的分布还是正确的,具体的绘图细节可自行设置。 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
        标签: 色标传感器wm03nct2a色标光电传感器红外线电梯光幕lmp331液压变送器

锐单商城拥有海量元器件数据手册IC替代型号,打造 电子元器件IC百科大全!

锐单商城 - 一站式电子元器件采购平台