Re: correlating two arrays

From: Erik Noble <enoble_at_nyahnyahspammersnyahnyah>
Date: Fri, 1 Feb 2008 16:37:11 -0500

Hi all.
Thank you for the advice last week about correlating two arrays in NCL.
My question now has to do with using values calculated inside a loop.

If I have to use a loop to calculate a statistic over many timesteps,
how can i "store" each value to plot later? are they automatically
stored in an array?
For instance, my code below uses the escorc command (last function in
code) to calculate correlation coefficients inside of a loop (I have
to use the do loop, sorry). I want to create an array of the r values
so that I can make an XY plot of the r values.
Sincerely,
Erik

; Script for creating correlating V wind from Reanalysis and WRF output
; The purpose of this program is to create files for interpolating
onto a regular grid
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
diri1 ="/Volumes/Data_and_Models/Model-Output/Athena/"
ifile1 ="WRF-SOP3_Athena_3_1_2.nc"
a=addfile(diri1+ifile1,"r")
b = addfile("/Volumes/Data_and_Models/NAMA-SOP3/met_em.d01.2006-09-13_12:00:00.nc","r")
; What times and how many time steps are in the data set?
  times = wrf_user_list_times(a) ; get times in the file
  ntimes = dimsizes(times) ; number of times in the file
 rank = dimsizes(filevardimsizes(a,"XLAT")) ; # of dimensions
; The specific pressure levels that we want the data interpolated to.
  pressure_levels = (/ 850., 700., 500., 300./) ; pressure levels to plot
  nlevels = dimsizes(pressure_levels) ; number of pressure levels
lat=a->XLAT(0,:,:)
lon=a->XLONG(0,:,:)
  do it =348,348 ; TIME LOOP
  ; First get the variables we will need

        p = wrf_user_getvar(a, "pressure",it) ; pressure is our
vertical coordinate
        v = wrf_user_getvar(a, "va",it) ; grid point variable
        u = wrf_user_getvar(a, "ua",it) ; grid point variable

        do level = 0,nlevels-1 ; LOOP OVER LEVELS
                pressure = pressure_levels(level)
                ;u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
                v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False)
                ;uv = table_attach_rows (u_plane, v_plane, 0) ;
Combine both U and V

                v_plane!0 = "lat"
                v_plane!1 = "lon"
               ;U = u_plane(37:55,11:89) ; only if you are specifying
lats and lons
               ;V = v_plane(37:55,11:89)
               ; printVarSummary(V)

                if ( pressure .eq. 700 ) then
                ;WRFv700 = V ; only if you are
specifying lats and lons
                ;FNLv700 = b->VV(0,10,37:55,11:89) ; 700mb =level 10 ;
only if you are specifying lats and lons
                WRFv700 = v_plane
                FNLv700 = b->VV(0,10,1::,:) ; 700mb =level 10

                ;Calculate the unwieghted pattern correlation between
WRF grid and Observation Grid
                r = escorc(ndtooned(WRFv700),ndtooned(FNLv700))
                print(r)
                end if
         end do
end do
end

On Jan 25, 2008 5:46 PM, Dennis Shea <shea_at_ucar.edu> wrote:
>
> Erik Noble wrote:
> > Dear Dennis.
> > Thank you for this.
> > If I am not concerned about weight, (i.e. Small region), then I should
> > proceed with example 5?
> >
> Yes
> > I notice that example 8 does not use the escorc function at all.
> >
> True, but I knew I'd be asked questions on areal weighted
> pattern correlation, so I put it under escorc so I could
> point people to something.
>
> > -Erik
> >
> > On 1/21/08 7:59 PM, "Dennis Shea" <shea_at_ucar.edu> wrote:
> >
> >
> >> Assuming you want a areal weighted pattern correlation, see Example 8
> >> http://www.ncl.ucar.edu/Document/Functions/Built-in/escorc.shtml
> >>
> >>
> >> The "wgt" variable in this example is assumed to be the cos latitude
> >>
> >> lat(lat)
> >> rad = 4.*atan(1.)/180.
> >> wgt = cos(rad*lat)
> >>
> >> Then proceed as in example 8
> >> ===============================
> >> If the lat/lon arrays are 2D, say, XLAT, XLONG
> >>
> >>
> >> wgt = cos(rad*XLAT) ; wgt[*][*]
> >>
> >>
> >> xAvgArea = x*wgt/sum(wgt) ; weighted area average
> >> yAvgArea = y*wgt/sum(wgt) **
> >> <http://www.ncl.ucar.edu/Document/Functions/Built-in/wgt_areaave.shtml>
> >>
> >> xAnom = x - xAvgArea ; anomalies
> >> yAnom = y - yAvgArea
> >>
> >> covxy = *sum*
> >> <http://www.ncl.ucar.edu/Document/Functions/Built-in/sum.shtml>(xAnom*yAnom*wg
> >> t) ; weighted anomaly covariance
> >>
> >> r = covxy/( *sqrt*
> >> <http://www.ncl.ucar.edu/Document/Functions/Built-in/sqrt.shtml>(*sum*
> >> <http://www.ncl.ucar.edu/Document/Functions/Built-in/sum.shtml>(wgt*xAnom^2))
> >> **sqrt* <http://www.ncl.ucar.edu/Document/Functions/Built-in/sqrt.shtml>(*sum*
> >> <http://www.ncl.ucar.edu/Document/Functions/Built-in/sum.shtml>(wgt*yAnom^2))
> >> )
> >>
> >>
> >> Erik Noble wrote:
> >>
> >>> Hi.
> >>>
> >>> I have two 2D arrays:
> >>>
> >>>
> >>>>> whos Model1
> >>>>>
> >>>>>
> >>> Name Size Bytes Class Attributes
> >>>
> >>> Model1 30x79 18960 double
> >>>
> >>>
> >>>
> >>>>> whos Obs
> >>>>>
> >>>>>
> >>> Name Size Bytes Class Attributes
> >>>
> >>> Obs 30x79 18960 double
> >>>
> >>> I want to correlate these two arrays. Which set of NCL functions does this?
> >>>
> >>>
> >>
> >>
> >
> >
>
>
> --
> ======================================================
> Dennis J. Shea tel: 303-497-1361 |
> P.O. Box 3000 fax: 303-497-1333 |
> Climate Analysis Section |
> Climate & Global Dynamics Div. |
> National Center for Atmospheric Research |
> Boulder, CO 80307 |
> USA email: shea 'at' ucar.edu |
> ======================================================
>
>
>
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Feb 01 2008 - 14:37:11 MST

This archive was generated by hypermail 2.2.0 : Mon Feb 04 2008 - 10:52:23 MST