怎么使用Python对NetCDF数据做空间相关分析

这篇文章主要介绍了怎么使用python对NetCDF数据做空间相关分析，具有一定借鉴价值，感兴趣的朋友可以参考下，希望大家阅读完这篇文章之后大有收获，下面让小编带着大家一起了解一下。

Python有哪些常用库

python常用的库：1.requesuts；2.scrapy；3.pillow；4.twisted；5.numpy；6.matplotlib；7.pygama；8.ipyhton等。

引言：我一直想理解空间相关分析的计算思维，于是今天又拿起Python脚本和数据来做练习。首先需要说明的是，这次实验的数据和Python脚本均来自于[好久不见]大佬，在跟大佬说明之后，允许我写到公众号来与大家共享，在此对大佬的指点表示感谢，这次实验的脚本可在气象家园或简书app(如果没记错的话)搜索到这次实验的相关内容，也可以微信或者后台发消息给我获取。在此之前我觉得自己还没理解这个方法的计算思维，检验的标准就是我能否迅速运用到其他方面。于是今天又重新回来温习一遍，我把自己的理解与大伙共同交流。

首先，数据的格式是NetCDF(.nc)数据，两个数据分别是[哈德来中心海温sst数据，pc数据是对东太平洋SSTA做的EOF获取]。知道数据信息之后我们就准备开始去运行程序。原始脚本包括了回归分析和相关分析两部分，但是今天我做了空间相关分析这一部分，有兴趣的可以到[好久不见]大佬的气象家园阅读喔！如果还没有安装Cartopy包的话请在后台联系我喔

为了方便理解每一步，我选择去Jupyter运行，因为可以一段一段程序的运行，这是比较方便的。绘图部分并不是很难，关键还是在于数据预处理部分。

空间相关分析的脚本如下：

import numpy as np #数值计算用，如相关系数import xarray as xr #读取.nc文件用from sklearn.feature_selection import f_regression #做显著性检验import matplotlib.pyplot as plt #绘制和展示图形用import cartopy.crs as ccrs #绘制地图用，如果没有安装好的话，请在后台联系我import cartopy.feature as cfeature #添加一些矢量用，这里没用到，因为我没数据from cartopy.mpl.ticker import LongitudeFORMatter, LatitudeFormatter #经纬度格式设置import cmaps #ncl的color，如果没有的话，请联系我，也可以在气象家园找到#使用上下文管理器读取.nc数据，并提取数据中的变量，可以提前用NASA的panoply这个软件查看.nc信息with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1:      pre = f1['sst_anom'][:-1, :, :]  # 三维数据全取，时间，纬度+经度      lat, lon = f1['lat'], f1['lon'] #提取经纬度，后面格网化需要用到pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])#0表示行个数，1列代表的个数，2经度代表个数with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2:      pc = f2['pc'][0, :]# 相关系数计算pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))# 做显著性检验pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat), len(lon))#用0代替NaNarea = np.where(pre_cor_sig < 0.05)# numpy的作用又来了 nx, ny = np.meshgrid(lon, lat)  # 格网化经纬度，打印出来看看就知道为什么要这么做了plt.figure(figsize=(16, 8)) #创建一个空画布#让colorbar字体设置为新罗马字符plt.rcParams['font.family'] = 'Times New Roman'plt.rcParams['font.size'] = 16ax2 = plt.subplot(projection=ccrs.PlateCarree(central_longitude=180))# 在画布上绘图，这个叫axes，这不是坐标轴喔ax2.coastlines(lw=0.4)ax2.set_global()c2 = ax2.contourf(nx, ny, pre_cor, extend='both', cmap=cmaps.nrl_sirkes, transform=ccrs.PlateCarree())plt.colorbar(c2,fraction=0.05,orientation='horizontal', shrink=0.4, pad=0.06)# extend关键字设置colorbar的形状，both为两端尖的，pad是距离主图的距离，其他参数WEB搜索# 显著性打点sig2 = ax2.scatter(nx[area], ny[area], marker='+', s=1, c='k', alpha=0.6, transform=ccrs.PlateCarree())# 凸显显著性区域plt.title('Correlation Analysis', fontdict={'family' : 'Times New Roman', 'size'   : 16})#标题字体也修改为新罗马字符，数字和因为建议都用新罗马字符ax2.set_xticks(np.arange(0, 361, 30),crs=ccrs.PlateCarree())# 经度范围设置，nunpy的作用这不就又来了嘛plt.xticks(fontproperties = 'Times New Roman',size=16) #修改xy刻度字体为新罗马字符plt.yticks(fontproperties = 'Times New Roman',size=16)ax2.set_yticks(np.arange(-90, 90, 15),crs=ccrs.PlateCarree())# 设置yax2.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label = False))#经度0度不加东西ax2.yaxis.set_major_formatter(LatitudeFormatter())# 设置经纬度格式，就是多少度显示那样的，而不是一些数字ax2.set_extent([-178, 178, -70, 70], crs=ccrs.PlateCarree())# 设置空间范围plt.grid(color='k')# 画一个网格吧plt.show()# 显示出图形

那么就运行看看效果吧

如果觉得这个color不喜欢的话，就换一下ncl的来吧，ncl的颜色多而漂亮，喜欢啥就换啥

想要理解这个方法的计算思维，有必要观察原始数据和数据处理之后的样式，理解了数据样式之后可能更有助于我们理解整个程序

import numpy as npimport xarray as xrfrom sklearn.feature_selection import f_regressionimport matplotlib.pyplot as pltimport cartopy.crs as ccrsimport cartopy.feature as cfeaturefrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterimport cmapswith xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\sst.DJF.mean.anom.nc') as f1:      pre = f1['sst_anom'][:-1, :, :]  # 三维数据全取，时间，纬度+经度      lat, lon = f1['lat'], f1['lon']pre2d = np.array(pre).reshape(pre.shape[0], pre.shape[1]*pre.shape[2])#0行代表的个数，1纬度，2经度#pre2d.shape是一个39行，16020列的矩阵，T之后就变为了16020行，39列with xr.open_dataset(r'D:\inuyasha\codeX\codeLEARN\pc.DJF.sst.nc') as f2:      pc = f2['pc'][0, :]#pc是一个39行的数组# # 相关系数pre_cor = np.corrcoef(pre2d.T, pc)[:-1, -1].reshape(len(lat), len(lon))#pre_cor.shape,(16020,)->reshape(89,180)# # 显著性检验# pre_cor_sig = f_regression(np.nan_to_num(pre2d), pc)[1].reshape(len(lat), len(lon))#用0代替NaN# area = np.where(pre_cor_sig < 0.05)nx, ny = np.meshgrid(lon, lat)  # 格网化nx,ny

看看格网化后的经纬度多规范啊。画张图来看看可能也会直观一些。

感谢你能够认真阅读完这篇文章，希望小编分享的“怎么使用Python对NetCDF数据做空间相关分析”这篇文章对大家有帮助，同时也希望大家多多支持编程网，关注编程网Python频道，更多相关知识等着你来学习!