[gradsusr] Model Verification

[Andrew Friedman]

Jeff,
Thanks for posting this really helpful response. I see how I can use this method in my research.
Andrew

On May 1, 2014, at 10:39 AM, Jeff Duda <jeffduda319 at gmail.com> wrote:

> If your forecast and verification data are on the same grid, you can use the maskout and const functions to construct fields of ones and zeros and then use the asum function to add them up to give you a contingency table.  I believe it can be done in a relatively small number of lines.  For example, to compute the correct yes forecasts (i.e., hits), you might use
> 
> *For each threshold
> 'define binforecast = const(const(maskout(rain,rain-threshold),1),0,-u)'
> 'define binobs = const(const(maskout(prec,prec-threshold,1),0,-u)'
> 'define a = const(binforecast+binobs,1)'
> 'define hits = asum(a,[lon/lat range])'
> 
> The first two lines do the same thing, one for the forecast rain, the other for the verification (which I've renamed observed) rain.  Since that first line is a bear, I'll go through it piece by piece:
> ***maskout(rain,rain-threshold)
> The maskout function returns the rain field, except values that fall below a threshold (that you pick) are converted to missing data.  We want missing data because those will be ignored in final calculations.
> ***const(maskout(rain-rain-threshold),1)
> The inner nested const function takes the output from the maskout function and converts all values that are NOT missing to 1.  It does not touch the undefined data.
> ***const(const(maskout(rain,rain-threshold),1),0,-u)
> The outer const function isn't technically necessary, but it completes the construction of a binary field (just ones and zeros) of the precip field where the value is 1 at all grid points where the forecast rain exceeds the threshold and zero otherwise.  This const function takes the output from the previous const function and converts all missing values (the flag -u) to zero.
> 
> We do the same thing for the observed precip field.
> 
> Now for the third line.  I want to create a field that is basically 1 for hits and zero (or undefined) otherwise.
> Since I have two binary fields, I can just add them.  Where there is a hit, I will get 2.  When adding fields, adding a real value to an undefined value at a grid point will turn the whole grid point undefined.  So this field really computes an intersection of areas.  I use the const function to convert the remaining valid values to 1.  Finally, I can compute the total number of hits in the domain by using the asum function.  Make sure you don't use atot() because the latter function incorporates latitude-weighting, which will result in a non-integer sum that isn't representative of your forecast result.
> 
> I wanted to create a true binary field (strictly ones and zeros).  If I changed the valid values to 0.5 instead of to 1 in the first two lines, I wouldn't really need the third line of code.  But this was just to illustrate the process.
> 
> Now use the same general algorithm to code misses, false alarams, and correct "no"s.
> 
> Good luck.
> 
> Jeff Duda
> 
> 
> On Thu, May 1, 2014 at 5:39 AM, Amulya Chevuturi <amulya.chevuturi at gmail.com> wrote:
> Dear Grads Users,
> 
> Recently I asked a question about correction in a script building contingency tables (Asnwered by Dr. Eric Altshuler). I am attaching the script at the end of the mail.
> 
> But this script is for when the model and observation temporal and spatial extent are same. Even the model and observations resolution is same for this script to run. Can anyone tell me how to build 2*2 contingency tables for when the model and observations datasets are at different resolutions and extent. 
> 
> With regards,
> Amulya Chevuturi
> School of Environmental Sciences
> Jawaharlal Nehru University
> New Delhi, India
> 
> 
> ******Script for contingency table*************
> 'reinit'
> 
> 'open model.ctl'
> 'set t 3'
> 'set y 57 208'
> 'set x 72 228'
> 'define a = rain'
> 'close 1'
> 
> 'open observation.ctl'
> 'set t 3'
> 'set y 57 208'
> 'set x 72 228'
> 'define b = prec'
> 
> n1=0
> n2=0
> n3=0
> n4=0
> j=72
> while (j<=228)
> 'set y 'j
> i=57
> while (i<=208)
> 'set x 'i
> 'd a'
> c=subwrd(result,4)
> 'd b'
> d=subwrd(result,4)
> if(c>0&d>0)
> n1=n1+1
> else
> if(c<=0&d<=0)
> n2=n2+1
> else
> if(c<=0&d>0)
> n3=n3+1
> else
> if(c>0&d<=0)
> n4=n4+1
> endif
> endif
> endif
> endif
> i=i+1
> endwhile
> j=j+1
> endwhile
> say n1
> say n2
> say n3
> say n4
> 'close 1'
> 
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
> 
> 
> 
> 
> -- 
> Jeff Duda
> Graduate research assistant
> University of Oklahoma School of Meteorology
> Center for Analysis and Prediction of Storms
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr