引言:我一直想理解空间相关分析的计算思维,于是今天又拿起python脚本和数据来做练习。首先需要说明的是,这次实验的数据和python脚本均来自于[好久不见]大佬,在跟大佬说明之后,允许我写到公众号来与大家共享,在此对大佬的指点表示感谢,这次实验的脚本可在气象家园或简书app(如果没记错的话)搜索到这次实验的相关内容,也可以微信或者后台发消息给我获取。在此之前我觉得自己还没理解这个方法的计算思维,检验的标准就是我能否迅速运用到其他方面。于是今天又重新回来温习一遍,我把自己的理解与大伙共同交流。
首先,数据的格式是netcdf(.nc)数据,两个数据分别是[哈德来中心海温sst数据,pc数据是对东太平洋ssta做的eof获取]。知道数据信息之后我们就准备开始去运行程序。原始脚本包括了回归分析和相关分析两部分,但是今天我做了空间相关分析这一部分,有兴趣的可以到[好久不见]大佬的气象家园阅读喔!如果还没有安装cartopy包的话请在后台联系我喔
为了方便理解每一步,我选择去jupyter运行,因为可以一段一段程序的运行,这是比较方便的。绘图部分并不是很难,关键还是在于数据预处理部分。
空间相关分析的脚本如下:
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
|
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代替nan
area = 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' ] = 16
ax2 = 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())
# 设置y
ax2.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的颜色多而漂亮,喜欢啥就换啥
想要理解这个方法的计算思维,有必要观察原始数据和数据处理之后的样式,理解了数据样式之后可能更有助于我们理解整个程序
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
|
import numpy as np
import xarray as xr
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
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经度
#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
|
看看格网化后的经纬度多规范啊。画张图来看看可能也会直观一些。
好吧,今天的分享就到这里了,理解了这个计算思维,能更好地迁移运用到其他研究方面,如果还没有安装cartopy包的话请在后台联系我喔,如果需要测试数据和脚本请在后台联系我,当然也可以去[好久不见]大佬的主页。如果觉得这次分享不错的话,还请老铁们点个赞,多多分享,欢迎交流学习,感谢各位!
原始资料:
http://bbs.06climate.com/forum.php?mod=viewthread&tid=92816&highlight=%cf%d4%d6%f8%d0%d4%bc%ec%d1%e9%2b%cf%e0%b9%d8%b7%d6%ce%f6
以上就是如何使用python对netcdf数据做空间相关分析的详细内容,更多关于python对netcdf数据做空间分析的资料请关注服务器之家其它相关文章!
原文链接:https://mp.weixin.qq.com/s/sGLf_3rCFbjtJ2atUkBa2Q