第四章-数据处理

第一章前言及简介第二章NCL变量及基本语法第三章文件输入输出第四章数据处理第五章绘图纲要•代数/逻辑表达式运算符•手动和自动的数据生成•内置函数和程序•自定义函数和程序•计算与元数据•几个常用函数的举例•调用Fortran外部子程序•全域变量•与其他命令结合使用代数运算符代数运算符+加-减＊乘/除^幂#数组乘大于小于用(…)调整计算顺序+是一种特殊的算符－代数算符:5.3+7.9513.25－连接字符串:“pine”+“apple”“pineapple”－代数算符以及连接字符串alpha”+”_+(5.3+7)alpha_12.3”逻辑表达式.le.(小于等于).lt.(小于).ge.(大于等于).gt.(大于).ne.(不等于).eq.(等于).and.(且).xor.(异或).or.(或).not.(非)手动生成数组•用(/…/)–a_integer=(/1,2,3/)–a_float=(/1.0,2.0,3.0/),a_double=(/1,2.0,3.2d/)–a_string=(/abc,”12345,”hello,world/)–a_logical=(/True,False,True/)–a_2darray=(/(/1,2,3/),(/4,5,6/),(/7,8,9/)/)•用new函数[Fortrandimension,allocate;Cmalloc]–x=new(array_size/shape,type,_FillValue)–_FillValue是可选的[若不指定则为默认]–“No_FillValue”指无缺省值–a=new(3,float)–b=new(10,double,1d20)–c=new((/5,6,7/),integer)–d=new(dimsizes(U),string)–e=new(dimsizes(ndtooned(U)),logical)•new和(/…/)可在脚本中任何地方使用–new并不常用自动生成数组•变量之间的赋值y=xy与x有同样的大小，类型和metadata–不需提前分配y数组•赋值（支持的数据格式）u=f-U–数据子集也可:u=f-U(:,3:9:2,:,10:20)•函数返回的数组:不需提前分配–grido=f2fsh(gridi,(/64,128/))–gridi(10,30,73,144)grido(10,30,64,128)grido=f2fsh_Wrap(gridi,(/64,128/));contributed.ncl数组降维•比如T(12,64,128)–Tjan=T(0,:,:)Tjan(64,128)–Tjan自动变为2D:Tjan(64,128)–维数减少–所有的metadata将被复制•可以保留维–Tjan=T(0:0,:,:)Tjan(1,64,128)–TJAN=new((/1,64,128/),typeof(T),T@_FillValue)TJAN(0,:,:)=T(0,:,:)•数组降维是一个“特色”[真的]数组的算符•类似于f90/f95,Matlab,IDL•数组必须同样大小:conform•令A和B为(10,30,64,128)=conform–C=A+B–D=A-B–E=A*B–C,D,E自动生成•令T和P为(10,30,64,128)–theta=T*(1000/P)^0.286theta(10,30,64,128)•令SST为(100,72,144)，且SICE=-1.8–SST=SSTSICE[f90:where(sst.lt.sice)sst=sice]–和的运算(通常)被称为clipping•用built-in函数–令Tbe(30,30,64,128)且P为(30)，则–theta=T*(1000/conform(T,P,1))^0.286•所有数组在运算中遇缺省值即被忽略内置函数和程序;计算变量的zonal和time平均;T(time,lev,lat,lon)T(0,1,2,3);无metadata传输Tzon=dim_avg_n(T,3);Tzon(ntim,klev,nlat)Tavg=dim_avg_n(T,0);Tavg(klev,nlat,mlon)•尽量用•掌握常用的函数–all,any,conform,ind,ind_resolve,dimsizes–fspan,ispan,ndtooned,onedtond,–mask,ismissing,where–system,systemfunc•函数可能需要行列调整–*必须*使用nameddimensions来调整dimsizes(x)•返回数组变量的维数•如果是多维数组，则返回值为1D整型数组fin=addfile(“in.nc”,”r”)t=fin-Tdimt=dimsizes(t)print(dimt)rank=dimsizes(dimt)print(rank=+rank)Variable:dimtType:integerTotalSize:16bytes4valuesNumberofdimensions:1Dimensionsandsizes:(4)(0)12(1)25(2)116(3)100(0)rank=4ispan(start:integer,finish:integer,stride:integer)•返回值则为1D整型数组–以start开始，以finish结束.time=ispan(1990,2001,2)print(time)Variable:timeType:integerNumberofDimensions:1Dimensionsandsizes:(6)(0)1990(1)1992(2)1994(3)1996(4)1998(5)2000fspan(start,finish,npts)b=fspan(-89.125,9.3,100)print(b)Variableb:Type:floatNumberofDimensions:1Dimensionsandsizes:(100)(0)-89.125(1)-88.13081(2)-87.13662(…)….(97)7.311615(98)8.305809(99)9.3•返回值则为1D等间距的浮点型或双精度型数组npts是指包括start和finish所有点的个数ismissingif(any(ismissing(xOrig)))then….else….endifx=(/1,2,-99,4,-99,-99,7/);x@_FillValue=-99xmsg=ismissing(x);print(xmsg)(/False,False,True,False,True,True,False/)•通常和数组函数一起用•if(all(ismissing(x)))then…[else…]endif•nFill=num(ismissing(x))•nVal=num(.not.ismissing(x))maskload$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.nclload$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl”in=addfile(“atmos.nc,r)ts=in-TS(0,:,:)oro=in-ORO(0,:,:);maskocean;[ocean=0,land=1,sea_ice=2]ts=mask(ts,oro,1)•NCL有一个1度的陆地－海洋mask–load$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl”–flagsforocean,land,lake,smallisland,iceshelfwhere;q数组;q0=q=q+256;f90:where(q.lt.0)q=q+256q=where(q.lt.0,q+256,q)•根据条件数组对数组进行赋值•函数where(conditional_expression\•,true_value(s)\•,false_value(s))•类似于f90的“where”•通过数组的条件（True或则False）对各个数进行操作x=where(T.ge.0.and.ismissing(Z),a+25,1.8*b)salinity=where(sst.lt.5.and.ice.gt.icemax\,salinity*0.9,salinity)不可:y=where(y.eq.0,y@_FillValue,1./y)可:y=1.0/where(y.eq.0,y@_FillValue,y)dim_*_n[dim_*]比如:x(time,lat,lon)=x(0,1,2)•函数dim_avg_n(x,n)•xZon=dim_avg_n(x,2)=xZon(ntim,nlat)•xTim=dim_avg_n(x,0)=xTim(nlat,mlon)•函数dim_avg(x)•xZon=dim_avg(x)=xZon(ntim,nlat)•xTim=dim_avg(x(lat|:,lon|:,time|:))=xTim(nlat,mlon)•对数组的行列dimension执行操作•avg,var,sum,pqsort,median,rmsd,…..•dim_*_n函数是新的(v5.1.1以上有)•对指定的行列进行操作•内存使用较少•dim_*函数•仅对最右边的行列操作•可能需要行列调整conform,conform_dimsx(ntim,klev,nlat,mlon),w(nlat);x(0,1,2,3)wx=conform(x,w,2);wx(ntim,klev,nlat,mlon)q=x*wx;q=x*conform(x,w,2)wx=conform_dims((/ntim,klev,nlat,mlon/),w,2)wx=conform_dims(dimsizes(x),w,2)•函数conform(x,r,ndim)•函数conform_dims(dims,r,ndim)•扩展数组(r)使其与(x)数组维数一样(dims)•ndim:标量或数组用以表明x的哪一维或哪些维和r的维数一致•数组r的数值将被‘传播’(复制)至x大小的数组里数组的操作要求数组间conformT(ntim,klev,nlat,mlon),dp(klev);T(0,1,2,3)dpT=conform(T,dp,1);dpT(ntim,klev,nlat,mlon)T_wgtAve=dim_sum_n(T*dpT,1)/dim_sum_n(dpT,1)ind;比如x[*],y[*],z[*][z@_FillValue];生成三个数组仅含‘good’iGood=ind(.not.ismissing(z))xGood=x(iGood)yGood=y(iGood)zGood=z(iGood)•ind仅对1D数组操作–返回值为True的indices–类似于IDL“where”和Matlab“find”[返回indices];比如a[*]ii=ind(a.gt.500)a(ii)=3*a(ii)+2•应当检查返回的indices是否缺省–if(any(ismissing(ii)))then….endifind,ndtooned,onedtond;假如q和x是nD数组q1D=ndtooned(q)x1D=ndtooned(x)ii=ind(q1D.gt.0..and.q1D.lt.5)jj=ind(q1D.gt.25)kk=ind(q1D.lt.-50)x1D(ii)=sqrt(q1D(ii))x1D(jj)=72x1D(kk)=-x1D(kk)*3.14159x=onedtond(x1D,dimsizes(x))•ind仅对1D数组操作–若是nD…则用ndtooned;配合使用