[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

20060123: bugs in gdprof and gdpoint



David,

It is true that the field of interpolated squared values is not the same
as the square of the grid values interpolated to a point. Eg, the
operation is not commutative (linear interpolation of a quadratic field).

Looking at GDPROF,  it obtains the vector component grids, and interpolates
the components independently (using bilinear interpolation) for the GVECT 
field, 
and these components are output in the text listing.

Your assertion is that the GFUNC scalar computation of direction and/or 
magnitude
is also correct in that the program composes the scalar field first, eg
for mag(wnd), the scalar magnitude of the wind field is calculated
and the resultant field is interpolated. To be the same as the constitution
of the GVECT components, each vector would have to be computed and interpolated 
first.
This would be more like mag(V) where V = squo(2,vadd(wnd,sm5v(wnd))) where the
first order derivative of the components is now implicit.

So, at grid points, the composed grid from GVECT components does equal the
interpolated GFUNC value, while in between grid points, they do not.
The GR_INTP function only implements bilinear interpolation at present. A higher
order scheme using cubic splines would be possible, but wouldn't solve the 
problem
at higher orders.

One way to work around this is to use gddiag to create a new grid of 
interpolated
components to the desired grid spacing. The current distribution allows GDDIAG 
to 
create the new gridfile, and specify the grids to be created as vector 
components 
using the GRDTYP parameter. That will allow your GFUNC interpolation to match 
the
components of the GVECT interpolation at the new points. You can then compose 
the
vectors from the interpolated points rather than interpolating the
scalar derived field as is done in GFUNC for GDPROF.

Steve Chiswell
Unidata User Support


>From: David Ovens <address@hidden>
>Organization: UCAR/Unidata
>Keywords: 200601231935.k0NJZV7s017961

>Hello,
>
>I have been asked by another one of our group to follow up on this bug
>report.  Have you been able to verify these errors?  Have these been
>passed on to be worked on by anyone?  Is there any timeline
>established for when they might be fixed? 
>
>I realize that my original emails may not have stated the errors very
>well.  In summary, wind direction calculations (and simple math
>functions involving 2 wind components) at interpolated grid
>points in gdprof and gdpoint are WRONG.  The reason is that wind
>directions (and math involving 2 wind components) are calculated at
>gridpoints first and are then interpolated without regard to vector
>componets (see original email below for examples of the erroneous
>output).  The correct way to do the calculation is to interpolate 
>components and then do the all math on the interpolated values.
>Though one could argue against this approach for something like 
>  GFUNC = MUL(urel,urel)
>
>you really cannot take the other approach for
>  GFUNC = dirr(wnd)
>  GFUNC = dirr(vecr(urel,vrel))
>
>I have only made a quick, single point check in gdthgt and have not
>tried gdtser.  The bug may exist in those as well.
>
>Thanks for any info,
>
>David
>
>
>On Thu, Nov 24, 2005 at 11:56:12PM -0800, ovens wrote:
>> Hello,
>> 
>> Turns out the math functions are occurring on the data BEFORE it is
>> interpolated.  Hence, at point (x,y) = (8.458,92.63) even though 
>>     UREL = -1.116
>> GFUNC  = MUL(urel,urel) gives    
>>     MUL(urel,urel) = 3.763 
>> instead of 1.25 because the surrounding 4 gridpoints:
>>      (x,y)                   UREL                  VREL
>>  -------------    ------------------       --------------
>>  (8,93) (9,93)     1.0634    -2.4366       0.379    1.129
>>  (8,92) (9,92)    -1.8116    -2.4366       1.879   -0.371
>> 
>> are each being squared and THEN interpolated.  At an exact grid point,
>> you get the correct/expected answer always.
>> 
>> Though I understand what is going on in the programs, I am still
>> confused about the vector output.  For the same point (8.458,92.63),
>> the  
>>     UREL = UWND = -1.116
>>     VREL = VWND =  0.769
>> and the program says that 
>>     GFUNC = dirn(wnd) =  159.70622
>> instead of 124.6 degrees (which you get from urel and vrel by using 
>> dirn = 270 - atan2(urel,vrel)*RTD).
>> 
>> What do you think the wind direction should be?
>> 
>> 
>> David
>> 
>> 
>> 
>> On Thu, Nov 24, 2005 at 03:01:07AM -0800, ovens wrote:
>> > Hello,
>> > 
>> > I have discovered the following bugs in GEMPAK 5.8.3a (Unidata binary
>> > version and my compiled and patched versions).
>> > 
>> > In gdprof and gdpoint the use of MUL(urel,urel) and MAG(WND) and
>> > DIRN(WND) are giving totally incorrect answers (for ECMWF, GFS, and
>> > ETA 104 grids at least).  When I run MUL(urel,const) I get the
>> > correct number.  This looks like some kind of memory error when both
>> > arguments are the same item.
>> > 
>> > Here are the commands and results from gdprof, gdpoint results follow:
>> > ====================  gdprof results ======================
>> > GEMPAK-GDPROF>li
>> >  GPOINT   = @8.458;92.63
>> >  GDATTIM  = f0
>> >  GVCORD   = hagl
>> >  GVECT    = WND
>> >  GDFILE   = /home/disk/data/gempak/model/ecmwf/2005112300_ecmwf.gem
>> >  LINE     = 3
>> >  MARKER   = 0
>> >  BORDER   = 1
>> >  PTYPE    = lin
>> >  SCALE    = 0
>> >  XAXIS    = 0
>> >  YAXIS    = 0/10
>> >  WIND     = BM1
>> >  REFVEC   =  
>> >  WINPOS   = 1
>> >  FILTER   = YES
>> >  TITLE    = 1
>> >  PANEL    = 0
>> >  CLEAR    = YES
>> >  TEXT     = 1/21//hw
>> >  DEVICE   = XW
>> >  OUTPUT   = T
>> >  THTALN   = 0
>> >  THTELN   = 0
>> >  MIXRLN   = 0
>> >  GFUNC    = urel
>> > r
>> > 
>> > 
>> >  Grid function:             urel
>> >  X,Y grid point:              8.46   92.63
>> >  Lat/lon:                    47.44 -122.31
>> >  Number of levels:            1
>> >  Data range:                 -1.12   -1.12
>> >  Scaling factor:             10** 0
>> >        051123/0000V000  UREL   @8.458;92.63
>> >        LATITUDE, LONGITUDE = (  47.44, -122.31 )
>> >                   HAGL                  UREL
>> >                   10.00               -1.116
>> >                   HAGL        U         UREL   V         VREL
>> >                   10.00               -1.116            0.769
>> > 
>> >  GEMPAK-GDPROF>gfunc = mul(urel,urel)
>> >  GEMPAK-GDPROF>r
>> >  GDPROF PARAMETERS
>> >  Grid function:             mul(urel,urel)
>> >  Grid vector function:      WND
>> >  X,Y grid point:              8.46   92.63
>> >  Lat/lon:                    47.44 -122.31
>> >  Number of levels:            1
>> >  Data range:                  3.76    3.76    ***** WRONG!!!!
>> >  Scaling factor:             10** 0
>> >  Wind display:            BM1                                             
>                                                                              
>    
>> >        051123/0000V000  MULRELUREL   @8.458;92.63
>> >        LATITUDE, LONGITUDE = (  47.44, -122.31 )
>> >                   HAGL            MULRELUREL
>> >                   10.00                3.763    ***** WRONG!!!!
>> >                   HAGL        U         UREL   V         VREL
>> >                   10.00               -1.116            0.769
>> > 
>> >  GEMPAK-GDPROF>gfunc = mul(urel,-1.116)
>> >  GEMPAK-GDPROF>r
>> >  Grid function:             mul(urel,-1.116)
>> >  Data range:                  1.25    1.25
>> >        051123/0000V000  MULREL-1.1   @8.458;92.63
>> >        LATITUDE, LONGITUDE = (  47.44, -122.31 )
>> >                   HAGL            MULREL-1.1
>> >                   10.00                1.246
>> >                   HAGL        U         UREL   V         VREL
>> >                   10.00               -1.116            0.769
>> > 
>> > 
>> > ====================  gdpoint results ======================
>> >  GEMPAK-GDPOINT>li
>> >  GDATTIM  = f0
>> >  GDFILE   = /home/disk/data/gempak/model/ecmwf/2005112300_ecmwf.gem
>> >  GLEVEL   = 10
>> >  GPOINT   = @8.458;92.63
>> >  GVCORD   = hagl
>> >  GFUNC    = mul(urel,-1.116)
>> >  SCALE    = 0
>> >  GEMPAK-GDPOINT>gfun = urel
>> >  GEMPAK-GDPOINT>r
>> >  
>> > GDFILE: /home/disk/data/gempak/model/ecmwf/2005112300_ecmw
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : urel                                              
>> > SCALE : 0                                                 
>> >      051123/0000F000 :     -1.11616
>> > 
>> >  GEMPAK-GDPOINT>gfu = mul(urel,urel)
>> >  GEMPAK-GDPOINT>r
>> > GDFILE: /home/disk/data/gempak/model/ecmwf/2005112300_ecmw
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : mul(urel,urel)                                    
>> > SCALE : 0                                                 
>> >      051123/0000F000 :      3.76344    ***** WRONG!!!!
>> >  GEMPAK-GDPOINT>gfu = mul(urel,-1.11616)
>> >  GEMPAK-GDPOINT>r
>> >  
>> > GDFILE: /home/disk/data/gempak/model/ecmwf/2005112300_ecmw
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : mul(urel,-1.11616)                                
>> > SCALE : 0                                                 
>> >      051123/0000F000 :      1.24581
>> >  GEMPAK-GDPOINT>gfun = vrel
>> >  GEMPAK-GDPOINT>r  
>> >  
> > GDFILE: /home/disk/data/gempak/model/ecmwf/2005112300_ecmw
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : vrel                                              
>> > SCALE : 0                                                 
>> >      051123/0000F000 :      0.76872
>> > 
>> >  GEMPAK-GDPOINT>gfu = mag(wnd)
>> >  GEMPAK-GDPOINT>r
>> >  
>> > GDFILE: /home/disk/data/gempak/model/ecmwf/2005112300_ecmw
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : mag(wnd)                                          
>> > SCALE : 0                                                 
>> >      051123/0000F000 :      2.10129    ***** WRONG!!!! s/b 1.355
>> > 
>> >  GEMPAK-GDPOINT>gfu = dirn(wnd)
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : dirn(wnd)                                         
>> >      051123/0000F000 :    159.70622    ***** WRONG!!!! s/b 124.6 degrees
>> > 
>> > 
>> > And, lest you think this is confined to the ecmwf file, here's the
>> > same bug with the 1-degree GFS file
>> >  GEMPAK-GDPOINT>r
>> > GDFILE: $HDS/gfs/2005112312_gfs.gem                       
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : urel                                              
>> > SCALE : 0                                                 
>> >      051123/1200F000 :     -0.83057
>> >  GEMPAK-GDPOINT>gfu = mul(urel,-0.83057)
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : mul(urel,-0.83057)                                
>> >      051123/1200F000 :      0.68985
>> >  GEMPAK-GDPOINT>gfu = mul(urel,urel)
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : mul(urel,urel)                                    
>> >      051123/1200F000 :      0.85965    ***** WRONG!!!! s/b 0.689
>> > 
>> > 
>> > and with the eta104 grids:
>> >  GEMPAK-GDPOINT>gdfi = $HDS/eta/2005112312_eta104.gem
>> >  GEMPAK-GDPOINT>gfu = urel
>> >  GEMPAK-GDPOINT>r
>> > GDFILE: $HDS/eta/2005112312_eta104.gem                    
>> > GPOINT: @8.458;92.63                                      
>> > GVCORD: hagl                                              
>> > GLEVEL: 10                                                
>> > GFUNC : urel                                              
>> > SCALE : 0                                                 
>> >      051123/1200F000 :     -2.27911
>> >  GEMPAK-GDPOINT>gfu= mul(urel,-2.27911)
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : mul(urel,-2.27911)                                
>> >      051123/1200F000 :      5.19433
>> >  Parameters requested: GDATTIM,GDFILE,GLEVEL,GPOINT,GVCORD,GFUNC,SCALE.
>> >  GEMPAK-GDPOINT>gfu = mul(urel,urel)
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : mul(urel,urel)                                    
>> >      051123/1200F000 :      5.28800    ***** WRONG!!!! s/b 5.194
>> > 
>> >  GEMPAK-GDPOINT>gfun = vrel
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : vrel                                              
>> >      051123/1200F000 :      2.67976
>> > 
>> >  GEMPAK-GDPOINT>gfu = mag(vecr(urel,vrel))
>> >  GEMPAK-GDPOINT>r
>> > GFUNC : mag(vecr(urel,vrel))                              
>> >      051123/1200F000 :      3.58182    ***** WRONG!!!! s/b 3.518
>> >                                              
>> > ======================================================================
>> > 
>> > Thanks in advance for taking a look at this.
>> > 
>> > 
>> > David
>> > 
>> > -- 
>> > David Ovens                 e-mail: address@hidden
>> > Research Meteorologist    phone: (206) 685-8108
>> > Dept of Atm. Sciences      plan: Real-time MM5 forecasting for the
>> > Box 351640                        Pacific Northwest
>> > University of Washington          http://www.atmos.washington.edu/mm5rt
>> > Seattle, WA  98195               Weather Graphics and Loops
>> >                                   http://www.atmos.washington.edu/~ovens/l
> oops
>> 
>> -- 
>> David Ovens           e-mail: address@hidden
>> Research Meteorologist    phone: (206) 685-8108
>> Dept of Atm. Sciences      plan: Real-time MM5 forecasting for the
>> Box 351640                        Pacific Northwest
>> University of Washington          http://www.atmos.washington.edu/mm5rt
>> Seattle, WA  98195               Weather Graphics and Loops
>>                                   http://www.atmos.washington.edu/~ovens/loo
> ps
>
>-- 
>David Ovens             e-mail: address@hidden
>Research Meteorologist    phone: (206) 685-8108
>Dept of Atm. Sciences      plan: Real-time MM5 forecasting for the
>Box 351640                        Pacific Northwest
>University of Washington          http://www.atmos.washington.edu/mm5rt
>Seattle, WA  98195               Weather Graphics and Loops
>                                  http://www.atmos.washington.edu/~ovens/loops
>
--
NOTE: All email exchanges with Unidata User Support are recorded in the
Unidata inquiry tracking system and then made publicly available
through the web.  If you do not want to have your interactions made
available in this way, you must let us know in each email you send to us.