|
+ h: W: o$ W1 T
2021.10.24更新：两个小Tips：1️⃣这个包我现在不用了，毕竟主力平台还是Win，现在推荐使用ProPlot来画图（ProPlot — ProPlot documentation），真的很好用（导师看了都直呼内行，开个玩笑）！2️⃣文中提到的EOF的代码就别看了（属于是早年的腊鸡代码），现在推荐直接使用eofs（Package description — eofs 1.4.0 documentation (ajdawson.github.io)）来完成EOF的计算。
! b! t: n" u1 w+ C# r$ l, b 2021.01.31更新：
: A( W; ~/ t) q2 F m2 r 创建新环境时请将Python版本指定为3.6，否则会冲突导致无法安装。 0 {0 f0 E- s! g
一点感想：PyNgl和PyNio也不再更新了，美帝又去更新自己的Geocat包了，可惜啊Geocat我自己初次尝试，连安装都很大困难，还是学学人家cartopy吧，要不是cartopy画图太丑。
8 x$ }' @! w& @1 E2 K; c 更新： + [5 W. Q5 |) |6 f% M
conda create --name pyncl --channel conda-forge/label/cf201901 pynio pyngl netcdf4 xarray ncl scipy G( K( j" K6 m2 A
此安装代码不建议在使用，速度慢，容易掉线。
4 i# u! k. X' J: v/ B 输入命令：vi ~/.condarc
2 i1 W+ l; Z8 U 讲最后一行的 -default 直接删除(前提是已经增加了清华或中科大源)
$ L: B/ k% s1 j$ R 然后在新环境下使用conda install xxx 依次安装pynio pyngl netcdf4 xarray ncl scipy即可快速安装 & }% B8 m, c7 \) z
如果换源之后还是不行 $ G/ ]$ v. y# t* ]) L! p
数据可视化的软件和库有很多，但是对于大气和海洋数据的可视化，因为涉及到投影、海岸线的绘制，支持的还非常少。 5 D/ H4 Y6 ^! T( S- d7 v: e3 ~; R
主要的有： + t v4 D/ t/ d8 A
①Grads：易用性很差，常常是数据处理计算出来再由Grads读取来绘制。同时语法也比较奇怪。对Win10的支持也不是很好。 & b3 k: a4 R( ^; t9 [, A* i
②NCL：NCL(The NCAR Command Language)是一种专门为大气海洋科学数据处理以及数据可视化设计的高级语言。其实是PyNgl和PyNio的前代产品。 " `* \# _, n! ?$ L6 u6 r: F" s
③MATLAB：很好用但是主攻方向不在地球科学，绘制出来的图片稍微欠缺没敢。(个人感觉) W9 ^ y8 I+ a i5 M' S
④Cartopy：Python下的一个库，由英国气象局开发，Basemap的继承产品(Basemap在2020年不再维护)，主要基于matplotlib实现可视化。目前还不是很完善，容易出现奇奇怪怪的问题。
5 _9 [; q/ W. z; Q ⑤PyNgl：NCL的续作由NCAR负责开发。绘制美观，功能强大，绘图时允许高度自定义。放一张官网扒的图
. D' Y3 M( g9 h; v" A . C" p D4 u9 [. K- v& d
还是挺美观的。 - Y ~- C$ j- y4 Y7 q1 X
为什么选择Python？因为 人生苦短,我用Python！！！ 它语法简单高效，一般来说上手三天即可(但是精通永无止境)，学习曲线很平滑。 - G, @- r8 n3 `9 G0 u# c
既然PyNgl是Python下面的一个库，那么直接用conda安装不就完了吗，理论上来说是的，但是关键在于： - ^! v7 }( r0 W
它不支持Win10！！！
' P9 M) G: p. e# @, q 它不支持Win10！！！
6 f* o5 r0 q; L4 D 它不支持Win10！！！
- Z% n2 f9 R/ b* r ] 这对于不常用Linux的科学工作者或者学生来说是很麻烦的事情，大多数人的主力系统都是Win10，当然macos系统也不需要参考此教程，可以直接利用conda安装。
4 O- G3 [- v0 o+ @8 c0 a 幸运的是微软家有一个很好用的东西：子系统(WSL)，Windows Subsystem for Linux（简称 WSL）是一个为在 Windows 10 上能够原生运行 Linux 二进制可执行文件（ELF 格式）的兼容层。它是由微软与 Canonical 公司合作开发，目标是使纯正的 Ubuntu 14.04 “Trusty Tahr”映像能下载和解压到用户的本地计算机，并且映像内的工具和实用工具能在此子系统上原生运行。 7 W" \/ n1 I2 v5 R
相比于 VMware 等虚拟机，WSL 占用内存和 CPU 资源更少，在 WSL 上运行软件的消耗和直接在 Windows 上差不多。而且，Windows 下可以直接访问 WSL 的环境。 3 V% o0 z3 Y) A% R+ a8 ~; f
安装子系统步骤:
$ E; _6 E/ u8 q/ v$ k. n/ k" S ①开启 WSL 可选特性
+ d2 h% ~4 z9 l8 p- `& q4 g5 _ 在控制面板的“启动或关闭 Windows 功能”中勾选“适用于 Linux 的 Windows 子系统”。 ' |$ k& E# _5 ]9 x( {
; m* u* U& e! [6 Q- @ 或在 PowerShell 中运行下述代码：
8 y9 d! H1 L% f | Enable-WindowsOptionalFeature -Online -FeatureName Microsoft-Windows-Subsystem-Linux
+ L( |6 _, ?7 ?* Q ②下载安装
: b7 v( X, b( {" q; e5 t* r' ]7 G7 _ 打开 Microsoft Store，搜索 Linux，就会显示 Ubuntu、suse 等几个发行版，点击进行安装即可。这里选择了 Ubuntu(第一个)。 - w" N$ ]% y+ U6 @
4 W8 F- [5 L$ T$ |. ^ 下载之后启动菜单里就会出现Ubuntu的图标了。让我们启动它，按照上面的提示等待几分钟，就可以进入初次登陆设置账号的界面。
9 e3 z2 O, e% Z. y- O( l5 X w 图片来自https://wu-kan.cn/_posts/2018-12-14-Windows-Subsystem-for-Linux/设置账号密码就完成了安装 t9 b9 O1 D' l5 T$ a* A2 O
图片来自https://wu-kan.cn/_posts/2018-12-14-Windows-Subsystem-for-Linux/关于Linux的命令可以自行查阅，只需要会cd，ls，pwd等简单命令即可。 1 W# ^# Z0 \6 b+ p. y+ k6 v
需要注意:如无必要不要输入下面命令，因为更新后貌似会有些问题，导致PyNgl无法使用。
4 \7 F7 d, q' Q( B" { sudo apt update -y
! o& A3 V0 U! ~6 M sudo apt upgrade -y # t3 A, G7 e( G2 _
子系统下anaconda的安装：
; w1 j/ Z1 r" b7 r9 L r ①首先从清华镜像Tsinghua Open Source Mirror下载ananconda在linux下的镜像，并放到Win10系统的任意位置(这里假设是D盘) ( ~# z3 ? z w+ A6 d
1 N/ {# c z3 z9 X* B8 _9 S: _ 接着在子系统中输入下面命令 ; P9 y* q/ o) e/ ]+ W4 \
bash /mnt/d/Anaconda3-2020.02-Linux-x86_64.sh#mnt其实是挂载到了本地的硬盘 / w, R; H$ M& J& r
一路回车，注意在正式安装前会弹出安装地址，记下这个位置。然后他会开始安装，等就行。
- @0 a$ G3 Q+ x i' ^ ②环境变量的设置(保证在子系统的任意位置均可使用Python和conda命令)
, N( l) ?- ]+ E( F# L% U7 p$ i7 C) e6 b 1，输入 vi ~/.bashrc (打开了一个编辑器)
7 e# n% y/ k1 e4 _6 K4 ` 2，按一下i键，进入编辑模式
/ {$ |0 m9 W8 S" x9 Q 3，在最后输入export PATH=/home/lishanliao/anaconda3/bin PATH 要换成你自己的目录，如果你是一路回车安装的一般只要把lishanliao换成你自己的用户名即可 . p, X+ c. }9 I
4，按下shift+ i旁边那个冒号键)，输入wq 保存并退出。
6 |, i2 n7 ~# q4 D+ S% Q 5，输入source ~/.bashrc 重新加载环境变量
5 y0 t- m# H; |2 s N 6，如果你害怕不会用这个编辑器，可以直接在命令行中输入下面命令，在执行5即可 7 n0 @# P" h3 I5 l: f1 b
echo export PATH ="/home/liang_yu/anaconda3/bin PATH" >> ~/.bashrc#还是要换成你自己的路径 ③conda下载源更换(加快库的下载速度) & I1 Q' ^/ s" I2 B9 x/ L# G
直接在linux命令行中输入下面命令(复制粘贴即可)
$ w7 R/ |6 j( \3 X% J. g( X- W: T conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/main/
+ d; P5 d# @( Y conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/free/6 B8 z& Q: S" b( F0 r. O- B
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/conda-forge/
( G( t u) A+ o- t/ ~$ V conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/msys2/
- g# N' Y% Q; [1 d2 a5 Z2 C conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/bioconda/3 U8 P' s; ?0 u3 y [
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/menpo/! ~7 A6 d# S8 \4 n
conda config --set show_channel_urls yes
/ S; w1 p w* J- J ④PyNgl的安装 % _* _- o* y1 W$ s# l0 J; K
在linux命令行中输入下面命令，他会自行安装，可能会出现网络错误，要么多试几遍要么挂梯子 4 T9 s8 j6 b, S- G4 ~! D
conda create --name pyncl --channel conda-forge/label/cf201901 pynio pyngl netcdf4 xarray ncl scipy
7 }, x- d: |6 W# M 我并没把库安装在base环境下，而是创建了一个新的环境叫pyncl，因为好像有某些冲突。 * w$ ], A3 P& Z7 E4 V, z* y
安装完成后，输入python，在输入 import Ngl，如果没有报错就是成功了。
: M& W2 f3 e5 w; m8 b+ V: M: @ ⑤接下来就可以画图啦
- R' ?3 U, w* a; Z 现在win10下写好代码，保存成.py文件 # n# f+ Q7 N4 C
首先激活环境，输入下面命令 ; \4 s/ F% Y6 s" r+ y1 _
conda activate pyncl
% X5 v6 ?& e1 g 接着输入,它会自动执行。 " [$ Q6 H# W% s- p/ J0 s
python /mnt/d/file.py 3 m# |, H% Z) I6 e2 ^8 w9 ?5 N
最后展示一下示例图和代码 , F8 S* G& ]" z, K/ M" r- S
- Q8 X0 T C+ w: A) i( ]3 Y 代码示例(前半部分是EOF的计算，可以不看)
, t6 i0 g7 ]+ G z1 y) X$ ^, n """实现了自由调整子图的位置,没有使用panel的面板绘图,使用了新的地图投影"""* R: f4 K2 l$ K- Y
import Ngl
$ `7 ?8 N. h6 `6 b import numpy as np- m+ ]9 ~, i5 z- R8 e4 h/ `& A
import math
5 R) C7 g6 I% l1 P# t# t# Z import netCDF4 as nc
* V3 L- \+ I! ~ Z, m$ Z, ?6 a0 _ def EOF_3d(data_3d,lat,lon,K):
& A3 N8 O. }1 E, p! } """data_3d表示三维数据+ |( r3 @$ |& }' f3 Q* X
例如某一时间段,某一区域范围
S. j: G3 U* c5 B 的SST数据,需先做拉直运算
% a( P* Q8 K9 R K是int表示去前K个模态
5 U+ E5 i& |% c8 T) s8 ` 同时考虑到时空转换K不能大于n! }, u& @1 b+ k4 Z
因为L_Q的形状为[n,n],看表达式就可以看出来"""% G' G$ S& {, Q$ n
lat_n,lon_n,time=data_3d.shape" k5 b; T2 Y. T& O: S+ F l- X5 F; a
m=lat_n*lon_n;n=time#m为个点数,n为观测次数
( G6 p3 }. [3 }% [' c7 W6 E m=m-np.count_nonzero(data_3d[:,:,0]!=data_3d[:,:,0])
+ f" j, x9 c8 A9 O6 e, @: W X=np.empty((m,n))4 t; j- d8 Z/ ?
k=07 v, z; |6 x) b, s+ X
for i in range(lat_n):( I" _9 q1 b. h5 |; C
for j in range(lon_n):2 Q6 U1 W% n0 }: [/ S7 O
if (np.isnan(data_3d[i,j,0])):#考虑陆地
9 g G9 \+ T: x& |' e( U6 t print (land)" @# `. }8 L# d. y( G
else:% Z" I2 Q' a3 e9 d$ w# }8 V3 X
#考虑地球因素乘于cos(纬度)
1 |, _9 e- _8 s4 I9 i- y" B X[k,:]=data_3d[i,j,:]*(math.cos(math.pi*lat[i]/180))**0.5; j8 y) U$ K* v* v/ i1 u* T) Y
k=k+1
8 k7 C; ?$ q# v6 Q! g' X0 l X_mean=X.mean(axis=1).reshape((m,1))+ \" o1 d4 D" n: ^1 M! ?: c
X_d=X-X_mean#利用距平阵计算
6 K y) }* W1 h( E4 n( a if (m<=n):#不使用时空转换, Q& K7 o! w* R, |0 S, h& V* Q
S=np.dot(X_d,X_d.T)/n
6 _7 V* x/ c5 g4 ~% h5 D, o lambda_,L=np.linalg.eig(S)
+ J! |/ T1 b7 r; j; N/ Z$ N b else:4 A+ ^+ _* r& [& w" o5 m
S_Q=np.dot(X_d.T,X_d)/m
! l$ z. ?. g% _$ j' @" j lambda_Q,L_Q=np.linalg.eig(S_Q)#注意L_Q的形状' N/ w& ~9 ?" |# S$ X, C1 ] n
lambda_=lambda_Q*(m/n)( k+ ]& |# T3 Q K4 G
L=np.empty((m,K))
! e1 h# n* T2 Y; d; u+ _' H for i in range(K):% A# i6 o: j( d* r7 j
l_Q=L_Q.real[:,i]
1 \9 B; [ R1 N6 G1 W% @ L[:,i]=np.dot(X_d,l_Q.T)/(m*lambda_Q[i])**0.59 E0 n0 [) {% j( f. U: I: b& s' u
lambda_=lambda_.real[0:K];L=L.real[:,0:K]#避免复数的问题
8 y; a9 H7 }9 Y; F2 E6 T7 M+ K9 V Y=np.dot(L.T,X_d)
6 C. u- j; y3 F& u; M2 [ p=np.empty((K,n))9 R' l" } a8 K7 d, p
e=np.empty((m,K))
* M. q& O: o, V, d! k' d/ t for i in range(K):) ?. S4 f9 \7 B/ a5 N6 E
e[:,i]=lambda_[i]**0.5*L[:,i]" s) e, J" ]9 ?( K% V3 i
p[i,:]=Y[i,:]/lambda_[i]**0.56 Y5 g9 Z( [5 z# E, _6 V% d Y
#将e恢复为空间性+ x1 @- \0 o" K
e_2d=np.empty((lat_n,lon_n,K))! ^$ m9 i& B* X5 ~! e
for m in range(K):
( @3 f6 C. o% a k=0
' o2 v. R% |% ^ for i in range(lat_n):
- L; K: q+ I% P/ Q' T9 r' S for j in range(lon_n):: l* C9 p5 x# Z% R0 N1 N& @
if (np.isnan(data_3d[i,j,0])):3 ?% G0 ^/ C' z& W( z1 K
e_2d[i,j,m]=np.nan! Y. d, K C1 S1 b
else:2 ]8 G' ^/ r: D7 O& D
e_2d[i,j,m]=e[k,m]/(math.cos(math.pi*lat[i]/180))**0.5
% X z/ Y$ I7 g6 H" ~+ ] k=k+1
: w9 j. k) N5 z" }0 q5 T5 Y. _& l3 Q return lambda_,e_2d,p,Y
; h, e% O6 X+ M #------------------
* b& O+ r3 j- P( b' W #时间:2000-2018年每年1月(624开始19个月)
; o8 d8 n5 T3 X" ~: ]8 w3 X5 ^ #区域:0-85,120-2453 l) k! J0 c: V: `; `
#lat[2:37] 35' B5 x: g( V+ j8 o
#lon[48:99] 51
* c5 W! r' v$ ` file="/mnt/c/Users/59799/Desktop/slp.mon.mean.nc"
9 ^9 d4 L& j( P. K2 d! |0 n% i( F #file=rC:\Users\59799\Desktop\slp.mon.mean.nc k8 \4 U b0 m% a9 u
data = nc.Dataset(file)
! j! @% a+ q* t lat = np.array(data.variables[lat])" Q% y3 B' K* I$ @
lon = np.array(data.variables[lon])
! `/ X- a0 A: V' n5 p( m time = np.array(data.variables[time])9 k; d3 n. `5 N, z% N: G1 D Z
slp = np.array(data.variables[slp]); i( B4 L r8 B2 i
data_3d=np.empty((35,51,19))
4 G% [4 L* S: P$ N& a9 @ time=6246 W" E; P: O( x+ E& P. ^
for k in range(19):% I9 R1 D0 U8 t( P! I Y- W) n: d; M
for i in range(35):
' C: m2 E1 p, n" c. D for j in range(51):
6 D6 K" h7 Q5 X! Y! G data_3d[i,j,k]=slp[time,2+i,48+j]% @+ T" _' u: W2 z4 K! Z4 K
time=time+12+ m! n# f7 z+ ?; s
lad,e_2d,p_t,y=EOF_3d(data_3d,lat[2:37],lon[48:99],10)' K. o) Q8 [: p6 w9 o. }* h
6 g9 G- S" Z$ Y #绘图
$ z, n9 T: l4 Q #BlueWhiteOrangeRed
+ O( g7 Q/ |+ X U$ | plot=[]3 a2 i3 v* g! S
wks_type = "png"#子系统没有可视化模块,输出成png图片在打开来看2 j5 h) b, D7 ~, Q1 Z, [
wks = Ngl.open_wks(wks_type,"/mnt/c/Users/59799/Desktop/EOF_first_model")#图片保存到本地的桌面. U$ k1 k$ t9 T" X
cnres = Ngl.Resources(), ?# h! _5 ~5 V4 w* H- r
cnres = Ngl.Resources()
+ h( ?4 C' D! P) j* d# t' y cnres.nglDraw = False#因为要绘制叠加图先阻止画图
9 W i6 }" _& t- N cnres.nglFrame = False#关闭框架与上一条搭配使用
$ [6 U) C; P* w* w1 m1 A& q cmap = Ngl.read_colormap_file("BlueWhiteOrangeRed")#读取colorbar,官网可查啊
! C) R1 B% o" o, A8 g% O cnres.cnFillPalette = cmap[:,:]#不完整使用整个色图,只是用颜色定义值的第五行开始
1 m( [! F& ?+ X0 L% g3 a3 \6 t cnres.cnFillOn = True #等值线填充
3 d# o; {7 | ~1 q+ l cnres.cnLinesOn = False#不画等值线
, @/ y; }3 m& T f cnres.cnLineLabelsOn = False#不画等值线标签(等值线数值)* ^2 s: q2 T( x) {7 E
cnres.lbOrientation = "horizontal"#水平放置colorbar* K& r' @2 c- t. f/ r$ ^8 l
cnres.mpProjection = "Orthographic"# |- f' p5 B. B7 w6 d! o
cnres.sfXArray = lon[48:99]#经度对应为x轴
" L1 v8 T: O& p- ~& T* X7 O2 H7 | cnres.sfYArray = lat[2:37]#纬度对应为y轴
* T' |' B3 W- L cnres.mpCenterLonF = lon[48:99].mean()/ L5 B' O( J; H0 z
cnres.mpCenterLatF = lat[2:37].mean()
7 e, G) R% S) J" Q: r6 D5 x6 l #cnres.mpLimitMode = "LatLon"#规定经纬度范围时需要使用,否则规定无效
& z& k' {0 W+ C/ }) O5 ^) V. l #cnres.mpMinLonF = lon[48:99].min()#不加范围限制则图片为所规定中心纬度经度显示的全球视角(数据只在规定范围有的话会出现无数据的空白区域(海岸线依旧存在))- T" p9 Z) Y, J* M- d, M' V4 c
#cnres.mpMaxLonF = lon[48:99].max()#不加范围限制则图片为所规定中心纬度经度显示的全球视角(数据只在规定范围有的话会出现无数据的空白区域)
" c4 x# X& j9 B9 j% m- W #cnres.mpMinLatF = lat[2:37].min()#不加范围限制则图片为所规定中心纬度经度显示的全球视角(数据只在规定范围有的话会出现无数据的空白区域)
7 R4 _5 Z( r0 _5 w #cnres.mpMaxLatF = lat[2:37].max()#不加范围限制则图片为所规定中心纬度经度显示的全球视角(数据只在规定范围有的话会出现无数据的空白区域)6 X: Q: q+ \6 F4 N M
cnres.nglMaximize = False#必须关闭最大化绘图不然下面的调整plot box的参数是无效的
6 o. T! N9 j. t1 F! n cnres.vpXF = 0.08 # The left side of the box
! C& t8 J ]- l0 h. h; _! f0 _ cnres.vpYF = 0.75 # The top side of the plot box
. H. F% x0 o6 K cnres.vpWidthF = 0.4 # The Width of the plot box+ M' D* r/ @4 E+ F
cnres.vpHeightF = 0.5 # The height of the plot box9 o0 r6 e+ z+ G8 _
cnres.mpGeophysicalLineThicknessF =3.3 #海岸线宽度(粗细程度) 默认的太小了. f" Z/ {& y1 a5 J d
cnres.mpFillOn = False#不开启地图填充,关于陆地海洋内陆水的填充自然不需要
' X% ^5 Z" y! T1 r cnres.tiMainString = "e_first_model"2 o( O/ ?5 w8 @
p1= Ngl.contour_map(wks,e_2d[:,:,0],cnres)
$ B- O# @' ], ^, v
+ Q" k. O: l2 [ xyres = Ngl.Resources()
1 d W' }) z p+ R xyres.nglDraw = False#因为要绘制叠加图先阻止画图
: O2 w2 B/ Z# W xyres.nglFrame = False#关闭框架与上一条搭配使用
. x1 r* k; R; E& Z: a1 ?9 O xyres.nglMaximize = False#必须关闭最大化绘图不然下面的调整plot box的参数是无效的
1 ^" D0 F: k1 @& m H xyres.vpXF = 0.59 # The left side of the box
. s1 X B4 h% R2 K# Z xyres.vpYF = 0.645 # The top side of the plot box
8 \4 q F+ R0 f7 e xyres.vpWidthF = 0.4 # The Width of the plot box O. Z3 Y# E) b( u1 l; F( ~5 S! G" J
xyres.vpHeightF = 0.32 # The height of the plot box
& Q3 j) r; B+ V! Q xyres.tiXAxisString = "Jan of Each Year"0 h4 @( y% Y/ e5 g7 E( L
xyres.tmXBMode = "Explicit" # Define your own tick mark labels.
) M" O# l* R5 u/ w, i7 N, C xyres.tmXBValues = list(range(1,20,1))# Values from 1 to 19." T0 W! \5 o N) s9 L4 F
xyres.tmXBLabels = ["00","01","02","03","04","05","06","07","08","09","10","11","12","13","14","15","16","17","18"]
1 U9 c3 J9 M* e$ p( T! f1 c- j xyres.tiMainString = "p_first_model"
5 E+ B4 x4 T3 b. L+ E x=np.arange(1,21,1)
$ \) ]& r6 p9 C p2 = Ngl.xy(wks,x,p_t[0],xyres)
+ n! c/ T5 N V* H& V' p& B$ d
( R, w! ?6 S+ H8 \/ J, M
" [% r& ^* G& Q+ ^2 f Ngl.draw(p1)
( h- p- V0 \/ t6 c* t# y. l Ngl.draw(p2)
3 G! L2 H/ U( R Ngl.frame(wks)
( f [- @7 m0 j3 `$ I( Q Ngl.end() ! A4 O {$ [ u& J5 P
" ~+ h! M3 ^6 M6 l; r* X( C- O/ v' q, k$ I
" F7 l' y* a9 g! u; @! I
. D7 a+ Y! ]4 L p A5 [
|
