[ncl-talk] Help with netcdf output
Rick Brownrigg
brownrig at ucar.edu
Tue Sep 19 09:22:12 MDT 2017
Hi,
I don't really understand your script, but if eof_regres is your regression
variable, I would think instead of this:
file_create->eof_regres = x
you'd want this to write the 1st regression to NetCDF:
file_create->eof_regres = eof_regress(0,:,:)
HTH...
Rick
On Tue, Sep 19, 2017 at 8:00 AM, Melanie O' hanoly <mel.hanoly at gmail.com>
wrote:
> Hello Ncl Users
>
> I'm trying to extract a variabcle from my routine from ncl to netcdf file.
> The routine calculates the EOF of climate index and then does the
> regression with another variable of my interest and I would like to extract
> to nc this result of regression (The routine is below).
> I know that "n" is the number of EOF components, how can I extract only
> the first?
>
> I'm tryng to use this but the final result is only the first eof Component
> and not the result of regression..
>
> ";=== output to netcdf file =====
>
> newfile = "output-EOF.nc"
> system("rm -f "+newfile)
>
> file_create = addfile(newfile, "c") ;creating file
> file_create->eof_regres = x
>
> EOF ROUTINE
>
> ;==============================================================
> ; eof_4.ncl
> ;==============================================================
> ; Match the AAO pattern at:
> ; http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/aao/aao.loading.shtml
> ;
> ; Use method:
> ; http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/aao/aao.shtml
> ;==============================================================
> ; Data source was Reanalysis-2 geopotential height.
> ; http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.pressure.html
> ;==============================================================
> ; NCL V6.4.0 has new functions eofunc_n_Wrap and eofunc_ts_n_Wrap
> ; that allow you to calculate the EOFs without first having to
> ; first reorder the data. See eof_4_640.ncl.
> ;==============================================================
> ; These files are loaded by default in NCL V6.2.0 and newer
> ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> ; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>
> ; ==============================================================
> ; User defined parameters that specify region of globe and
> ; ==============================================================
> latS = -90.
> latN = -20.
> plev = 700
>
> yrStrt = 1979
> yrLast = 2000
>
> var = "hgt"
> title = str_upper(var)+": "+plev+"hPa "
>
> ymStrt = yrStrt*100 + 1
> ymLast = yrLast*100 + 12
>
> neof = 1 ; Leading EOF only
> optEOF = True
> optETS = False
>
> ; ==============================================================
> ; Open the file: Read only the user specified period and level
> ; ==============================================================
> f = addfile ("./hgt.mon.mean.nc", "r")
>
> TIME = f->time
> YYYY = cd_calendar(TIME,-1)/100 ; entire file
> iYYYY = ind(YYYY.ge.yrStrt .and. YYYY.le.yrLast)
>
> x = short2flt( f->$var$(iYYYY,{plev},{latS:latN},:) )
> printVarSummary(x) ; (time, lat,lon)
>
> ; ==============================================================
> ; compute climatology and Anomalies
> ; ==============================================================
> xClm = clmMonTLL(x) ; (12,lat,lon)
> printVarSummary(xClm)
>
> xAnom = calcMonAnomTLL ( x, xClm) ; (time, lat,lon)
> printVarSummary(xAnom)
> printMinMax(xAnom, True)
>
> ; =================================================================
> ; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
> ; =================================================================
> rad = 4.*atan(1.)/180.
> clat = xAnom&lat
> clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
> printVarSummary(clat)
>
> ; =================================================================
> ; weight all observations
> ; =================================================================
> xw = xAnom*conform(xAnom, clat, 1)
> copy_VarMeta(x, xw)
> xw at long_name = "Wgt: "+x at long_name
>
> ; =================================================================
> ; Reorder (lat,lon,time) the *weighted* input data
> ; Compute EOFs & Standardize time series
> ; =================================================================
> wx = xw(lat|:,lon|:,time|:) ; convenience, cleaner code
> delete(xw)
>
> eof = eofunc_Wrap(wx, neof, optEOF)
> eof = -1*eof ; *special* match sign of CPC
>
> eof_ts = eofunc_ts_Wrap (wx, eof, optETS)
>
> printVarSummary( eof ) ; examine EOF variables
> printVarSummary( eof_ts )
>
> eof_ts = dim_standardize_n( eof_ts, 0, 1) ; normalize
>
> ; =================================================================
> ; Regress
> ; =================================================================
>
> eof_regres = eof ; create an array w meta data
> do ne=0,neof-1
> eof_regres(ne,:,:) = (/ regCoef(eof_ts(ne,:), xAnom(lat|:,lon|:,time|:)) /)
> end do
>
> ; =================================================================
> ; Extract the YYYYMM from the time coordinate
> ; associated with eof_ts [same as x&time]
> ; =================================================================
>
> yyyymm = cd_calendar(eof_ts&time,-1)
> yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
>
> ;============================================================
> ; PLOTS
> ;============================================================
> wks = gsn_open_wks("png","eof") ; send graphics to PNG file
>
> plot = new(neof,graphic) ; create graphic array
> ; only needed if paneling
> ; EOF patterns
>
> res = True
> res at gsnDraw = False ; don't draw yet
> res at gsnFrame = False ; don't advance frame yet
> res at gsnPolar = "SH"
>
> res at mpFillOn = False ; turn off map fill
> res at mpMaxLatF = latN
> res at mpCenterLonF = 180
>
> res at cnFillOn = True ; turn on color fill
> res at cnFillPalette = "BlueWhiteOrangeRed"
> res at cnLinesOn = False ; True is default
> res at cnLineLabelsOn = False ; True is default
> res at lbLabelBarOn = False ; turn off individual lb's
>
> ; set symmetric plot min/max
> symMinMaxPlt(eof_regres, 16, False, res) ; contributed.ncl
> res at cnLevelSpacingF = 5.0 ; *special* match CPC
>
> ; panel plot only resources
>
> resP = True ; modify the panel plot
> resP at gsnMaximize = True ; large format
> resP at gsnPanelLabelBar = True ; add common colorbar
> resP at gsnPaperOrientation = "portrait" ; force portrait
>
> resP at txString = title+": "+yrStrt+"-"+yrLast
>
> ;*******************************************
> ; first plot
> ;*******************************************
> do n=0,neof-1
> res at gsnLeftString = "EOF "+(n+1)
> res at gsnRightString = sprintf("%5.1f", eof at pcvar(n)) +"%"
> plot(n)=gsn_csm_contour_map_polar(wks,eof_regres(n,:,:),res)
> end do
> gsn_panel(wks,plot,(/neof,1/),resP) ; now draw as one plot
>
> ;*******************************************
> ; second plot
> ;*******************************************
> ; EOF time series [bar form]
>
> rts = True
> rts at gsnDraw = False ; don't draw yet
> rts at gsnFrame = False ; don't advance frame yet
> rts at gsnScale = True ; force text scaling
>
> ; these four rtsources allow the user to stretch the plot size, and
> ; decide exactly where on the page to draw it.
>
> rts at vpHeightF = 0.40 ; Changes the aspect ratio
> rts at vpWidthF = 0.85
> rts at vpXF = 0.10 ; change start locations
> rts at vpYF = 0.75 ; the plot
>
>
> rts at tiYAxisString = "Standardized" ; y-axis label
>
> rts at gsnYRefLine = 0. ; reference line
> rts at gsnXYBarChart = True ; create bar chart
> rts at gsnAboveYRefLineColor = "red" ; above ref line fill red
> rts at gsnBelowYRefLineColor = "blue" ; below ref line fill blue
>
> ; panel plot only resources
> rtsP = True ; modify the panel plot
> rtsP at gsnMaximize = True ; large format
> rtsP at txString = title+": "+yrStrt+"-"+yrLast
>
> ; create individual plots
> do n=0,neof-1
> rts at gsnLeftString = "EOF "+(n+1)
> rts at gsnRightString = sprintf("%5.1f", eof_regres at pcvar(n)) +"%"
> plot(n) = gsn_csm_xy (wks,yrfrac,eof_ts(n,:),rts)
> end do
> gsn_panel(wks,plot,(/neof,1/),rtsP) ; now draw as one plot
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170919/bed9e690/attachment.html>
More information about the ncl-talk
mailing list