<div dir="ltr"><div><div><div><div>Hi,<br><br></div>I don't really understand your script, but if eof_regres is your regression variable, I would think instead of this:<br><br> file_create->eof_regres = x<br><br></div>you'd want this to write the 1st regression to NetCDF:<br><br> file_create->eof_regres = eof_regress(0,:,:)<br><br></div>HTH...<br></div>Rick<br><br><div><div><br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Sep 19, 2017 at 8:00 AM, Melanie O' hanoly <span dir="ltr"><<a href="mailto:mel.hanoly@gmail.com" target="_blank">mel.hanoly@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hello Ncl Users<div><br><div>I'm trying to extract a variabcle from my routine from ncl to netcdf file.</div><div>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).</div><div>I know that "n" is the number of EOF components, how can I extract only the first?</div><div><br>I'm tryng to use this but the final result is only the first eof Component and not the result of regression..</div><div><br></div><div>";=== output to netcdf file =====</div><div><br></div><div>newfile = "output-EOF.nc"</div><div>system("rm -f "+newfile)</div><div><br></div><div>file_create = addfile(newfile, "c") ;creating file</div><div>file_create->eof_regres = x</div><div><br></div><div>EOF ROUTINE</div><div><pre style="word-wrap:break-word;white-space:pre-wrap">;=============================<wbr>==============================<wbr>===
; eof_4.ncl
;=============================<wbr>==============================<wbr>===
; Match the AAO pattern at:
; <a href="http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/aao/aao.loading.shtml" target="_blank">http://www.cpc.ncep.noaa.gov/<wbr>products/precip/CWlink/daily_<wbr>ao_index/aao/aao.loading.shtml</a>
;
; Use method:
; <a href="http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/aao/aao.shtml" target="_blank">http://www.cpc.ncep.noaa.gov/<wbr>products/precip/CWlink/daily_<wbr>ao_index/aao/aao.shtml</a>
;=============================<wbr>==============================<wbr>===
; Data source was Reanalysis-2 geopotential height.
; <a href="http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.pressure.html" target="_blank">http://www.esrl.noaa.gov/psd/<wbr>data/gridded/data.ncep.<wbr>reanalysis2.pressure.html</a>
;=============================<wbr>==============================<wbr>===
; 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.
;=============================<wbr>==============================<wbr>===
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/<wbr>nclscripts/csm/contributed.<wbr>ncl"
; ==============================<wbr>==============================<wbr>==
; User defined parameters that specify region of globe and
; ==============================<wbr>==============================<wbr>==
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
; ==============================<wbr>==============================<wbr>==
; Open the file: Read only the user specified period and level
; ==============================<wbr>==============================<wbr>==
f = addfile ("./<a href="http://hgt.mon.mean.nc" target="_blank">hgt.mon.mean.nc</a>", "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:<wbr>latN},:) )
printVarSummary(x) ; (time, lat,lon)
; ==============================<wbr>==============================<wbr>==
; compute climatology and Anomalies
; ==============================<wbr>==============================<wbr>==
xClm = clmMonTLL(x) ; (12,lat,lon)
printVarSummary(xClm)
xAnom = calcMonAnomTLL ( x, xClm) ; (time, lat,lon)
printVarSummary(xAnom)
printMinMax(xAnom, True)
; ==============================<wbr>==============================<wbr>=====
; create weights: sqrt(cos(lat)) [or sqrt(gw) ]
; ==============================<wbr>==============================<wbr>=====
rad = 4.*atan(1.)/180.
clat = xAnom&lat
clat = sqrt( cos(rad*clat) ) ; gw for gaussian grid
printVarSummary(clat)
; ==============================<wbr>==============================<wbr>=====
; weight all observations
; ==============================<wbr>==============================<wbr>=====
xw = xAnom*conform(xAnom, clat, 1)
copy_VarMeta(x, xw)
xw@long_name = "Wgt: "+x@long_name
; ==============================<wbr>==============================<wbr>=====
; Reorder (lat,lon,time) the *weighted* input data
; Compute EOFs & Standardize time series
; ==============================<wbr>==============================<wbr>=====
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
; ==============================<wbr>==============================<wbr>=====
; Regress
; ==============================<wbr>==============================<wbr>=====
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
; ==============================<wbr>==============================<wbr>=====
; Extract the YYYYMM from the time coordinate
; associated with eof_ts [same as x&time]
; ==============================<wbr>==============================<wbr>=====
yyyymm = cd_calendar(eof_ts&time,-1)
yrfrac = yyyymm_to_yyyyfrac(yyyymm, 0.0); not used here
;=============================<wbr>==============================<wbr>=
; PLOTS
;=============================<wbr>==============================<wbr>=
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@gsnDraw = False ; don't draw yet
res@gsnFrame = False ; don't advance frame yet
res@gsnPolar = "SH"
res@mpFillOn = False ; turn off map fill
res@mpMaxLatF = latN
res@mpCenterLonF = 180
res@cnFillOn = True ; turn on color fill
res@cnFillPalette = "BlueWhiteOrangeRed"
res@cnLinesOn = False ; True is default
res@cnLineLabelsOn = False ; True is default
res@lbLabelBarOn = False ; turn off individual lb's
; set symmetric plot min/max
symMinMaxPlt(eof_regres, 16, False, res) ; contributed.ncl
res@cnLevelSpacingF = 5.0 ; *special* match CPC
; panel plot only resources
resP = True ; modify the panel plot
resP@gsnMaximize = True ; large format
resP@gsnPanelLabelBar = True ; add common colorbar
resP@gsnPaperOrientation = "portrait" ; force portrait
resP@txString = title+": "+yrStrt+"-"+yrLast
;*****************************<wbr>**************
; first plot
;*****************************<wbr>**************
do n=0,neof-1
res@gsnLeftString = "EOF "+(n+1)
res@gsnRightString = sprintf("%5.1f", eof@pcvar(n)) +"%"
plot(n)=gsn_csm_contour_map_<wbr>polar(wks,eof_regres(n,:,:),<wbr>res)
end do
gsn_panel(wks,plot,(/neof,1/),<wbr>resP) ; now draw as one plot
;*****************************<wbr>**************
; second plot
;*****************************<wbr>**************
; EOF time series [bar form]
rts = True
rts@gsnDraw = False ; don't draw yet
rts@gsnFrame = False ; don't advance frame yet
rts@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@vpHeightF = 0.40 ; Changes the aspect ratio
rts@vpWidthF = 0.85
rts@vpXF = 0.10 ; change start locations
rts@vpYF = 0.75 ; the plot
rts@tiYAxisString = "Standardized" ; y-axis label
rts@gsnYRefLine = 0. ; reference line
rts@gsnXYBarChart = True ; create bar chart
rts@gsnAboveYRefLineColor = "red" ; above ref line fill red
rts@gsnBelowYRefLineColor = "blue" ; below ref line fill blue
; panel plot only resources
rtsP = True ; modify the panel plot
rtsP@gsnMaximize = True ; large format
rtsP@txString = title+": "+yrStrt+"-"+yrLast
; create individual plots
do n=0,neof-1
rts@gsnLeftString = "EOF "+(n+1)
rts@gsnRightString = sprintf("%5.1f", eof_regres@pcvar(n)) +"%"
plot(n) = gsn_csm_xy (wks,yrfrac,eof_ts(n,:),rts)
end do
gsn_panel(wks,plot,(/neof,1/),<wbr>rtsP) ; now draw as one plot
</pre></div></div></div>
<br>______________________________<wbr>_________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@ucar.edu">ncl-talk@ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="http://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">http://mailman.ucar.edu/<wbr>mailman/listinfo/ncl-talk</a><br>
<br></blockquote></div><br></div>