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

20021108: AWIPS 220 grid



Stonie,

After looking at the gbpolr.c code, and the OESA88 KWBM product,
I believe that:

1) The OESA88 KWBM product is being created with the wrong
   I & J scanning direction flags in byte 28 of the GDS.

   As seen in the NOAAPORT data stream, byte 28 is:

   64 = 01000000

     which in table 8 of the WMO388 grim manual provides
     that the I scanning mode=0 is the +i direction,

     and the J scanning mode=1 is the +j direction.

I believe that this should be:

  128 = 10000000


2) In order to make the grid decode corectly in dcgrib2 or
   nagrib (using cpyfil=gds), I was able to hack the gbpolr.c code
   to flip the sense of the I and J scanning mode flags when the
   projection plane is the south pole. I don't beleive that this
   is the correct approach...and would have advese consequences
   if any other S. Hemisphere polar stereographic grids exist
   (grid #28 comes to mind as a touchstone, but I don't know
   if any of these grids are around).



I will provide you with an attatchment of gbpolr.c that has a hemflag
check and flips the sense of I and J direction for the southern hemisphere.
You can try this at your own discretion. You should see that the OESA88 KWBM
grid is correctly plotted at this point....but no guarantees as to how long
this product will exist in this manner.

I would like to get someone at NCEP to comment of the creation of the AWIPS grid
#220. Since AWIPS generally hard codes the bounds in the software rather than
reading the values from the actual data products, they may not have encountered
the problem. If we can clarify the Byte 28 settings before committing to
code changes, it would be better in the long run.

Steve Chiswell

#include "gbcmn.h"

