试了下用python处理并绘制北极雪水当量数据（来源：北极雪水当量格网数据集，，以往数据处理与图像绘制我习惯于使用matlab或R，绘制使用ArcGIS。不过python毕竟是万金油语言，试一试如何处理，将来以备不时之需也好，就当编程复建了。
网上python处理与绘制代码很多，我看了眼，感觉都比较杂乱简洁，且缺少极地投影绘制，对于可能遇到的问题以及整个流程的介绍目前我没找到让我满意的资料，决定现写一篇笔记留着自用。
nc数据算是气象里最常用的数据类型，无论是再分析资料还是模式结果都以这类数据储存，其处理与可视化过程也是基本操作了，希望这篇博文能给初学者一些帮助，毕竟绘图与数据处理是最基本的内容，所涉及的代码也很简单，我这里简单汇总一下。
本文内容包括：
1、netcdf4库对nc数据进行简单读写与处理
2、matplotlib和Cartopy对数据进行可视化
3、matplotlib库绘图细节
4、补充内容：使用numpy库进行经纬度换算
可根据自己需要至不同部分查看。
nc数据读写与处理
首先对对下载好的nc数据进行读取，并查看nc信息，以便后续处理
import numpy as npimport netCDF4 as ncd=nc.Dataset('E:/arctic/XDA19070302_154/Arctic_SWE_.nc')#读取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) #显示nc数据信息
得到的nc数据信息为：
<class 'dict_items'>[('time', <class 'netCDF4._netCDF4.Variable'>float64 time(time)standard_name: timelong_name: timeunits: days since -01-01calendar: standardaxis: Tunlimited dimensions: timecurrent shape = (365,)filling off), ('x', <class 'netCDF4._netCDF4.Variable'>float32 x(x)standard_name: projection_x_coordinatelong_name: x coordinate of projectionunits: maxis: Xunlimited dimensions: current shape = (978,)filling off), ('y', <class 'netCDF4._netCDF4.Variable'>float32 y(y)standard_name: projection_y_coordinatelong_name: y coordinate of projectionunits: maxis: Yunlimited 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_areafalse_easting: 0.0false_northing: 0.0latitude_of_projection_origin: 90.0longitude_of_projection_origin: 0.0long_name: CRS definitionlongitude_of_prime_meridian: 0.0semi_major_axis: 6378137.0inverse_flattening: 298.257223563spatial_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.0missing_value: -9999.0unlimited dimensions: timecurrent shape = (365, 978, 978)filling off)]
阅读nc信息时，我们应当注意以下信息：变量、变量维度、地理信息。
在本个nc数据中，我们共有4个变量:time（时间，每日）、x（投影x坐标，经度）y(投影y坐标）、swe(雪水当量，它是时间、x和y共同描述的函数），shape则为变量的维度。
下面，我们开始提取各变量，并进行简单计算，我们需要做的是：根据每日雪水当量数据，计算年均雪水当量数据。
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#重新将年雪水当量填充，便于输出
好了，我们计算出了的泛北极（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 npimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport matplotlib.ticker as mtickerimport netCDF4 as ncimport matplotlib.path as mpathimport cmapsmpl.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.5verts = np.vstack([np.sin(theta), np.cos(theta)]).Tcircle = 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/SWE.jpg')#保存
最后，出的图：
泛北极地区年均雪水当量
对比一下数据文档给出的：
集成的泛北极地区年均雪水当量（2002-）
大致的分布还是正确的，具体的绘图细节可自行设置。
色标
python自带的色标并不美观，气象家园有大神将NCL的色标移至了python中，库为cmaps，我在绘图中使用的便是这个库的色标，详见cmaps
坐标添加
Cartopy本身在绘制极地投影时一个bug，由于自带的边界是方形的，在绘制时，我们要给它绘制上圆形边界，此时经纬度变无法添加，只能自己根据经纬度大致位置，当作文本添加至图中，添加坐标可见：
How to Add More Lines of Longitude/latitude in NorthPolarStereo Projection
在这里我没有添加，偷个懒……
三、 经纬度换算(XY坐标换至经纬度坐标）
这个nc数据的经纬度用XY坐标进行描述的，虽然对于图形绘制没有影响，但考虑到经纬度坐标对于经纬度的设置更加直观方便，这里贴上我将此数据坐标换算的代码：
import numpy as npimport 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=lony=latx_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 npimport netCDF4 as ncimport matplotlib.pyplot as plt###引入库包import numpy as npimport matplotlib as mplimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTERimport matplotlib.ticker as mtickerimport netCDF4 as ncimport matplotlib.path as mpathimport cmapsd=nc.Dataset('E:/arctic/XDA19070302_154/Arctic_SWE_.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']._FillValueprint(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_lonf_w.variables['y'][:] = d_latf_w.createVariable( 'swe_year', np.float64, ('x','y'))f_w.variables['swe_year'][:] = swe_yearf_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.5verts = np.vstack([np.sin(theta), np.cos(theta)]).Tcircle = 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/SWE.jpg')
WRF后处理/Python处理nc数据与可视化/极坐标网格绘制(Cartopy netcdf4）——以北极雪水当量数据为例
