[gradsusr] Model Verification

[Jeff Duda]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20140501/8b9136eb/attachment.html