void gb_polr ( unsigned char *ptarray )
/************************************************************************
 * gb_polr                                                              *
 *                                                                      *
 * This routine will get the GDS information for a Polar Stereographic  *
 * grid. All of the information is put into the data structure for      *
 * the GDS.                                                             *
 *                                                                      *
 * gb_polr ( ptarray )                                                  *
 *                                                                      *
 * Input parameters:                                                    *
 *      *ptarray        unsigned char   Data buffer                     *
 **                                                                     *
 * Log:                                                                 *
 * J. Chou/EAI          02/93                                           *
 * S. Jacobs/EAI        10/93           Copied from GBUTIL              *
 * S. Jacobs/EAI         1/94           Clean up; Rename variables      *
 * L. Williams/EAI       7/94           Reformat header                 *
 * D.W.Plummer/NCEP      2/96           Cleanup GBDIAGs and comments    *
 * S. Jacobs/NCEP       12/00   Added prototypes                        *
 * S. Chiswell/Unidata  11/02   Apparanet problem with Awips grid 220   *
 *                              i/j scanning direction for S. Hemisphere*
 *                              appears to be against convention        *
 ***********************************************************************/
{
int     Dx, Dy, flag1, flag2, mode, Nx, Ny, La1, Lo1, LoV, dummy, indx;

int     hemflag;

double  loncnt, lat1, lon1, rtemp, X1, X2, Y1, Y2,
        TDx, TDy, Xll, Xur, Yll, Yur;

/*---------------------------------------------------------------------*/
        /*
         * BYTES 7-8
         * Nx - number of points along X-axis
         */
        indx = 6;
        Nx = gb_btoi(ptarray, indx, 2, FALSE );

        /*
         * BYTES 9-10
         * Ny - number of points along Y-axis
         */
        indx = 8;
        Ny = gb_btoi(ptarray, indx, 2, FALSE );

        /*
         * BYTES 11-13
         * La1 - latitude of first grid point
         */
        indx = 10;
        La1 = gb_btoi(ptarray, indx, 3, TRUE );
        lat1 = ( La1 / 1000.0 ) * DEG_TO_RAD;

        /*
         * BYTES 14-16
         * Lo1 - longitude of first grid point
         */
        indx = 13;
        Lo1 = gb_btoi(ptarray, indx, 3, TRUE );
        lon1 = ( Lo1 / 1000.0 ) * DEG_TO_RAD;

        /*
         * BYTE 17
         * Resolution and component flags
         */
        indx = 16;
        flag1 = gb_btoi(ptarray, indx, 1, FALSE );

        /*
         * BYTES 18-20
         * Lov - orientation of the grid (center longitude)
         */
        indx = 17;
        LoV = gb_btoi(ptarray, indx, 3, TRUE );
        loncnt = ( LoV / 1000.00 ) * DEG_TO_RAD;

        /*
         * BYTES 21-23
         * Dx - X-direction grid length
         */
        indx = 20;
        Dx = gb_btoi(ptarray, indx, 3, TRUE );

        /*
         * BYTES 24-26
         * Dy - Y-dirction grid length
         */
        indx = 23;
        Dy = gb_btoi(ptarray, indx, 3, TRUE );

        /*
         * BYTE 27
         * Projection center flag  [0=NP, 1=SP]
         */
        indx = 26;
        flag2 = gb_btoi(ptarray, indx, 1, FALSE );

        /*
         * BYTE 28
         * Scanning mode
         */
        indx = 27;
        mode = gb_btoi(ptarray, indx, 1, FALSE );

        /*
         * BYTES 29-32
         * Reserved - set to zero
         */
        indx = 28;
        dummy = gb_btoi(ptarray, indx, 4, FALSE );

        /*
         * Compute the linear coordinates of the first point.
         */
        if ( flag2 == 0 ) {
                /*
                 * Northern Hemisphere
                 */
                X1 =  RADIUS * tan ( PI4TH - lat1/2 ) * sin ( lon1-loncnt );
                Y1 = -RADIUS * tan ( PI4TH - lat1/2 ) * cos ( lon1-loncnt );
                hemflag = 1;
        }
        else {
                /*
                 * Southern Hemisphere
                 */
                X1 = RADIUS * tan ( PI4TH + lat1/2 ) * sin ( lon1-loncnt );
                Y1 = RADIUS * tan ( PI4TH + lat1/2 ) * cos ( lon1-loncnt );
                hemflag = -1;
        }

        /*
         * Compute the grid spacing
         */
        TDx = Dx / ( 1 + sin ( PI3RD ) );
        TDy = Dy / ( 1 + sin ( PI3RD ) );

        /*
         * Compute the linear coordinates of the second point.
         * Check the scanning mode to determine +/- grid spacing.
         */
        if ( ( mode >> 7 ) == 1 )
            X2 = X1 + ( Nx - 1 ) * ( -TDx ) * hemflag;
        else
            X2 = X1 + ( Nx - 1 ) * TDx * hemflag;

        if ( ( ( mode >> 6 ) & 1 ) == 1 )
            Y2 = Y1 + ( Ny - 1 ) * TDy * hemflag;
        else
            Y2 = Y1 + ( Ny - 1 ) * ( -TDy ) * hemflag;

        /*
         * Get the point coordinates paired as Lower Left and Upper Right.
         */
        Xll = G_MIN( X1, X2 );
        Yll = G_MIN( Y1, Y2 );
        Xur = G_MAX( X1, X2 );
        Yur = G_MAX( Y1, Y2 );

        /*
         * Compute the lat/lon of the first point.
         */
        if ( flag2 == 0 ) {
                /*
                 * Northern Hemisphere
                 */
                gds.latll = ( HALFPI - 2 *
                      atan2 ( sqrt ( pow (Xll,2.0) + pow (Yll,2.0) ),
                              RADIUS ) ) * RAD_TO_DEG;
                rtemp = loncnt + atan2 ( Xll, -Yll );
        }
        else {
                /*
                 * Southern Hemisphere
                 */
                gds.latll = - ( HALFPI - 2 *
                      atan2 ( sqrt ( pow (Xll,2.0) + pow (Yll,2.0) ),
                              RADIUS ) ) * RAD_TO_DEG;
                rtemp = loncnt + atan2 ( Xll,  Yll );
        }

        if ( rtemp > PI )
            gds.lonll = ( rtemp - TWOPI ) * RAD_TO_DEG;
        else if ( rtemp  < -PI )
            gds.lonll = ( rtemp + TWOPI ) * RAD_TO_DEG;
        else
            gds.lonll = rtemp * RAD_TO_DEG;

        /*
         * Compute the lat/lon of the second point.
         */
        if ( flag2 == 0 ) {
                /*
                 * Northern Hemisphere
                 */
                gds.latur = ( HALFPI - 2 *
                      atan2 ( sqrt ( pow (Xur,2.0) + pow (Yur,2.0) ),
                              RADIUS ) ) * RAD_TO_DEG;
                rtemp = loncnt + atan2 ( Xur, -Yur );
        }
        else {
                /*
                 * Southern Hemisphere
                 */
                gds.latur = - ( HALFPI - 2 *
                      atan2 ( sqrt ( pow (Xur,2.0) + pow (Yur,2.0) ),
                              RADIUS ) ) * RAD_TO_DEG;
                rtemp = loncnt + atan2 ( Xur,  Yur );
        }

        if ( rtemp > PI )
            gds.lonur = ( rtemp - TWOPI ) * RAD_TO_DEG;
        else if ( rtemp < -PI )
            gds.lonur = ( rtemp + TWOPI ) * RAD_TO_DEG;
        else
            gds.lonur = rtemp * RAD_TO_DEG;

                /*
                 * Set the proj angle information, the number of x,y points,
                 * and the coded flags.
                 */
        if ( flag2 == 0 ) {
                /*
                 * North Pole
                 */
                gds.angle1 = 90.0;
        }
        else {
                /*
                 * South Pole
                 */
                gds.angle1 = -90.0;
        }
        gds.angle2 = (float) ( LoV / 1000.0 );
        gds.angle3 = 0.0;

        gds.kx = Nx;
        gds.ky = Ny;

        gds.flag1 = flag1;
        gds.flag2 = flag2;
        gds.scan_mode = mode;

        if ( GBDIAG_GDS == TRUE )  {
            printf ( " GDS bytes  7 -  8 (Nx)            = %d\n", Nx );
            printf ( " GDS bytes  9 - 10 (Ny)            = %d\n", Ny );
            printf ( " GDS bytes 11 - 13 (La1)           = %d\n", La1 );
            printf ( " GDS bytes 14 - 16 (Lo1)           = %d\n", Lo1 );
            printf ( " GDS byte       17 (flag1)         = %d\n", flag1 );
            printf ( " GDS bytes 18 - 20 (LoV)           = %d\n", LoV );
            printf ( " GDS bytes 21 - 23 (Dx)            = %d\n", Dx );
            printf ( " GDS bytes 24 - 26 (Dy)            = %d\n", Dy );
            printf ( " GDS byte       27 (flag2)         = %d\n", flag2 );
            printf ( " GDS byte       28 (scan_mode)     = %d\n", mode );
            printf ( " GDS bytes 29 - 32 (skipped)\n" );
        }

}