[ncl-talk] ANN Taylor diagram SUCCESSFUL
Sri Nandini
snandini at marum.de
Wed Mar 15 04:56:11 MDT 2017
Thank you.
After i updated to new version of NCL 6.4.0, my script worked for ANN, DJF, JJA for 2meter temperature and precipitation using taylor stats.ncl and taylor_diagram_cam.ncl.
I have a last questions:
I wish to do "validation or comparison" from my different experiments involving different CESM coupled resolutions (2 and 1degree), and with observations.
Inorder to input more variables, and cases, i use the same script but run this function twice or 3 times depending on the comparision.
stat_taylor = taylor_stats(t2Clm, aClm, clat, 0) ; this gives my cc, ratio and bias for one set of comparisons say model-model
cc(0,:) = stat_taylor(0)
ratio(0,:) = stat_taylor(1)
bias(0,:) = stat_taylor(2)
Then i just repeat for another set of comparison into the loop or outside.
stat_taylor2 = taylor_stats(t3Clm, aClm, clat, 0)
cc(1,:) = stat_taylor2(0)
ratio(1,:) = stat_taylor2(1)
bias(1,:) = stat_taylor2(2)
And pass them into same array for plotting of the cc, ratio and bias.
This works. I just wanted to clarify the method.
Deeply appreciated for all the help :)
On Mar 14, 2017 4:52:30 AM, Dennis Shea wrote:
> Please look at
> > http://www.ncl.ucar.edu/Applications/taylor.shtml
> There have been some changes.
> Read the top and see Examples 7 & 8
>
> ---
> Also, there is a new function called
'taylor_stats' which will be officially added in a future release.
>
> http://www.ncl.ucar.edu/Document/Functions/Contributed/taylor_stats.shtml
> This calculates statistics needed for the Taylor Diagram: pattern_correlation, ratio and bias.
> Note that this is prototyped for two-dimensional arrays.
>
> You will have to explicitly load this:
>
> load "./taylor_stats.ncl"
>
> S\\
>
> On Fri, Mar 10, 2017 at 11:25 AM, Dennis Shea > <> shea at ucar.edu> >> wrote:
> > A Taylor Diagram 'primer' is here:
> >
> > > > http://www-pcmdi.llnl.gov/about/staff/Taylor/CV/Taylor_diagram_primer.pdf
> >
> > I suggest reading it carefully. The summary
> >
> > "In general, the Taylor diagram characterizes the statistical relationship between two fields,
> > a "test" field (often representing a field simulated by a model) and a "reference" field
> > (usually representing “truth”, based on observations). Note that the means of the fields are
> > subtracted out before computing their second - order statistics , so the diagram does not
> > provide information about overall biases, but solely characterizes the centered pattern error.
> >
> > ====
> >
> > ---
> > [0] Every time you use printVarsummary
> >
> > printVarSummary(pcor)
> > printMinMax(pcor,0) ; use this ... are values reasonable
> > print("---")
> > [1] Typically, the ***long-term climatology*** is used NOT yearly data.
> >
> > ; t2(ntim,nlat,mlon), t2(100,360,720)
> > t2Clm = dim_avg_n_Wrap( t2, 0)
> > printVarSummary(t2Clm) ; (lat,lon)
> > printMinMax(t2Clm, 0)
> >
> > [2] Perform the pattern-correlations with the ***long term climatologies.***
> >
> > [3] The ratios are just that... RATIOs: numerator/denominator
> > You did not compute ratios at all.
> > In fact, you created the arrays which are initialized to a deafault _FillValu
> > --
> > ratio = new ((/nCase, nVar/), typeof(p_rat))
> > printVarSummary(ratio);===>ratio[3 3]
> > print (ratio)
> > cc = new ((/nCase, nVar/), typeof(p_cc) )
> > printVarSummary(cc);===>ratio[3 3]
> > print (cc)
> > --
> > But you did not populate the entries with any values. Hence, all _FillValue
> >
> > ==============
> >
> > An excellent. albeit long, NCL tutorial is here. Please loom at this:
> >
> > http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/
> >
> > Good luck
> >
> >
> >
> >
> >
> >
> >
> > On Thu, Mar 9, 2017 at 8:40 AM, Sri Nandini > > <> > snandini at marum.de> > >> > wrote:
> > > Thank you.
> > >
> > > I have now changed the code and modified for annual climatology using month_to_annual ncl function.
> > > I have also inserted 3 model cases of different experiments and only used 2meter temperature as an example to run before inserting other variables.
> > > They are all on same grid and resolution now so interpolation is done on [100 360 720]
> > > I have slight confusion over putting the ratios and pattern correlations into arrays for plotting inside the taylor diagram.
> > > ;==================================================================weighted patterned correlations
> > > pcor = pattern_cor( AIR, t2, clat, 0)
> > > printVarSummary(pcor)
> > > pcor1 = pattern_cor( AIR, t3, clat, 0)
> > > printVarSummary(pcor1)
> > > pcor2 = pattern_cor( AIR, t4, clat, 0)
> > > printVarSummary(pcor2)
> > > print (pcor2)
> > > ;===========================================================ratio of the standardized variances
> > > std1= wgt_arearmse(AIR, t2,clat, 1.0, 0)
> > > printVarSummary(std1)
> > > print (std1)
> > > std2= wgt_arearmse(AIR, t3,clat, 1.0, 0)
> > > printVarSummary(std2)
> > > std3= wgt_arearmse(AIR, t3,clat, 1.0, 0)
> > > printVarSummary(std3)
> > > ;;================================ Put the ratios and pattern correlations into arrays for plotting
> > > nDataSets = 3 ; number of datasets
> > >
> > > p_rat = (/std1,std2,std3/); ===>p_rat[3 100]
> > > printVarSummary(p_rat)
> > > print (p_rat); ===>p_rat[0 247]
> > > p_cc = (/pcor,pcor1,pcor2/); ===>p_cc[3 100]
> > > printVarSummary(p_cc)
> > > print (p_cc); ===>p_cc[0 1];======i have values up to here, after this the missing values occur
> > > ratio = new ((/nCase, nVar/), typeof(p_rat))
> > > printVarSummary(ratio);===>ratio[3 3]
> > > print (ratio)
> > > cc = new ((/nCase, nVar/), typeof(p_cc) )
> > > printVarSummary(cc);===>ratio[3 3]
> > > print (cc)
> > >
> > > All my ratio and cc are with missing values 9.96921e+36 hence apart from the framework plot of taylor diagram the values arent being plotted..i assume i am not passing the arrays well enough or i should average my annual climatology?
> > > Could someone please confirm on this method?
> > >
> > > On Mar 8, 2017 9:45:34 PM, Dennis Shea wrote:
> > > > Overview:
> > > >
> > > > Most commonly, the Taylor diagram is used to compare one-or more data sets (variables) to variables from one or more reference data sets. For example: results from a new model configuration compared to previous model configurations and/or to observational data sets. The title of Karl Taylor's paper is:
> > > >
> > > > Karl E. Taylor > > > >
> > > > Summarizing multiple aspects of model performance in a single diagram> > > >
JGR, vol 106, no. D7, 7183-7192, April 16, 2001> > > > . > > > > The keyword is "multiple". Further, most commonly, *climatologies* of multiple variables from different data sets are compared ... *not* (say) DJF for each year.
> > > >
> > > > --
> > > > NCL's Taylor Diagram was initially developed to compare *multiple* atmospheric model experiments ('cases') with a (say) control run) or observational datasets(say) ECMWF Reanalysis).
> > > >
> > > >
> > > > I made some documentation changes to
> > > > > > > > http://www.ncl.ucar.edu/Applications/taylor.shtml
> > > >
> > > > Specifically, the section:
> > > >
> > > > ======
> > > > The > > > > RATIO> > > > and > > > > CC> > > > arguments have two dimensions.
The left dimension refers to the number of data sets used (eg: model experiments).
The right dimension holds the computedl values. If only one comparison dataset is used, the user
must create a two dimensional array prior to calling the function. EG:
> > > >
CC = conform_dims( (/> > > > 1,1> > > > /), cc)
RATIO = conform_dims( (/> > > > 1,1> > > > /), ratio)
====
> > > > In your case ...
> > > > (a) read 2 data sets with different resolutions and regrid
to same
> > > > t2 = linint2_Wrap(T2w&lon, T2w&lat, T2w, True, wAIR&lon, wAIR&lat, 0)
> > > > wAIR[100 360 720] t2[100 360 720]
> > > > (b) calculate a pattern correlation
> > > > pcor = pattern_cor( wAIR, t2, clat, 0)
> > > > print(pcor) ; 1 dimension of [100](a) read 2 data sets with different resolutions and regrid
to same
> > > > t2 = linint2_Wrap(T2w&lon, T2w&lat, T2w, True, wAIR&lon, wAIR&lat, 0)
> > > > wAIR[100 360 720] t2[100 360 720]
> > > >
> > > > You are comparing the pattern correlations & ratio for each year
> > > > Since there is only one data set with (100) pattern correlations (cc) and ratios (ratio), you must force this into a 2-dimensional array
> > > >
> > > > CC = conform_dims( (/> > > > 1,1> > > > /), cc, 1)
; ===> CC(1,100)
> > > > RATIO = conform_dims( (/> > > > 1,1/)> > > > , ratio, 1> > > > ) > > > > ; ===>RATIO(1,100)
> > > >
> > > > Also, > > > > I think you want an areally weighted root-mean-square-difference
> > > >
> > > > > > > > http://www.ncl.ucar.edu/Document/Functions/Built-in/wgt_arearmse.shtml
> > > >
> > > > =====
> > > >
> > > >
> > > > On Wed, Mar 8, 2017 at 7:49 AM, Sri Nandini > > > > <> > > > snandini at marum.de> > > > >> > > > wrote:
> > > > > Thank you
> > > > > The original data was downloaded from > > > > > https://www.esrl.noaa.gov/psd/data/gridded/data.UDel_AirT_Precip.html> > > > > and i renamed it a year ago. My mistake. It contains 1332monthly timeseries.
> > > > >
> > > > > Both have been used to extract only the needed same time.
> > > > >
> > > > > I have done:
> > > > > (a) read 2 data sets with different resolutions and regrid
to same
> > > > > t2 = linint2_Wrap(T2w&lon, T2w&lat, T2w, True, wAIR&lon, wAIR&lat, 0)
> > > > > wAIR[100 360 720] t2[100 360 720]
> > > > > (b) calculate a pattern correlation
> > > > > pcor = pattern_cor( wAIR, t2, clat, 0)
> > > > > print(pcor) ; 1 dimension of [100]
> > > > >
> > > > > (c) i also calculated the ratio of the standardized variances by:
> > > > > std1= dim_rmsd(wAIR(lat|:,lon|:,time|:), t2(lat|:,lon|:,time|:))
> > > > >
> > > > > printVarSummary(std1) ;[360] x [720]
> > > > > warning:dim_rmsd: 173406 rightmost sections of one or both of the input arrays contained all missing values
> > > > > Error: fatal:Dimension sizes on right hand side of assignment do not match dimension sizes of left hand side
> > > > >
> > > > > The taylor diagram function uses 2d arrays:
> > > > > , RATIO[*][*]:numeric \ ; 2d array: ratios
, CC[*][*]:numeric \ ; 2d array: pattern correlation
> > > > > How can i put the the ratios and pattern correlations into arrays for plotting when i only have 1 dimension for my pattern_corr?
> > > > > Any advice is appreciated.
> > > > >
> > > > >
> > > > >
> > > > > I would like to know if this method of pattern corr and standard deviation normalized is the right one for Taylor diagrams.
> > > > > Deeply appreciated.
> > > > >
> > > > > On Mar 6, 2017 4:56:34 PM, Dennis Shea wrote:
> > > > > > There are many things that just don't add up for me. The code you originally posted 100 time (?months?) steps. The model data shows 1800 months while the NCEP Reanalysis had 748 months. You said you wanted "seasonal (DJF) Taylor diagrams" .... so 1800/12=150 years (150 DJF) and the NCEP has 742/12=62.3 years (62 DJS). I am no sure where you get 100 time steps.
> > > > > > The code you posted indicate
> > > > > > ncl 0> f = addfile("T2M2.nc","r")
> > > > > > ncl 1> T2 = f->TREFHT
> > > > > > ncl 2> printVarSummary(T2)
> > > > > >
> > > > > > Variable: T2
> > > > > > Dimensions and sizes: [time | 1800] x [lat | 96] x [lon | 144] ==> #years= 1800/12=150 years ==> 150 DJF
> > > > > > Coordinates:
> > > > > > time: [ 31..54750]
> > > > > > lat: [ -90.. 90] <================ CAM: CESM ... South-to-North
> > > > > > lon: [ 0..357.5]
> > > > > >
> > > > > >
> > > > > > ncl 5> g = addfile("> > > > > > air.mon.mean.nc> > > > > > ","r")
> > > > > > ncl 6> print(g)
> > > > > > ncl 7> air = short2flt(g->air)
> > > > > > ncl 8> printVarSummary(air)
> > > > > >
> > > > > > Variable: air
> > > > > > Dimensions and sizes: [time | 748] x [level | 17] x [lat | 73] x [lon | 144] ==> #years = 748/12=62.3 ==> 62 DJF
> > > > > > Coordinates:
> > > > > > time: [17067072..17612736]
> > > > > > level: [1000..10]
> > > > > > lat: [90..-90] <============== NCEP: North-to-South
> > > > > > lon: [ 0..357.5]
> > > > > >
> > > > > > ========================
> > > > > > Further, the linint2 documentation states "... > > > > > > nyi> > > > > >
subsection of > > > > > > yi> > > > > > must be monotonically increasing,". The 'yi' are the latitudes. Again, the code you posted did not show that you made the NCEP Reanalysis increasing.
> > > > > >
> > > > > > I have no idea why you you are interpolating to a 5 degree grid. Why not interpolate to the 2.5 grid
> > > > > >
> > > > > >
> > > > > > f = addfile("T2M2.nc","r")
> > > > > > T2 = f->TREFHT
> > > > > > printVarSummary(T2) ; 1800 time steps
> > > > > >
> > > > > > g = addfile("> > > > > > air.mon.mean.nc> > > > > > ","r")
> > > > > > air = short2flt(g->air(:,0,:,:) ; 1000 mb only
> > > > > > printVarSummary(air) ; **look** at 'lat' 90 to -90 North-to_south
> > > > > >
> > > > > > air = air(:,::-1,:) ; make latitudes monotonically increasing (NCL syntax)
> > > > > > printVarSummary(air) ; **look** at 'lat' -90 to 90 South-to-North
> > > > > > ; you have 748 time steps
> > > > > >
> > > > > > ----
> > > > > > Note: These must have the same # of time steps *if* you have a time timension.
> > > > > >
> > > > > > t2 = linint2_Wrap(T2&lon, T2&lat, T2, True, air&lon, air&lat, 0)
> > > > > > printVarSummary(t2) ; 1880 time steps
> > > > > >
> > > > > >
> > > > > > Goos Luck
> > > > > >
> > > > > > On Wed, Mar 1, 2017 at 11:57 PM, Sri Nandini > > > > > > <> > > > > > snandini at marum.de> > > > > > >> > > > > > wrote:
> > > > > > > The files are on the ftp server except > > > > > > > air.mon.mean.nc
> > > > > > >
> > > > > > > I get these errors when uploading:
> > > > > > > local: > > > > > > > air.mon.mean.nc> > > > > > > remote: > > > > > > > air.mon.mean.nc
> > > > > > > 227 Entering Passive Mode (128,117,23,220,192,5).
> > > > > > > 553 Could not create file.
> > > > > > > the file > > > > > > > air.mon.mean.nc> > > > > > > can be found here:
> > > > > > > https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.pressure.html
> > > > > > >
> > > > > > > Thank you.
> > > > > > >
> > > > > > > On Mar 1, 2017 4:10:20 PM, Dennis Shea wrote:
> > > > > > > > re:
> > > > > > > > fo = linint2(airw&lon,airw&lat,airw(time|:,lat|:,lon|:), True, LON, LAT, 0)
> > > > > > > > Which gives an error:
> > > > > > > > warning:linint2: xi, yi, xo, and yo must be monotonically increasing
> > > > > > > > ====
> > > > > > > >
> > > > > > > > > > > > > > > > In your 1st post you had
> > > > > > > >
> > > > > > > > lon = (0., 2.8125, .... , 357,0125) <== ?typo? 357.0125 ??
lat = (-87.8638, ... ,87.8638)
LON = (0., 2.5, ... , 357.5) ; length 144
LAT = (-90.,87.5,...90.) ; length 73
> > > > > > > > If the > > > > > > > > print(airw&lon) and > > > > > > > > print(airw&lat) are increasing then I am not sure what the problem is.
> > > > > > > >
> > > > > > > > > > > > > > > > Really bi-linear interpolation is pretty simple.
> > > > > > > >
> > > > > > > > ===
> > > > > > > > > > > > > > > > Why are you use 'escorc' rather that 'pattern_cor' ?
> > > > > > > >
> > > > > > > > ===
> > > > > > > > > > > > > > > > You can send a *minimal, clean script* and the necessary files to
> > > > > > > >
> > > > > > > > > > > > > > > > ftp > > > > > > > > ftp.cgd.ucar.edu
> > > > > > > > > > > > > > > > anonymous
> > > > > > > > > > > > > > > > your_email
> > > > > > > > > > > > > > > > cd incoming
> > > > > > > > > > > > > > > > put ...
> > > > > > > > > > > > > > > > put ...
> > > > > > > > > > > > > > > > quit
> > > > > > > >
> > > > > > > >
> > > > > > > >
> > > > > > > > On Wed, Mar 1, 2017 at 12:09 AM, Sri Nandini > > > > > > > > <> > > > > > > > snandini at marum.de> > > > > > > > >> > > > > > > > wrote:
> > > > > > > > > Hello
> > > > > > > > >
> > > > > > > > > Thank you.
> > > > > > > > >
> > > > > > > > > Im still at the interpolation stage of the script, what i wish to do is change the model temperature grid (lat,lon) same to my observational temperature.
> > > > > > > > > To do that i followed:
> > > > > > > > > fo = linint2(airw&lon,airw&lat,airw(time|:,lat|:,lon|:), True, LON, LAT, 0)
> > > > > > > > > Which gives an error:
> > > > > > > > > warning:linint2: xi, yi, xo, and yo must be monotonically increasing
> > > > > > > > >
> > > > > > > > > My airw= [100 360 720] ;observed temp
> > > > > > > > > My t2mw= [100 96 144] ' cesm modelled temp
> > > > > > > > >
> > > > > > > > > My fo= [100 360 720] but the error comes when i apply the corr for taylor diag:
> > > > > > > > > error: fatal:escorc: The last dimension of x must be equal to the last dimension of y
> > > > > > > > >
> > > > > > > > > My code for corr is: re=escorc(airw,fo) where both datasets now have same grid.
> > > > > > > > >
> > > > > > > > > Would anyone have a sample code to share if they were successful in running a taylor diagram after interpolating the datasets please? and applied centered corr and rsmd?
> > > > > > > > > > > > > > > > > > Deeply appreciated
> > > > > > > > >
> > > > > > > > > On Feb 28, 2017 7:36:15 PM, Adam Phillips wrote:
> > > > > > > > > > Hi Sri,> > > > > > > > > > The error message is telling you what the issue is:
> > > > > > > > > > > > > > > > > > > > fatal:linint2: The rightmost dimensions of fi must be nyi x nxi, where nyi and nxi are the lengths of yi and xi respectively
> > > > > > > > > > > > > > > > > > > > fi is your input array, and the error message is saying that the rightmost dimensions of fi must be equal to the sizes of the lon and lat arrays that you pass in.
> > > > > > > > > >
> > > > > > > > > > > > > > > > > > > > Thus, if you pass in your t2mw array to linint2 like this:
> > > > > > > > > > > > > > > > > > > > fo = linint2(t2mw&lon,t2mw&lat,t2mw, True, LON, LAT, 0)
> > > > > > > > > > > > > > > > > > > > should work as lon and lat are the 2 rightmost dimensions of t2mw. However, if you pass in your airw array:
> > > > > > > > > > > > > > > > > > > > fo = linint2(airw&lon,airw&lat,airw, True, LON, LAT, 0)
> > > > > > > > > > > > > > > > > > > > You will get the referenced error message as time is the rightmost dimension in airw. The solution is to reorder the dimensions in airw:
> > > > > > > > > > > > > > > > > > > > fo = linint2(airw&lon,airw&lat,airw(time|:,lat|:,lon|:), True, LON, LAT, 0)
> > > > > > > > > >
> > > > > > > > > > > > > > > > > > > > Hope that helps. If you have any further questions please respond to the ncl-talk email list and not to me personally.
> > > > > > > > > > > > > > > > > > > > Adam
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > > On Tue, Feb 28, 2017 at 8:43 AM, Sri Nandini > > > > > > > > > > <> > > > > > > > > > snandini at marum.de> > > > > > > > > > >> > > > > > > > > > wrote:
> > > > > > > > > > > Thank you
> > > > > > > > > > >
> > > > > > > > > > > So i have cleaned out the commented statements, and the errors i am getting are these:
> > > > > > > > > > >
> > > > > > > > > > > fatal:linint2: The rightmost dimensions of fi must be nyi x nxi, where nyi and nxi are the lengths of yi and xi respectively
> > > > > > > > > > >
> > > > > > > > > > > printVarSummary(t2mw)
> > > > > > > > > > > Dimensions and sizes: [time | 100] x [lat | 192] x [lon | 288]
> > > > > > > > > > > printVarSummary(airw)
> > > > > > > > > > > Dimensions and sizes: [lat | 62] x [lon | 162] x [time | 100]
> > > > > > > > > > > It stops at the interpolation stage before even coming to the correlation section.
> > > > > > > > > > > Basically i am unsure of my interpolation method between observations and model output data.
> > > > > > > > > > >
> > > > > > > > > > > On Feb 24, 2017 4:28:11 PM, Dennis Shea wrote:
> > > > > > > > > > > > I started to look at your code. However, it has lots of commented statements and, to me, it is not clear what is happening. Really, you should send only clean scripts. We like to help but we don't have the time to > > > > > > > > > > > > decipher codes.
> > > > > > > > > > > >
> > > > > > > > > > > > The commented
> > > > > > > > > > > >
> > > > > > > > > > > > ;cor1 = dim_avg_n_Wrap(pattern_cor( t2mw, airw, clat, 0), 0)
> > > > > > > > > > > >
> > > > > > > > > > > > contains
> > > > > > > > > > > > > > > > > > > > > > > > http://www.ncl.ucar.edu/Document/Functions/Contributed/pattern_cor.shtml
> > > > > > > > > > > >
> > > > > > > > > > > > What is wrong with the following:
> > > > > > > > > > > >
> > > > > > > > > > > > printVarSummary(t2mw)
> > > > > > > > > > > > printVarSummary(airw)
> > > > > > > > > > > > pcor = pattern_cor( t2mw, airw, clat, 0)
> > > > > > > > > > > > print(pcor)
> > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > > > > > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > > > On Thu, Feb 23, 2017 at 9:03 AM, Sri Nandini > > > > > > > > > > > > <> > > > > > > > > > > > snandini at marum.de> > > > > > > > > > > > >> > > > > > > > > > > > wrote:
> > > > > > > > > > > > > Dear NCL community,
> > > > > > > > > > > > > Greetings!
> > > > > > > > > > > > >
> > > > > > > > > > > > > I am trying to plot seasonal (DJF) Taylor diagrams and have errors in interpolation my datasets on same grid.
> > > > > > > > > > > > > (a) read 2 data sets with different resolutions and regrid
(modelled and obs temperature for a test)
> > > > > > > > > > > > > (b) calculate a pattern correlation
> > > > > > > > > > > > >
> > > > > > > > > > > > > My script is attached below ::
> > > > > > > > > > > > >
> > > > > > > > > > > > > ;=============================> > > > > > > > > > > > > =====================================
> > > > > > > > > > > > > ;Taylor diagram calculations
> > > > > > > > > > > > > ;================================ interpolation onto common grid (of observational data)
> > > > > > > > > > > > > ;ncl pattern_cor between different size arrays
> > > > > > > > > > > > > t2mw = f->t2mw(0,0,:,{0:360}); remove cyclic point
> > > > > > > > > > > > > lon = f->lon({0:360})
> > > > > > > > > > > > > ;************************************************
> > > > > > > > > > > > > ; interpolate to new grid
> > > > > > > > > > > > > ;***********************************************
> > > > > > > > > > > > > newlat = fspan(-60.,60,24)
> > > > > > > > > > > > > newlon = fspan(0.,355.,72)
> > > > > > > > > > > > >
> > > > > > > > > > > > > newt2mw = linint2_Wrap(lon,t2mw&lat,t2mw,True,newlon,newlat,0)
> > > > > > > > > > > > > newt2mw!0 ="lat"
> > > > > > > > > > > > > newt2mw!1 = "lon"
> > > > > > > > > > > > > newt2mw&lat = newlat
> > > > > > > > > > > > > newt2mw&lon = newlon
> > > > > > > > > > > > > ;=============================> > > > > > > > > > > > > =========================================centered Pattern correlation (coslat weighting has been done previously above)
> > > > > > > > > > > > > re=escorc(airw,newt2mw)
> > > > > > > > > > > > > ;cor1 = dim_avg_n_Wrap(pattern_cor( t2mw, airw, clat, 0), 0);rc = pattern_cor(x, y,gw, 0) ; gaussian weighting, centered
> > > > > > > > > > > > > mmd= (/cor1/)
> > > > > > > > > > > > > printVarSummary(mmd)
> > > > > > > > > > > > > ;================================Standard deviation ;================================
> > > > > > > > > > > > >
> > > > > > > > > > > > > ;pre0_Std = dim_avg_n_Wrap( dim_stddev_n_Wrap( t2mw, (/1,2/)), 0)
> > > > > > > > > > > > > ;std1 = dim_rmsd_Wrap(airw,t2mw, 0);computes rootmean square difference
> > > > > > > > > > > > > ;std1 = dim_rmsd_n(t2mw, airw, 0);
> > > > > > > > > > > > > std1 = dim_rmsd( t2mw(lat|:,lon|:,time|:), airw(lat|:,lon|:,time|:) ) ; ==> rmsdTime(nlat,nlon)
> > > > > > > > > > > > >
> > > > > > > > > > > > > ;rmsdTime = dim_rmsd_n( x, y, 0 ) ; ==> no reordering needed
> > > > > > > > > > > > > ================================================================
> > > > > > > > > > > > > I had a look at other interpolation functions for correlations between different grids such a s:
> > > > > > > > > > > > >
Assume > > > > > > > > > > > > > fi> > > > > > > > > > > > > is a 4D array dimensioned > > > > > > > > > > > > > ntim> > > > > > > > > > > > > x
> > > > > > > > > > > > > nlvl> > > > > > > > > > > > > x > > > > > > > > > > > > > nlat> > > > > > > > > > > > > x > > > > > > > > > > > > > mlon> > > > > > > > > > > > > (> > > > > > > > > > > > > ntim> > > > > > > > > > > > > =50,
> > > > > > > > > > > > > nlvl> > > > > > > > > > > > > =30, > > > > > > > > > > > > > nlat> > > > > > > > > > > > > =64, > > > > > > > > > > > > > mlon> > > > > > > > > > > > > =128), and that
the
rightmost dimension is to be treated as cyclic (the user should
not
add a cyclic point for the rightmost dimension).> > > > > > > > > > > > >
> > > > > > > > > > > > >
All times and levels will be interpolated and returned in a new
array
> > > > > > > > > > > > > fo> > > > > > > > > > > > > dimensioned > > > > > > > > > > > > > ntim> > > > > > > > > > > > > x > > > > > > > > > > > > > nlvl> > > > > > > > > > > > > x > > > > > > > > > > > > > 73> > > > > > > > > > > > >
x
> > > > > > > > > > > > > 144> > > > > > > > > > > > > :
> > > > > > > > > > > > >
> > > > > > > > > > > > > lon = (0., 2.8125, .... , 357,0125)
lat = (-87.8638, ... ,87.8638)
LON = (0., 2.5, ... , 357.5) ; length 144
LAT = (-90.,87.5,...90.) ; length 73
fo = > > > > > > > > > > > > > linint2_Wrap> > > > > > > > > > > > > (lon,lat,fi, True, LON,LAT, 0)
> > > > > > > > > > > > > Error:
> > > > > > > > > > > > >
> > > > > > > > > > > > > Deeply appreciated
> > > > > > > > > > > > >
> > > > > > > > > > > > > _______________________________________________
> > > > > > > > > > > > >
ncl-talk mailing list
> > > > > > > > > > > > > ncl-talk at ucar.edu
> > > > > > > > > > > > >
List instructions, subscriber options, unsubscribe:
> > > > > > > > > > > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > > > > > > > > > > > >
> > > > > > > > > > > >
> > > > > > > > > > >
> > > > > > > > > > > _______________________________________________
> > > > > > > > > > >
ncl-talk mailing list
> > > > > > > > > > > ncl-talk at ucar.edu
> > > > > > > > > > >
List instructions, subscriber options, unsubscribe:
> > > > > > > > > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > > --
> > > > > > > > > > Adam Phillips
> > > > > > > > > > > > > > > > > > > > Associate Scientist, > > > > > > > > > > Climate and Global Dynamics Laboratory, NCAR
> > > > > > > > > >
> > > > > > > > > > > > > > > > > > > > www.cgd.ucar.edu/staff/asphilli/> > > > > > > > > > > > > > > > > > > > 303-497-1726> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > > > > > > > > > > > >
> > > > > > > > > >
> > > > > > > > > > > > > > > > > >
> > > > > > > > > _______________________________________________
> > > > > > > > >
ncl-talk mailing list
> > > > > > > > > ncl-talk at ucar.edu
> > > > > > > > >
List instructions, subscriber options, unsubscribe:
> > > > > > > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > > > > > > > >
> > > > > > > >
> > > > > > > > > > > > > >
> > > > > > > _______________________________________________
> > > > > > >
ncl-talk mailing list
> > > > > > > ncl-talk at ucar.edu
> > > > > > >
List instructions, subscriber options, unsubscribe:
> > > > > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > > > > > > > > >
> > > > > _______________________________________________
> > > > >
ncl-talk mailing list
> > > > > ncl-talk at ucar.edu
> > > > >
List instructions, subscriber options, unsubscribe:
> > > > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > > > > p_rat = (/std1,std2,std3/); ===>p_rat[3 100]9.96921e+36
> > > > > (0,0) 9.96921e+36
> > > > > (0,1) 9.96921e+36
> > > > > (0,2) 9.96921e+36
> > > > > (1,0) 9.96921e+36
> > > > > (1,1) 9.96921e+36
> > > > > (1,2) 9.96921e+36
> > > > > (2,0) 9.96921e+36
> > > > > (2,1) 9.96921e+36
> > > > > (2,2) 9.96921e+36
> > > > >
> > > >
> > > >
> > > >
> > > > > > > >
> > >
> > > _______________________________________________
> > >
ncl-talk mailing list
> > > ncl-talk at ucar.edu
> > >
List instructions, subscriber options, unsubscribe:
> > > http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > >
> >
> >
> >
>
undef("taylor_stats")
function taylor_stats(t[*][*]:numeric, r[*][*]:numeric, w:numeric, opt[1]:integer)
;
; Calculate statistics which are used by Taylor diagram
;
; Input:
; t - test array (lat,lon)
; r - reference array (eg: model control run; ECMWF reanalysis)
; same size, shape and ordering as 't'
; w - weights: latitudinal weights: gaussian, cos(lat), area(nlat)
; + w[1] - no weighting used
; + w[*] - rectilinear grid: w(lat)
; + w[*][*] - curvilinear grid; w(lat,lon); match 't' and 'r'
; opt - option to select what will be returned
; 0: return (/ pattern_correlation, ratio, bias /)
; 1: return (/ pattern_correlation, ratio, bias, tmean, rmean, tvar, rvar, rmse/)
;-------------------------------------------------------------------------------------------
; where:
; {x/r}mean - area weighted means for the 'test' and 'reference' arrays
; {x/r}var - area weighted variances for the 'test' and 'reference' arrays
; rmse - area weighted root mean square of grid-point differences
;
local dimt, dimr, dimw, rankr, rankt, rankw, tmean, rmean, tdiff, rdiff, rmse \
, rmsdiff, tw, wsum, tvar, rvar, ratio, bias, stats_taylor
begin
dimt = dimsizes(t)
dimr = dimsizes(r)
dimw = dimsizes(w)
rankt = dimsizes(dimt)
rankr = dimsizes(dimr)
rankw = dimsizes(dimw)
; error checking
if (rankr.lt.2 .or. rankt.lt.2) then
print("taylor_ratio: rank must be > 2: rank(r)="+rankr+" rank(t)="+rankt)
exit
end if
if (rankr.ne.rankt) then
print("taylor_ratio: rank(x), rank(r) do not match: rank(r)="+rankr+" rank(t)="+rankt)
exit
end if
if (any(dimt.ne.dimr)) then
print("taylor_ratio: dimension sizes must match")
print("---")
print("taylor_ratio: dimt="+dimt)
print("---")
print("taylor_ratio: dimr="+dimr)
exit
end if
if (rankw.gt.2) then
print("taylor_ratio: rank(w) must be <=2 : rank(w)="+rankw)
exit
end if
; centered pattern correlation: pc
; centered areal weighted means: tmean, rmean
; difference: (t-tmean), (t-rmean)
; squared difference: (t-tmean)^2, (r-rmean)^2
; weighted squared difference: w*(t-tmean)^2, w*(r-rmean)^2
; sum of weights ==> wsum
; weighted mean variance = SUM[w*(t-tmean)^2)]/wsum
; SUM[w*(r-rmean)^2)]/wsum
; All above are relative to a centered mean
; The following is a statistic that measure weighted *grid point differences*
; weighted grid-point root-mean-square-error: w*(t(j,i)-r(j,i))^2/sumw
pc = pattern_cor(t, r, w, 0) ; centered pattern correlation: opt_pc=0
if (rankw.eq.1) then ; RECTILINEAR GRID: w[1] or w[*]
tmean = wgt_areaave(t, w, 1.0, 0) ; area weighted means
rmean = wgt_areaave(r, w, 1.0, 0)
if (isscalar(w)) then
tw = conform(t, w,-1) ; w[1] = => tw[*][*]
else
tw = conform(t, w, 0) ; w[*] = => tw[*][*]
end if
wsum = sum(tw) ; sum of weights
tdiff = t-tmean ; diff from 'test' weighted centered mean
rdiff = r-rmean ; 'reference'
tvar = sum(tw*tdiff^2)/wsum ; mean area wgted difference: SUM[wgt*tdiff^2]/SUM[wgt]
rvar = sum(tw*rdiff^2)/wsum
rmse = wgt_arearmse(t,r, w, 1.0, 0) ; area weighted individual grid point differences
else ; CURVILINEAR GRID: w[*][*]
tmean = wgt_areaave2(t, w, 1) ; area weighted means
rmean = wgt_areaave2(r, w, 1)
wsum = sum(w)
tdiff = t-tmean ; difference from central weighted mean
rdiff = r-rmean ; " " " " "
tvar = sum(w*(tdiff^2))/wsum ; mean area wgted difference: SUM[wgt*tdiff^2]/SUM[wgt]
rvar = sum(w*(rdiff^2))/wsum ;
rmse = wgt_arearmse2(t,r, w, 0) ; area weighted individual grid-point differences
end if
bias = tmean-rmean ; test - reference
if (rmean.ne.0) then
bias = (bias/rmean)*100 ; bias [%]
else
bias = totype(-999, typeof(bias))
bias at _FillValue = bias
end if
bias at long_name = "bias: [(tmean-rmean)/rmean]*100"
bias at units = "%"
if (rvar.ne.0) then
ratio = sqrt(tvar/rvar)
else
ratio = totype(-999, typeof(ratio))
ratio at _FillValue = ratio
end if
ratio at long_name = "RATIO: Taylor Diagram"
if (opt.eq.0) then
stats_taylor = (/pc, ratio, bias/)
stats_taylor at long_name = "0-pattern_cor; 1-ratio; 2-bias (%)"
else
stats_taylor = (/pc, ratio, bias, tmean, rmean, tvar, rvar, rmse/)
stats_taylor at long_name = "0-pattern_cor; 1-ratio; 2-bias (%); " +\
"3-tmean; 4-rmean; 5-tvar; 6-rvar; 7-rmse"
end if
return(stats_taylor)
end
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20170315/5922f8e2/attachment.html
More information about the ncl-talk
mailing list