external/io_grib1/WGRIB/gds_grid.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include "grib.h"
   4 #include "bds.h"
   5 #include "gds.h"
   6
   7 /*
   8  * get grid size from GDS
   9  *
  10  * added calculation of nxny of spectral data and clean up of triangular
  11  * grid nnxny calculation     l. kornblueh
  12  * 7/25/03 wind fix Dusan Jovic
  13  * 9/17/03 fix scan mode
  14  */
  15
  16 int GDS_grid(unsigned char *gds, unsigned char *bds, int *nx, int *ny,
  17              long int *nxny) {
  18
  19     int i, d, ix, iy, pl;
  20     long int isum;
  21
  22     *nx = ix = GDS_LatLon_nx(gds);
  23     *ny = iy = GDS_LatLon_ny(gds);
  24     *nxny = ix * iy;
  25
  26     /* thin grid */
  27
  28     if (GDS_Gaussian(gds) || GDS_LatLon(gds)) {
  29         if (ix == 65535) {
  30             *nx = -1;
  31             /* reduced grid */
  32             isum = 0;
  33             pl = GDS_PL(gds);
  34             for (i = 0; i < iy; i++) {
  35                 isum += gds[pl+i*2]*256 + gds[pl+i*2+1];
  36             }
  37             *nxny = isum;
  38         }
  39         return 0;
  40     }
  41     if (GDS_Triangular(gds)) {
  42         i = GDS_Triangular_ni(gds);
  43         d = GDS_Triangular_nd(gds);
  44         *nx = *nxny = d * (i + 1) * (i + 1);
  45         *ny = 1;
  46         return 0;
  47     }
  48     if (GDS_Harmonic(gds)) {
  49         /* this code assumes j, k, m are consistent with bds */
  50         *nx = *nxny = (8*(BDS_LEN(bds)-15)-BDS_UnusedBits(bds))/
  51                 BDS_NumBits(bds)+1;
  52            if ((8*(BDS_LEN(bds)-15)-BDS_UnusedBits(bds)) % BDS_NumBits(bds)) {
  53                fprintf(stderr,"inconsistent harmonic BDS\n");
  54            }
  55         *ny = 1;
  56     }
  57     return 0;
  58 }
  59
  60 #define NCOL 15
  61 void GDS_prt_thin_lon(unsigned char *gds) {
  62     int iy, i, col, pl;
  63
  64     iy = GDS_LatLon_ny(gds);
  65     iy = (iy + 1) / 2;
  66     iy = GDS_LatLon_ny(gds);
  67
  68     if ((pl = GDS_PL(gds)) == -1) {
  69         fprintf(stderr,"\nprogram error: GDS_prt_thin\n");
  70         return;
  71     }
  72     for (col = i = 0; i < iy; i++) {
  73         if (col == 0) printf("   ");
  74         printf("%5d", (gds[pl+i*2] << 8) + gds[pl+i*2+1]);
  75         col++;
  76         if (col == NCOL) {
  77             col = 0;
  78             printf("\n");
  79         }
  80     }
  81     if (col != 0) printf("\n");
  82 }
  83
  84 /*
  85  * prints out wind rel to grid or earth
  86  */
  87
  88 static char *scan_mode[8] = {
  89         "WE:NS",
  90         "NS:WE",
  91
  92         "WE:SN",
  93         "SN:WE",
  94
  95         "EW:NS",
  96         "NS:EW",
  97
  98         "EW:SN",
  99         "SN:EW" };
 100
 101
 102 void GDS_winds(unsigned char *gds, int verbose) {
 103     int scan = -1, mode = -1;
 104
 105     if (gds != NULL) {
 106         if (GDS_LatLon(gds)) {
 107             scan = GDS_LatLon_scan(gds);
 108             mode = GDS_LatLon_mode(gds);
 109         }
 110         else if (GDS_Mercator(gds)) {
 111             scan =GDS_Merc_scan(gds);
 112             mode =GDS_Merc_mode(gds);
 113         }
 114         /* else if (GDS_Gnomonic(gds)) { */
 115         else if (GDS_Lambert(gds)) {
 116             scan = GDS_Lambert_scan(gds);
 117             mode = GDS_Lambert_mode(gds);
 118         }
 119         else if (GDS_Gaussian(gds)) {
 120             scan = GDS_LatLon_scan(gds);
 121             mode = GDS_LatLon_mode(gds);
 122         }
 123         else if (GDS_Polar(gds)) {
 124             scan = GDS_Polar_scan(gds);
 125             mode = GDS_Polar_mode(gds);
 126         }
 127         else if (GDS_RotLL(gds)) {
 128             scan = GDS_RotLL_scan(gds);
 129             mode = GDS_RotLL_mode(gds);
 130         }
 131         /* else if (GDS_Triangular(gds)) { */
 132         else if (GDS_ssEgrid(gds)) {
 133             scan = GDS_ssEgrid_scan(gds);
 134             mode = GDS_ssEgrid_mode(gds);
 135         }
 136         else if (GDS_fEgrid(gds)) {
 137             scan = GDS_fEgrid_scan(gds);
 138             mode = GDS_fEgrid_mode(gds);
 139         }
 140         else if (GDS_ss2dEgrid(gds)) {
 141             scan = GDS_ss2dEgrid_scan(gds);
 142             mode = GDS_ss2dEgrid_mode(gds);
 143         }
 144     }
 145     if (verbose == 1) {
 146         if (mode != -1) {
 147             if (mode & 8) printf("winds in grid direction:");
 148             else printf("winds are N/S:");
 149         }
 150     }
 151     else if (verbose == 2) {
 152         if (scan != -1) {
 153             printf(" scan: %s", scan_mode[(scan >> 5) & 7]);
 154         }
 155         if (mode != -1) {
 156             if (mode & 8) printf(" winds(grid) ");
 157             else printf(" winds(N/S) ");
 158         }
 159     }
 160 }
 161
 162