Benutzerspezifische Werkzeuge

##### Sektionen
Sie sind hier: NCL / DKRZ PyNGL example linear regression and movin average
Info
Alle Inhalte des Nutzerportal sind nur auf Englisch verfügbar.

# DKRZ PyNGL example linear regression and movin average

This example shows how to compute the linear regression and moving average with numpy. In addition the x-axis labels are changed to the YYYY-MM-DD format using the advantage of xarray's time handling. The legend is created and placed using the function create_legend.

Example:

`import numpy as npimport xarray as xrfrom   scipy.interpolate import interp1dfrom   scipy import stats,signalimport Ngl, sys,os#-----------------------------------------------------------# Function create_legend()# # create a legend from scratch#-----------------------------------------------------------def create_legend(wks,labels,colors,dx,dy):        nlabs = len(labels)        txres               =  Ngl.Resources()    txres.txFontHeightF =  0.014                #-- default size is HUGE!    txres.txJust        = "CenterLeft"          #-- puts text on top of bars        plres               =  Ngl.Resources()    plres.gsLineThicknessF =  2.0               #-- set line thickness                     x, y = np.array([0.0,0.02]),np.array([0.0,0.0])    dtx, ddy = 0.03, 0.02        for i in range(0,nlabs):        plres.gsLineColor   =  colors[i]        #-- set line color        lg1 = Ngl.polyline_ndc(wks, x+dx, y+dy-(i*ddy), plres)        txres.txFontColor  =  colors[i]         #-- set font color        Ngl.text_ndc(wks,labels[i], x+dx+dtx, y+dy-(i*ddy), txres)#-----------------------------------------------------------#     main program#-----------------------------------------------------------#-- input filediri  = os.path.join(os.environ['HOME'], "NCL/PyNGL/tests/")fname = "rectilinear_grid_2D.nc"                #-- data file name#-- open file and read variablesf = xr.open_dataset(os.path.join(diri,fname))t    = f['tsurf'][:,::-1,:]                     #-- first time step, reverse latitudetime = f['time'][:]                             #-- all timesteps#-- create the x-axis values and labels use xarrays time formatitime  = np.arange(0,len(time))                 #-- x-axis valuesdate   = []                                     #-- empty array to hold the date stringsfor i in range(0,len(time)):    s = str(time[i].values).split('T')    date.append(s[0])                           #-- x-axis labels    print('--> date: ',date)print('--------------------')#-- compute the mean of the 2d field (time,lat,lon)t_mean1 = t.mean()              #-- compute mean over all dimensions; returns 276.90332t_mean2 = np.mean(t)            #-- compute mean over all dimensions; returns 276.90332#-- compute the field mean for each timestep -> time seriet_mean2d = np.array(t.mean(dim=['lat','lon']))print('--> data: ',t_mean2d)print('--------------------------------------')#-- trend computation: create rx and calculate the regression coefficient ryrx  = np.vstack([itime, np.ones(len(itime))]).Tm,t = np.linalg.lstsq(rx, t_mean2d, rcond=None)[0] #-- returns gradient and y-intersectionry  = (m*rx)+tprint('--> trend: ',ry[:,0])print('--------------------------------------')#-- moving average; large window value smooths the linewindow  = 4weights = np.repeat(1.0, window)/windowt_mvavg = np.convolve(t_mean2d, weights, 'valid')print('--> moving avg: ',t_mvavg)print('--------------------------------------')#-------------------------------#--       graphics#-------------------------------lg_labels = ['data',\             'trend - np.linalg',\             'moving average - np.linalg.lstsq (w 4)']  #-- legend labels             lg_colors = ['black','red','blue']              #-- line/legend colors#-- open graphics outputwks = Ngl.open_wks('png','plot_interp_regress_stat')res = Ngl.Resources()                           #-- plot options desiredres.nglDraw           =  Falseres.nglFrame          =  Falseres.nglPointTickmarksOutward = True             #-- point tickmarks outwardres.caXMissingV       =  -999.                  #-- indicate missing valueres.caYMissingV       =  -999.                  #-- indicate missing valueres.trXMinF           =  0.res.trXMaxF           =  itime[-1]res.trYMinF           =  276.4res.trYMaxF           =  277.5res.xyLineThicknessF  =  3                      #-- line thicknessesres.xyLineColor       =  lg_colors[0]           #-- line colorres.tmXBMode          = 'Explicit'              #-- use explicit valuesres.tmXBValues        =  itime[::4]             #-- use the new x-values arrayres.tmXBLabels        =  date[::4]              #-- use the new x-values array as labelsres.tmXBLabelFontHeightF = 0.008res.tmXBLabelAngleF   =  45res.tmXBMinorOn       =  False                  #-- turn off minor tickmark#-- base plot time seriesplot0 = Ngl.xy(wks, itime, t_mean2d, res)       #-- create the plot             #-- plot trendres.xyLineColor       =  lg_colors[1]           #-- line colorres.xyLineThicknessF  =  2                      #-- line thicknessesres.xyDashPattern     =  0                      #-- line dash patternplot1 = Ngl.xy(wks,rx,ry,res)Ngl.overlay(plot0,plot1)#-- plot moving averageres.xyLineColor       =  lg_colors[2]           #-- line colorres.xyLineThicknessF  =  4                      #-- line thicknessesres.xyDashPattern     =  1                      #-- line dash patternplot2 = Ngl.xy(wks,itime,t_mvavg,res)Ngl.overlay(plot0,plot2)#-- draw the plotNgl.draw(plot0)#-- create a legend from scratchcreate_legend(wks,lg_labels,lg_colors,0.58,0.93)#-- advance the frameNgl.frame(wks)Ngl.end()`

Result:

Artikelaktionen