/*****************************************************************************
 * writeshp.c
 *
 * DESCRIPTION
 *    This file contains all the routines used to write the grid out to
 * .shp format (which can be used via ArcGIS or ArcExplorer)
 * Associated with the .shp file are a .dbf, and .shx file.

 * .flt format (which can be used via ArcGIS S.A. or by GrADS.
 * Associated with the .flt file are a .prj, .hdr, and .ave files.
 * Also calls gribWriteGradsCTL, to create a .ctl file.
 *
 * HISTORY
 *   9/2002 Arthur Taylor (MDL / RSIS): Created.
 *  12/2002 Rici Yu, Fangyu Chi, Mark Armstrong, & Tim Boyer
 *          (RY,FC,MA,&TB): Code Review 2.
 *   6/2003 AAT: Split write.c into 2 pieces (writeflt and writeshp).
 *
 * NOTES
 * 1) As far as I can tell .flt and .shp files don't allow for two missing
 *    values, so I use missPri for both, if need be.
 * 2) Look up "address" in dbf definition.
 *****************************************************************************
 */
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "mymapf.h"
#include "tendian.h"
#include "write.h"
#include "myerror.h"
#include "myutil.h"
#include "myassert.h"
#include <errno.h>
#include "type.h"
#include "weather.h"
#include "chain.h"
#include <math.h>

extern double POWERS_ONE[];

/*****************************************************************************
 * checkFileSize() -- Review 12/2002
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This attempts to find out if a file is the expected size.
 *
 * ARGUMENTS
 * filename = Name of file to check the size of. (Input)
 *      len = Expected length. (Input)
 *
 * FILES/DATABASES:
 * Looks at the length of a given file.
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = Problems opening the files.
 * -2 = File is either wrong size or errors occured while checking the size.
 *
 * HISTORY
 *  11/2002 Arthur Taylor (MDL/RSIS): Created.
 *  12/2002 (RY,FC,MA,&TB): Code Review.
 *   5/2003 AAT: Removed reliance on errno (since Tcl/Tk confuses the issue).
 *
 * NOTES
 *****************************************************************************
 */
static int checkFileSize (char *filename, sInt4 len)
{
   FILE *fp;
   if ((fp = fopen (filename, "rb")) == NULL) {
      errSprintf ("ERROR: After write, Problems opening %s for file length "
                  "check.", filename);
      return -1;
   }
   fseek (fp, 0L, SEEK_END);
   if (ftell (fp) != len) {
      errSprintf ("ERROR: file %s is %ld not %ld bytes long.\n",
                  filename, ftell (fp), len);
      fclose (fp);
      return -2;
   }
   fclose (fp);
   return 0;
}

/*****************************************************************************
 * CreatePrj() --
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This creates a .prj file.  The .prj file contains the information that
 * ArcGIS 8 needs to understand the map projection of the .shp files.  Since I
 * have already projected the .shp files to a geographic coordinate system
 * (From a lambert conformal projection), the only thing that I think ESRI
 * needs to know is the assumed shape of the earth.
 *   This could eventually replace the .ave file, but as long as ArcView 3.*
 * is out, I'm going to create both.
 *
 * ARGUMENTS
 * filename = Name of file to save to. (Input/Output)
 *      gds = Grid Definition (from the parsed GRIB message) to write. (Input)
 *
 * FILES/DATABASES:
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = Problems opening the files.
 * -2 = Problems writing entire .shp file
 * -3 = Problems writing entire .shx file
 *
 * HISTORY
 *   1/2005 Arthur Taylor (MDL): Created.
 *
 * NOTES
 *****************************************************************************
 */
static int CreatePrj (char *filename, gdsType *gds)
{
   FILE *fp;

   strncpy (filename + strlen (filename) - 3, "prj", 3);
   if ((fp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }
   if (gds->f_sphere) {
      if( (gds->majEarth > 6370.0) && (gds->majEarth < 6380.0) ){
         fprintf (fp, "GEOGCS[\"NCEP_SPHERE\",DATUM[\"NCEP_SPHERE\",SPHEROID"
                  "[\"NCEP_SPHERE\",6371200.0,0.0]],PRIMEM[\"Greenwich\",0.0]"
                  ",UNIT[\"Degree\",0.0174532925199432955]]\n");
      } else {
         fprintf (fp, "GEOGCS[\"NCEP_SPHERE\",DATUM[\"NCEP_SPHERE\",SPHEROID"
                  "[\"NCEP_SPHERE\",%f,0.0]],PRIMEM[\"Greenwich\",0.0],UNIT["
                  "\"Degree\",0.0174532925199432955]]\n", gds->majEarth);
      }
   } else {
      fprintf (fp, "GEOGCS[\"NCEP_SPHERE\",DATUM[\"NCEP_SPHERE\",SPHEROID"
               "[\"NCEP_SPHERE\",%f,%f]],PRIMEM[\"Greenwich\",0.0],UNIT["
               "\"Degree\",0.0174532925199432955]]\n", gds->majEarth,
               gds->minEarth);
   }
   fclose (fp);
   return 0;
}

/*****************************************************************************
 * CreateShpPoly() --
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This creates a .shp / .shx file.  The .shp / .shx file contains the
 * lat/lon values of the grid as polygons instead of as points.  These formats
 * are specific to Esri ArcView.
 *   The points are assumed to be the corner points of the grid, so there
 * should be 1 more row/column in the dp than given to the shpPnt command,
 * but Nx and Ny are the same as with shpPnt.
 *
 * ARGUMENTS
 *   filename = Name of file to save to. (Output)
 *         dp = Array of lat/lon pairs for the grid cells. (Input)
 *              dp is expected to have 1 extra set of x than Nx
 *              and 1 extra set of y than Ny.
 *         Nx = Number of x values in the grid. (Input)
 *         Ny = Number of y values in the grid. (Input)
 * f_nMissing = True if we should NOT store missing. (Input)
 *  grib_Data = Actual data so we can determine where missing values are. (In)
 *     attrib = Tells what type of missing values we used. (Input)
 *
 * FILES/DATABASES:
 *   Creates a .shp file, which is a binary file (Mixed Endian) consisting of
 *      a set of lat/lon polygon shapes.
 *   Also, creates a .shx file, which is a binary file (Mixed Endian) which is
 *      used as an index into the .shp file.
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = Problems opening the files.
 * -2 = Problems writing entire .shp file
 * -3 = Problems writing entire .shx file
 *
 * HISTORY
 *   4/2003 Arthur Taylor (MDL/RSIS): Created.
 *   5/2003 AAT: Modified to allow for f_nMissing.
 *   5/2003 AAT: Fixed orientation to polygon issue.
 *   5/2003 AAT: For poly spanning date line, allow outside of -180..180
 *   5/2003 AAT: Decided to have 1,1 be lower left corner in .shp files.
 *   5/2003 AAT: Removed reliance on errno (since Tcl/Tk confuses the issue).
 *   7/2003 AAT: 1,1 lower left affected orientation of polygons.
 *   3/2004 AAT: Updated to handle Alaska polygons.
 *
 * NOTES
 * 1) Doesn't inter-twine the write to files, as that could cause the hard
 *   drive to slow down.
 * 2) Assumes data goes from left to right before going up and down.
 *****************************************************************************
 */
static int CreateShpPoly (char *filename, LatLon *dp, int Nx, int Ny,
                          sChar f_nMissing, double *grib_Data,
                          gridAttribType *attrib, sChar LatLon_Decimal)
{
   FILE *sfp;           /* The open file pointer for .shp */
   FILE *xfp;           /* The open file pointer for .shx. */
   int i;               /* A counter used for the shx values. */
   sInt4 Head1[7];      /* The Big endian part of the Header. */
   sInt4 Head2[2];      /* The Little endian part of the Header. */
   double Bounds[] = {
      0., 0., 0., 0., 0., 0., 0., 0.
   };                   /* Spatial bounds of the data. minLon, minLat,
                         * maxLon, maxLat, ... */
   sInt4 dataType = 5;  /* Polygon shp type. */
   sInt4 curRec[2];     /* rec number, and content length. */
   LatLon *curDp;       /* current data point. */
   double *curData;     /* Current Grid cell data point (for f_nMissing) */
   int err;             /* Internal err number. */
   int recLen;          /* Length in bytes of a record in the .shp file. */
   int recLen2;         /* Length in bytes of a 2 poly record. */
   sInt4 dpLen;         /* Number of data polygons in output file. */
   sInt4 numRec;        /* The total number of records actually stored. */
   int x, y;            /* Counters used for traversing the grid. */
   double pts[10];      /* Holds the lon/lat points for the cur polygon. */
   double pts1[20];     /* The points for cur poly if split on dateline. */
   double pts2[20];     /* The points for cur poly if split on dateline. */
   int k;               /* Used for traversing the polygon record. */
   double *cur;         /* Used for traversing the polygon record. */
   double PolyBound[4]; /* Used for holding the polygon bounds. */
   sInt4 PolygonSpecs[] = {
      1, 5, 0
   };                   /* NumParts, NumPnts, Index of Ring are constant for
                         * This type of .shp polygon. */
   sChar *f_dateline;   /* array of flags if the poly crosses the dateline. */
   int indexLf = 0;     /* Where to add data to the left polygon chain. */
   int indexRt = 0;     /* Where to add data to the right polygon chain. */
   sChar f_left;        /* Flag to add to the left or right polygon */
   double delt;         /* change in lon of a given line segment. */
   double newLat;       /* latitude where line segment crosses the dateline */
   sInt4 secChain;      /* first index to the second chain. */
   int dateLineCnt;     /* A count of number of polys that cross dateline. */

   strncpy (filename + strlen (filename) - 3, "shp", 3);
   if ((sfp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }
   dpLen = Nx * Ny;
   /* Perform a simple file size check.  File has to be less than
    * 4,294,967,294 bytes. (2^31-1) * 2.  Reasoning: file size in 2 byte words
    * is stored as a sInt4 */
   if (!f_nMissing) {
      /* If include all cells (assuming no dateline issues), the size of the
       * .shp file is: 100 + dpLen * (8 +4 +4*8 +3*4 +10*8= 180bytes)
       * So dpLen <= (2^31-1) * 2 - 100 / 180 = 23,860,928.8 */
      if (dpLen > 23860928L) {
         errSprintf ("Trying to create a small poly shp file with %d cells.\n"
                     "This is > small polygon maximum of 23,860,928\n",
                     dpLen);
         return -1;
      }
   }

   /* Start Writing header in first 100 bytes. */
   Head1[0] = 9994;     /* ArcView identifier. */
   memset ((Head1 + 1), 0, 5 * sizeof (sInt4)); /* set 5 unused to 0 */
   recLen = sizeof (sInt4) + 4 * sizeof (double) + 3 * sizeof (sInt4) +
         10 * sizeof (double);
   recLen2 = sizeof (sInt4) + 4 * sizeof (double) +
         2 * sizeof (sInt4) + 2 * sizeof (sInt4) + 10 * sizeof (double) +
         10 * sizeof (double);
   /* .shp file size (in 2 byte words). */
   Head1[6] = (100 + (2 * sizeof (sInt4) + recLen) * dpLen) / 2;
   FWRITE_BIG (Head1, sizeof (sInt4), 7, sfp);
   Head2[0] = 1000;     /* ArcView version identifier. */
   Head2[1] = dataType; /* Signal that these are polygon data. */
   FWRITE_LIT (Head2, sizeof (sInt4), 2, sfp);
   /* Write out initial guess for bounds... Need to revisit. */
   FWRITE_LIT (Bounds, sizeof (double), 8, sfp);

   /* Start Writing data. */
   curRec[1] = recLen / 2; /* Content length in (2 byte words) */
   curRec[0] = 1;
   dateLineCnt = 0;
   f_dateline = (sChar *) malloc ((Nx * Ny) * sizeof (sChar));

   /* If the simple case where we don't have to worry about removing the
    * missing values.  */
   for (y = 0; y < Ny; y++) {
      curData = grib_Data + y * Nx;
      curDp = dp + y * (Nx + 1);
      for (x = 0; x < Nx; x++) {
         if ((!f_nMissing) || (attrib->f_miss == 0) ||
             ((*curData != attrib->missPri) &&
              ((attrib->f_miss != 2) || (*curData != attrib->missSec)))) {
            /* Get the current polygon. */
            cur = pts;
            /* Order matters here.  Must be clockwise !!! */
            *(cur++) = curDp->lon;
            *(cur++) = curDp->lat;
            /* Switched again on 7/14/2003 because of lower left adjustment
             * in 5/2003. */
            *(cur++) = curDp[Nx + 1].lon; /* 1 row up. */
            *(cur++) = curDp[Nx + 1].lat; /* 1 row up. */
            *(cur++) = curDp[Nx + 2].lon; /* 1 row up + 1 accross. */
            *(cur++) = curDp[Nx + 2].lat; /* 1 row up + 1 accross. */
            *(cur++) = curDp[1].lon;
            *(cur++) = curDp[1].lat;
            *(cur++) = curDp->lon;
            *cur = curDp->lat;

            /* Compute the bounds of this polygon. */
            cur = pts;
            PolyBound[0] = *cur; /* min x */
            PolyBound[2] = *cur; /* max x */
            cur++;
            PolyBound[1] = *cur; /* min y */
            PolyBound[3] = *cur; /* max y */
            cur++;
            for (k = 1; k <= 3; k++) {
               if (PolyBound[0] > *cur)
                  PolyBound[0] = *cur;
               else if (PolyBound[2] < *cur)
                  PolyBound[2] = *cur;
               cur++;
               if (PolyBound[1] > *cur)
                  PolyBound[1] = *cur;
               else if (PolyBound[3] < *cur)
                  PolyBound[3] = *cur;
               cur++;
            }
            f_dateline[curRec[0] - 1] = 0;
            if ((PolyBound[2] - PolyBound[0]) > 180) {
               f_dateline[curRec[0] - 1] = 1;
               dateLineCnt++;
               PolyBound[2] = 180;
               PolyBound[0] = -180;
               f_left = 1;
               indexLf = 0;
               indexRt = 0;
               pts1[indexLf++] = pts[0];
               pts1[indexLf++] = pts[1];
               for (k = 1; k <= 4; k++) {
                  delt = pts[k * 2] - pts[(k - 1) * 2];
                  if ((delt > 180) || (delt < -180)) {
                     DateLineLat (pts[k * 2], pts[k * 2 + 1],
                                  pts[(k - 1) * 2], pts[(k - 1) * 2 + 1],
                                  &newLat);
                     newLat = myRound (newLat, LatLon_Decimal);
                     if (f_left) {
                        if (pts1[0] < 0) {
                           pts1[indexLf++] = -180;
                           pts2[indexRt++] = 180;
                        } else {
                           pts1[indexLf++] = 180;
                           pts2[indexRt++] = -180;
                        }
                        pts1[indexLf++] = newLat;
                        pts2[indexRt++] = newLat;
                        pts2[indexRt++] = pts[k * 2];
                        pts2[indexRt++] = pts[k * 2 + 1];
                        f_left = 0;
                     } else {
                        if (pts1[0] < 0) {
                           pts1[indexLf++] = -180;
                           pts2[indexRt++] = 180;
                        } else {
                           pts1[indexLf++] = 180;
                           pts2[indexRt++] = -180;
                        }
                        pts2[indexRt++] = newLat;
                        pts1[indexLf++] = newLat;
                        pts1[indexLf++] = pts[k * 2];
                        pts1[indexLf++] = pts[k * 2 + 1];
                        f_left = 1;
                     }
                  } else {
                     if (f_left) {
                        pts1[indexLf++] = pts[k * 2];
                        pts1[indexLf++] = pts[k * 2 + 1];
                     } else {
                        pts2[indexRt++] = pts[k * 2];
                        pts2[indexRt++] = pts[k * 2 + 1];
                     }
                  }
               }
               pts2[indexRt++] = pts2[0];
               pts2[indexRt++] = pts2[1];
#ifdef DEBUG
/*
               for (k = 0; k < indexLf; k++) {
                  printf ("%f,", pts1[k]);
               }
               printf ("\n");
               for (k = 0; k < indexRt; k++) {
                  printf ("%f,", pts2[k]);
               }
               printf ("\n");
*/
#endif
            }

            /* Update Bounds of all data. */
            if (curRec[0] == 1) {
               Bounds[0] = PolyBound[0];
               Bounds[1] = PolyBound[1];
               Bounds[2] = PolyBound[2];
               Bounds[3] = PolyBound[3];
            } else {
               if (Bounds[0] > PolyBound[0])
                  Bounds[0] = PolyBound[0];
               if (Bounds[1] > PolyBound[1])
                  Bounds[1] = PolyBound[1];
               if (Bounds[2] < PolyBound[2])
                  Bounds[2] = PolyBound[2];
               if (Bounds[3] < PolyBound[3])
                  Bounds[3] = PolyBound[3];
            }
            if (f_dateline[curRec[0] - 1]) {
               /* Write record header. */
               curRec[1] = recLen2 / 2; /* Content length in (2 byte words) */
               FWRITE_BIG (curRec, sizeof (sInt4), 2, sfp);
               curRec[1] = recLen / 2; /* Content length in (2 byte words) */
               /* Write the data type. */
               FWRITE_LIT (&dataType, sizeof (sInt4), 1, sfp);
               /* Write polygons bounds */
               FWRITE_LIT (PolyBound, sizeof (double), 4, sfp);
               /* Write out the Polygon Specs. */
               PolygonSpecs[0] = 2;
               PolygonSpecs[1] = 10;
               FWRITE_LIT (PolygonSpecs, sizeof (sInt4), 3, sfp);
               PolygonSpecs[0] = 1;
               PolygonSpecs[1] = 5;
               /* Write out index of second chain. */
               secChain = indexLf / 2;
               FWRITE_LIT (&secChain, sizeof (sInt4), 1, sfp);
               /* Points ... 10 of them indexLf + indexRt = 20 (10 points) */
               FWRITE_LIT (pts1, sizeof (double), indexLf, sfp);
               FWRITE_LIT (pts2, sizeof (double), indexRt, sfp);
            } else {
               /* Write record header. */
               FWRITE_BIG (curRec, sizeof (sInt4), 2, sfp);
               /* Write the data type. */
               FWRITE_LIT (&dataType, sizeof (sInt4), 1, sfp);
               /* Write polygons bounds */
               FWRITE_LIT (PolyBound, sizeof (double), 4, sfp);
               /* Write out the Polygon Specs. */
               FWRITE_LIT (PolygonSpecs, sizeof (sInt4), 3, sfp);
               /* Points ... 5 of them */
               FWRITE_LIT (pts, sizeof (double), 10, sfp);
            }
            curRec[0]++;
            /* Assuming no dateline issues, the size of the .shp file is:
             * 100 + dpLen * (8 +4 +4*8 +3*4 +10*8= 180bytes)
             * So the max number of records allowed is <=
             * (2^31-1) * 2 - 100 / 180 = 23,860,928.8 */
            if (curRec[0] > 23860928) {
               errSprintf ("Trying to create a small poly shp file with %d cells.\n"
                           "This is > small polygon maximum of 23,860,928\n",
                           curRec[0]);
               return -1;
            }
         }
         curData++;
         curDp++;
      }
   }
   numRec = curRec[0] - 1;

   /* The dateline polys are bigger than the normal poly's by: 10 * sizeof
    * (double) + 1 * sizeof (sInt4); */

   /* Store the updated file length. */
   /* .shp file size (in 2 byte words). */
   Head1[6] = (100 + (2 * sizeof (sInt4) + recLen) * numRec) / 2;

   /* The dateline polys are bigger than the normal poly's by: 10 * sizeof
    * (double) + 1 * sizeof (sInt4); so we have to add that. */
   Head1[6] += (recLen2 - recLen) * dateLineCnt / 2;

   fseek (sfp, 24, SEEK_SET);
   FWRITE_BIG (&(Head1[6]), sizeof (sInt4), 1, sfp);

   /* Store the updated Bounds. */
   fseek (sfp, 36, SEEK_SET);
   /* FWRITE use 4 since we are only updating 4 bounds (not 8). */
   FWRITE_LIT (Bounds, sizeof (double), 4, sfp);
   fflush (sfp);
   fclose (sfp);

   /* Check that .shp is now the correct file size. */
   if ((err = checkFileSize (filename, Head1[6] * 2)) != 0) {
      free (f_dateline);
      return err;
   }

   filename[strlen (filename) - 1] = 'x';
   if ((xfp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      free (f_dateline);
      return -1;
   }

   /* Write ArcView header. */
   Head1[6] = (100 + 8 * numRec) / 2; /* shx file size (in words). */
   FWRITE_BIG (Head1, sizeof (sInt4), 7, xfp);
   FWRITE_LIT (Head2, sizeof (sInt4), 2, xfp);
   FWRITE_LIT (Bounds, sizeof (double), 8, xfp);

   curRec[0] = 50;      /* 100 bytes / 2 = 50 words */
   for (i = 0; i < numRec; i++) {
      if (f_dateline[i]) {
         curRec[1] = recLen2 / 2; /* Content length in words (2 bytes) */
      } else {
         curRec[1] = recLen / 2; /* Content length in words (2 bytes) */
      }
      FWRITE_BIG (curRec, sizeof (sInt4), 2, xfp);
      if (f_dateline[i]) {
         /* X / 2 because of (2 byte words) */
         curRec[0] += (recLen2 + 2 * sizeof (sInt4)) / 2;
      } else {
         /* X / 2 because of (2 byte words) */
         curRec[0] += (recLen + 2 * sizeof (sInt4)) / 2;
      }
   }
   fflush (xfp);
   fclose (xfp);

   /* Check that .shx is now the correct file size. */
   free (f_dateline);
   return checkFileSize (filename, 100 + 8 * numRec);
}

/*****************************************************************************
 * CreateShpPnt() -- Review 12/2002
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This creates a .shp / .shx file.  The .shp / .shx file contains the
 * lat/lon values of the grid as points instead of as polygons.  These formats
 * are specific to Esri ArcView.
 *
 * ARGUMENTS
 *       filename = Name of file to save to. (Output)
 *            map = Holds the current map projection info to compute from. (In)
 *            gds = Grid Definiton from the parsed GRIB msg to write. (Input)
 * LatLon_Decimal = Number of decimals to round lat/lon's to. (Input)
 *     f_nMissing = True if we should NOT store missing. (Input)
 *      grib_Data = Actual data to determine where missing values are. (In)
 *         attrib = Tells what type of missing values we used. (Input)
 *
 * FILES/DATABASES:
 *   Creates a .shp file, which is a binary file (Mixed Endian) consisting of
 *      a set of lat/lon point shapes.
 *   Also, creates a .shx file, which is a binary file (Mixed Endian) which is
 *      used as an index into the .shp file.
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = Problems opening the files.
 * -2 = Problems writing entire .shp file
 * -3 = Problems writing entire .shx file
 *
 * HISTORY
 *   9/2002 Arthur Taylor (MDL/RSIS): Created.
 *  11/2002 AAT: Noticed on HP that it didn't write entire .shp file
 *              (due to nfs mount timeouts).. Added a number of error checks.
 *  12/2002 (RY,FC,MA,&TB): Code Review.
 *   5/2003 AAT: Modified to allow us to not store missing values.
 *   5/2003 AAT: Decided to have 1,1 be lower left corner in .shp files.
 *   5/2003 AAT: Removed reliance on errno (since Tcl/Tk confuses the issue).
 *   8/2007 AAT: Modified to internally compute dp as needed, so we don't
 *               have to store massive quantities of lat/lons.
 *
 * NOTES
 * Doesn't inter-twine the write to files, as that could cause the hard
 *   drive to slow down.
 *****************************************************************************
 */
static int CreateShpPnt (char *filename, myMaparam *map, gdsType *gds,
                         sChar LatLon_Decimal, sChar f_nMissing,
                         double *grib_Data, gridAttribType *attrib)
{
   FILE *sfp;           /* The open file pointer for .shp */
   FILE *xfp;           /* The open file pointer for .shx. */
   int i;               /* A counter used for the shx values. */
   int x, y;            /* Counters used for traversing the grid. */
   sInt4 Head1[7];      /* The Big endian part of the Header. */
   sInt4 Head2[2];      /* The Little endian part of the Header. */
   double Bounds[] = {
      0., 0., 0., 0., 0., 0., 0., 0.
   };                   /* Spatial bounds of the data. minLon, minLat,
                         * maxLon, maxLat, ... */
   sInt4 dataType = 1;  /* Point shp type data. */
   sInt4 curRec[2];     /* rec number, and content length. */
   LatLon curDp;        /* Current lat/lon data point. */
   double *curData;     /* Current Grid cell data point (for f_nMissing) */
   int err;             /* Internal err number. */
   int recLen;          /* Length in bytes of a record in the .shp file. */
   sInt4 numRec;        /* The total number of records actually stored. */

   strncpy (filename + strlen (filename) - 3, "shp", 3);
   if ((sfp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }
   /* Perform a simple file size check.  File has to be less than
    * 4,294,967,294 bytes. (2^31-1) * 2.  Reasoning: file size in 2 byte words
    * is stored as a sInt4 */
   if ((!f_nMissing) || (attrib->f_miss == 0)) {
      /* If include all cells, the size of the
       * .shp file is: 100 + dpLen * (8 +4 +2*8= 28bytes)
       * So gds->Ny * gds->Nx <= (2^31-1) * 2 - 100 / 28 = 153,391,685.5 */
      if (gds->Ny * gds->Nx > 153391685) {
         errSprintf ("Trying to create a point shp file with %d cells.\n"
                     "This is > point maximum of 153,391,685\n",
                     gds->Ny * gds->Nx);
         return -1;
      }
   }

   /* Start Writing header in first 100 bytes. */
   Head1[0] = 9994;     /* ArcView identifier. */
   memset ((Head1 + 1), 0, 5 * sizeof (sInt4)); /* set 5 unused to 0 */
   recLen = sizeof (sInt4) + 2 * sizeof (double);
   /* .shp file size (in 2 byte words). */
   Head1[6] = (100 + (2 * sizeof (sInt4) + recLen) * gds->Nx * gds->Ny) / 2;
   FWRITE_BIG (Head1, sizeof (sInt4), 7, sfp);
   Head2[0] = 1000;     /* ArcView version identifier. */
   Head2[1] = dataType; /* Signal that these are point data. */
   FWRITE_LIT (Head2, sizeof (sInt4), 2, sfp);
   /* Write out initial guess for bounds... Need to revisit. */
   FWRITE_LIT (Bounds, sizeof (double), 8, sfp);

   /* Start Writing data. */
   curRec[1] = recLen / 2; /* Content length in (2 byte words) */
/*   curDp = dp;*/
   if ((!f_nMissing) || (attrib->f_miss == 0)) {
      myCxy2ll (map, 1, 1, &(curDp.lat), &(curDp.lon));
      curDp.lat = myRound (curDp.lat, LatLon_Decimal);
      curDp.lon = myRound (curDp.lon, LatLon_Decimal);
      Bounds[0] = Bounds[2] = curDp.lon;
      Bounds[1] = Bounds[3] = curDp.lat;
      for (y = 0; y < gds->Ny; y++) {
         for (x = 0; x < gds->Nx; x++) {
            curRec[0] = 1 + x + y * gds->Nx;
            myCxy2ll (map, x + 1, y + 1, &(curDp.lat), &(curDp.lon));
            curDp.lat = myRound (curDp.lat, LatLon_Decimal);
            curDp.lon = myRound (curDp.lon, LatLon_Decimal);

            FWRITE_BIG (curRec, sizeof (sInt4), 2, sfp);
            FWRITE_LIT (&dataType, sizeof (sInt4), 1, sfp);
            FWRITE_LIT (&(curDp.lon), sizeof (double), 1, sfp);
            FWRITE_LIT (&(curDp.lat), sizeof (double), 1, sfp);
            /* Update Bounds. */
            if (curDp.lon < Bounds[0]) {
               Bounds[0] = curDp.lon;
            } else if (curDp.lon > Bounds[2]) {
               Bounds[2] = curDp.lon;
            }
            if (curDp.lat < Bounds[1]) {
               Bounds[1] = curDp.lat;
            } else if (curDp.lat > Bounds[3]) {
               Bounds[3] = curDp.lat;
            }
         }
      }
      numRec = gds->Nx * gds->Ny;
   } else {
      curRec[0] = 1;
      for (y = 0; y < gds->Ny; y++) {
         curData = grib_Data + y * gds->Nx;
         for (x = 0; x < gds->Nx; x++) {
            if ((*curData != attrib->missPri) &&
                ((attrib->f_miss != 2) || (*curData != attrib->missSec))) {
               myCxy2ll (map, x + 1, y + 1, &(curDp.lat), &(curDp.lon));
               curDp.lat = myRound (curDp.lat, LatLon_Decimal);
               curDp.lon = myRound (curDp.lon, LatLon_Decimal);

               FWRITE_BIG (curRec, sizeof (sInt4), 2, sfp);
               FWRITE_LIT (&dataType, sizeof (sInt4), 1, sfp);
               FWRITE_LIT (&(curDp.lon), sizeof (double), 1, sfp);
               FWRITE_LIT (&(curDp.lat), sizeof (double), 1, sfp);
               /* Update Bounds. */
               if (curRec[0] == 1) {
                  Bounds[0] = Bounds[2] = curDp.lon;
                  Bounds[1] = Bounds[3] = curDp.lat;
               } else {
                  if (curDp.lon < Bounds[0]) {
                     Bounds[0] = curDp.lon;
                  } else if (curDp.lon > Bounds[2]) {
                     Bounds[2] = curDp.lon;
                  }
                  if (curDp.lat < Bounds[1]) {
                     Bounds[1] = curDp.lat;
                  } else if (curDp.lat > Bounds[3]) {
                     Bounds[3] = curDp.lat;
                  }
               }
               curRec[0]++;
               /* The size of the .shp file is:
                * 100 + dpLen * (8 +4 +2*8= 28bytes)
                * So curRec[0] <= (2^31-1) * 2 - 100 / 28 = 153,391,685.5 */
               if (curRec[0] > 153391685) {
                  errSprintf ("Trying to create a point shp file with %d cells.\n"
                              "This is > point maximum of 153,391,685\n",
                              curRec[0]);
                  return -1;
               }
            }
            curData++;
/*            curDp++;*/
         }
      }
      numRec = curRec[0] - 1;
      /* Store the updated file length. */
      /* .shp file size (in 2 byte words). */
      Head1[6] = (100 + (2 * sizeof (sInt4) + recLen) * numRec) / 2;
      fseek (sfp, 24, SEEK_SET);
      FWRITE_BIG (&(Head1[6]), sizeof (sInt4), 1, sfp);
   }
   /* Store the updated Bounds. */
   fseek (sfp, 36, SEEK_SET);
   /* FWRITE - 4 since we are only updating 4 bounds (not 8). */
   FWRITE_LIT (Bounds, sizeof (double), 4, sfp);
   fflush (sfp);
   fclose (sfp);
   /* Check that .shp is now the correct file size. */
   if ((err = checkFileSize (filename, Head1[6] * 2)) != 0) {
      return err;
   }

   filename[strlen (filename) - 1] = 'x';
   if ((xfp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }

   /* Write ArcView header (.shx file). */
   Head1[6] = (100 + 8 * numRec) / 2; /* shx file size (in words). */
   FWRITE_BIG (Head1, sizeof (sInt4), 7, xfp);
   FWRITE_LIT (Head2, sizeof (sInt4), 2, xfp);
   FWRITE_LIT (Bounds, sizeof (double), 8, xfp);

   curRec[0] = 50;      /* 100 bytes / 2 = 50 words */
   curRec[1] = recLen / 2; /* Content length in words (2 bytes) */
   for (i = 0; i < numRec; i++) {
      FWRITE_BIG (curRec, sizeof (sInt4), 2, xfp);
      curRec[0] += (recLen + 2 * sizeof (sInt4)) / 2; /* (2 byte words) */
   }
   fflush (xfp);
   fclose (xfp);

   /* Check that .shx is now the correct file size. */
   return checkFileSize (filename, 100 + 8 * numRec);
}

/*****************************************************************************
 * CreateDbf() -- Review 12/2002
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This creates a .dbf file.  The .dbf file contains the actual info, which
 * is associated with the points.  The .dbf file format is very flexible but
 * unfortunately is mostly in ASCII.
 *
 * ARGUMENTS
 *   filename = Name of file to save to. (Output)
 *     Nx, Ny = The dimensions of the grid. (Input)
 *  grib_Data = The actual data (assumed parsed by ParseGrid) (Input)
 *    element = Name of current variable. (Input)
 *   fieldLen = How many char to store the # in (includes decLen). (Input)
 *     decLen = How many char to store decimal part of # in. (Input)
 * f_nMissing = True if we should NOT store missing. (Input)
 *     attrib = What kind of missing value management is used. (Input)
 *  f_verbose = True if we want the verbose form. (Input)
 *        map = Holds the current map projection info to compute from. (In)
 *
 * FILES/DATABASES:
 *   Creates a .dbf file, which is a binary file (Little Endian) consisting of
 *      the info to associate with the lat/lon point shapes stored in the .shp
 *      and .shx files.
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = Problems opening the files.
 *
 * HISTORY
 *   9/2002 Arthur Taylor (MDL/RSIS): Created.
 *  12/2002 (RY,FC,MA,&TB): Code Review.
 *   4/2003 AAT: Removed "include lat/lon option", because it was not being
 *          used and because it is difficult to implement for polygon .shp
 *          files.
 *   5/2003 AAT Switched from using " %08ld%10.3f" to having the fieldLen
 *          and decLen passed in because GRIB1 default missing is so large
 *          that it couldn't be stored.
 *   5/2003 AAT Added f_nMissing option.
 *   5/2003 AAT: Added rounding to decimal.
 *   5/2003 AAT: Decided to have 1,1 be lower left corner in .shp files.
 *   1/2005 AAT: Modified for verbose output
 *
 * NOTES
 * 1) Look up "address" in dbf definition.
 *****************************************************************************
 */
static int CreateDbf (char *filename, sInt4 Nx, sInt4 Ny,
                      double *grib_Data, const char *element, uChar fieldLen,
                      sChar decLen, sChar f_nMissing,
                      gridAttribType *attrib, char f_verbose, myMaparam *map)
{
   uChar header[] = { 3, 101, 4, 20 }; /* Header info for dbf. */
   int numCol;          /* The number of columns we use. */
   char names[6][11];   /* Field (column) names. */
   char type[6];        /* Type of data (number in all columns.) */
   uChar fldLen[6];     /* field len for columns. */
   uChar fldDec[6];     /* field decimal lens for columns. */
   sInt4 reserved[] = { 0, 0, 0, 0, 0 }; /* need 20 bytes of 0. */
   sInt4 address = 0;   /* Address to find data? */
   sInt4 NxNy = Nx * Ny; /* The total number of values. */
   FILE *fp;            /* The open .dbf file pointer. */
   int x, y;            /* Current grid cell location. */
   short int recLen;    /* Size in bytes of one cell of data */
   sInt4 id;            /* A unique decimal id for the current cell. */
   double *curData;     /* A pointer to Grib data for current cell. */
   sInt4 totSize;       /* Total size of the .dbf file. */
   uChar uc_temp;       /* Temp. storage of type unsigned char. */
   short int si_temp;   /* Temp. storage of type short int. */
   char formatBuf[20];  /* Used to format float output sent to .dbf */
   sInt4 numRec;        /* The total number of records actually stored. */
   double shift;        /* power of 10 used in rounding. */
   double lat;          /* Latitude of the current point for f_verbose */
   double lon;          /* Longitude of the current point for f_verbose */

   if (decLen > 17)
      decLen = 17;
   if (decLen < 0)
      decLen = 0;
   shift = POWERS_ONE[decLen];
   strncpy (filename + strlen (filename) - 3, "dbf", 3);
   if ((fp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }
   if (f_verbose) {
      numCol = 6;
      strncpy (names[0], "POINTID", 11);
      type[0] = 'N';
      fldLen[0] = 8;
      fldDec[0] = 0;
      strncpy (names[1], "X", 11);
      type[1] = 'N';
      fldLen[1] = 4;
      fldDec[1] = 0;
      strncpy (names[2], "Y", 11);
      type[2] = 'N';
      fldLen[2] = 4;
      fldDec[2] = 0;
      strncpy (names[3], "LON", 11);
      type[3] = 'N';
      fldLen[3] = 10;
      fldDec[3] = 5;
      strncpy (names[4], "LAT", 11);
      type[4] = 'N';
      fldLen[4] = 9;
      fldDec[4] = 5;
      strncpy (names[5], element, 11);
      type[5] = 'N';
      fldLen[5] = fieldLen;
      fldDec[5] = decLen;
      recLen = (short int) (1 + fldLen[0] + fldLen[1] + fldLen[2] + fldLen[3]
                            + fldLen[4] + fldLen[5]);
   } else {
      numCol = 2;
      strncpy (names[0], "POINTID", 11);
      type[0] = 'N';
      fldLen[0] = 8;
      fldDec[0] = 0;
      strncpy (names[1], element, 11);
      type[1] = 'N';
      fldLen[1] = fieldLen;
      fldDec[1] = decLen;
      recLen = (short int) (1 + fldLen[0] + fldLen[1]);
   }
   totSize = 1 + 32 + 32 * numCol + recLen * NxNy;

   /* Start writing the header. */
   /* Write the ID and date of last update. */
   fwrite (header, sizeof (char), 4, fp);
   /* Write number of records, size of header, then length of record. */
   FWRITE_LIT (&NxNy, sizeof (sInt4), 1, fp);
   si_temp = (short int) (32 + 32 * numCol + 1);
   FWRITE_LIT (&si_temp, sizeof (short int), 1, fp);
   FWRITE_LIT (&recLen, sizeof (short int), 1, fp);
   /* Write Reserved bytes.. 20 bytes of them.. all 0. */
   fwrite (reserved, sizeof (char), 20, fp);

   /* Write field descriptor array. */
   for (x = 0; x < numCol; x++) {
      /* Already taken care of padding with 0's since we did a strncpy */
      fwrite (names[x], sizeof (char), 11, fp);
      fputc (type[x], fp);
      /* Should be FWRITE_LIT, but doesn't matter. 11/8/2004 */
      fwrite (&(address), sizeof (sInt4), 1, fp); /* with address 0 */
      fputc (fldLen[x], fp);
      fputc (fldDec[x], fp);
      fwrite (reserved, sizeof (char), 14, fp); /* reserved already has 0's */
   }
   /* Write trailing header character. */
   uc_temp = 13;
   fputc (uc_temp, fp);

   /* Start writing the body. */
   if (f_verbose) {
      sprintf (formatBuf, " %%08ld%%0%dd%%0%dd%%%d.%df%%%d.%df%%%d.%df",
               fldLen[1], fldLen[2], fldLen[3], fldDec[3], fldLen[4],
               fldDec[4], fldLen[5], fldDec[5]);
   } else {
      sprintf (formatBuf, " %%08ld%%%d.%df", fldLen[1], fldDec[1]);
   }
   if ((!f_nMissing) || (attrib->f_miss == 0)) {
      for (y = 0; y < Ny; y++) {
         curData = grib_Data + y * Nx;
         /* Ny is 1 if we are creating BigPolyDbf files. */
         if (Ny == 1) {
            id = y + 1;
         } else {
            id = 10000 + y + 1;
         }
         /* May have 2 missing values (ie 9999 && 9998) in the .dbf file. */
         for (x = 0; x < Nx; x++) {
            if (f_verbose) {
               myCxy2ll (map, x + 1, y + 1, &lat, &lon);
               lat = myRound (lat, 5);
               lon = myRound (lon, 5);
               fprintf (fp, formatBuf, id, x + 1, y + 1, lon, lat,
                        ((floor (*curData * shift + .5)) / shift));
            } else {
               fprintf (fp, formatBuf, id, ((floor (*curData * shift + .5))
                                            / shift));
            }
            curData++;
            /* Ny is 1 if we are creating BigPolyDbf files. */
            if (Ny == 1) {
               id++;
            } else {
               id += 10000;
            }
         }
      }
   } else {
      numRec = 0;
      for (y = 0; y < Ny; y++) {
         curData = grib_Data + y * Nx;
         /* Ny is 1 if we are creating BigPolyDbf files. */
         if (Ny == 1) {
            id = y + 1;
         } else {
            id = 10000 + y + 1;
         }
         for (x = 0; x < Nx; x++) {
            if ((*curData != attrib->missPri) &&
                ((attrib->f_miss != 2) || (*curData != attrib->missSec))) {
               if (f_verbose) {
                  myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                  lat = myRound (lat, 5);
                  lon = myRound (lon, 5);
                  fprintf (fp, formatBuf, id, x + 1, y + 1, lon, lat,
                           ((floor (*curData * shift + .5)) / shift));
               } else {
                  fprintf (fp, formatBuf, id,
                           ((floor (*curData * shift + .5)) / shift));
               }
               numRec++;
            }
            curData++;
            /* Ny is 1 if we are creating BigPolyDbf files. */
            if (Ny == 1) {
               id++;
            } else {
               id += 10000;
            }
         }
      }
      /* Update file total # of records. */
      fseek (fp, 4, SEEK_SET);
      FWRITE_LIT (&numRec, sizeof (sInt4), 1, fp);
      totSize = 1 + 32 + 32 * numCol + recLen * numRec;
   }
   fclose (fp);

   /* Check that .dbf is now the correct file size. */
   return checkFileSize (filename, totSize);
}


/*****************************************************************************
 * CreateWxDbf() --
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This creates a .dbf file.  The .dbf file contains the actual info, which
 * is associated with the points.  The .dbf file format is very flexible but
 * unfortunately is mostly in ASCII.
 *   The reason this is separated from CreateDbf() is because the weather
 * data is a lot more complex than regular elements, so we need to provide
 * more columns to the user.
 *
 * ARGUMENTS
 *   filename = Name of file to save to. (Output)
 *     Nx, Ny = The dimensions of the grid. (Input)
 *  grib_Data = The actual data (assumed parsed by ParseGrid) (Input)
 *    element = Name of current variable. (Input)
 * f_nMissing = True if we should NOT store missing. (Input)
 *     attrib = What kind of missing value management is used. (Input)
 *     wxType = The lookup table associated with the weather. (Input)
 *  f_verbose = True if we want the verbose form. (Input)
 *        map = Holds the current map projection info to compute from. (In)
 *
 * FILES/DATABASES:
 *   Creates a .dbf file, which is a binary file (Little Endian) consisting of
 *      the info to associate with the lat/lon point shapes stored in the .shp
 *      and .shx files.
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = Problems opening the files.
 *
 * HISTORY
 *   5/2003 Arthur Taylor (MDL/RSIS): Created.
 *   5/2003 AAT: Decided to have 1,1 be lower left corner in .shp files.
 *   7/2003 AAT: Added "SimpleWx" column.
 *   1/2005 AAT: Modified for verbose output
 *
 * NOTES
 *****************************************************************************
 */
static int CreateWxDbf (char *filename, sInt4 Nx, sInt4 Ny,
                        double *grib_Data, const char *element,
                        sChar f_nMissing, gridAttribType *attrib,
                        sect2_WxType *WxType, char f_verbose, myMaparam *map)
{
   uChar header[] = { 3, 101, 4, 20 }; /* Header info for dbf. */
   int numColN = 25;    /* The number of columns we use. */
   int numColV = 29;    /* The number of columns we use. */
   char namesN[][11] = {
      "POINTID", "element", "WX-INDEX", "Visibility",
      "NDFDWxCode",
      "Weather_1", "WX-Inten_1", "Cover_1", "Hazard_1",
      "Weather_2", "WX-Inten_2", "Cover_2", "Hazard_2",
      "Weather_3", "WX-Inten_3", "Cover_3", "Hazard_3",
      "Weather_4", "WX-Inten_4", "Cover_4", "Hazard_4",
      "Weather_5", "WX-Inten_5", "Cover_5", "Hazard_5",
   };                   /* Field (column) names. */
   char namesV[][11] = {
      "POINTID", "X", "Y", "LON", "LAT",
      "element", "WX-INDEX", "Visibility", "NDFDWxCode",
      "Weather_1", "WX-Inten_1", "Cover_1", "Hazard_1",
      "Weather_2", "WX-Inten_2", "Cover_2", "Hazard_2",
      "Weather_3", "WX-Inten_3", "Cover_3", "Hazard_3",
      "Weather_4", "WX-Inten_4", "Cover_4", "Hazard_4",
      "Weather_5", "WX-Inten_5", "Cover_5", "Hazard_5",
   };                   /* Field (column) names. */
   char typeN[] = {
      'N', 'C', 'N', 'N',
      'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
   };                   /* Type of data for columns. */
   char typeV[] = {
      'N', 'N', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
      'C', 'N', 'N', 'N',
   };                   /* Type of data for columns. */
   uChar fldLenN[] = {
      8, 0, 8, 8,
      4,
      0, 3, 2, 10,
      0, 3, 2, 10,
      0, 3, 2, 10,
      0, 3, 2, 10,
      0, 3, 2, 10,
   };                   /* field len for columns. */
   uChar fldLenV[] = {
      8, 4, 4, 10, 9,
      0, 8, 8, 4,
      0, 3, 2, 10,
      0, 3, 2, 10,
      0, 3, 2, 10,
      0, 3, 2, 10,
      0, 3, 2, 10,
   };                   /* field len for columns. */
   uChar fldDecN[] = {
      0, 0, 0, 2,
      0,
      0, 0, 0, 0,
      0, 0, 0, 0,
      0, 0, 0, 0,
      0, 0, 0, 0,
      0, 0, 0, 0
   };                   /* field decimal lens for columns. */
   uChar fldDecV[] = {
      0, 0, 0, 5, 5,
      0, 0, 2, 0,
      0, 0, 0, 0,
      0, 0, 0, 0,
      0, 0, 0, 0,
      0, 0, 0, 0,
      0, 0, 0, 0
   };                   /* field decimal lens for columns. */
   int cnt;             /* Used to loop over fldLen[] to compute recLen */
   sInt4 reserved[] = { 0, 0, 0, 0, 0 }; /* need 20 bytes of 0. */
   sInt4 address = 0;   /* Address to find data? */
   sInt4 NxNy = Nx * Ny; /* The total number of values. */
   FILE *fp;            /* The open .dbf file pointer. */
   int x, y;            /* Current grid cell location. */
   short int recLen;    /* Size in bytes of one cell of data */
   sInt4 id;            /* A unique decimal id for the current cell. */
   double *curData;     /* A pointer to Grib data for current cell. */
   sInt4 totSize;       /* Total size of the .dbf file. */
   uChar uc_temp;       /* Temp. storage of type unsigned char. */
   short int si_temp;   /* Temp. storage of type short int. */
   char formBuf1[100];  /* Used to format float output sent to .dbf */
   char formBuf[NUM_UGLY_WORD][100]; /* Used to format weather columns */
   sInt4 numRec;        /* The total number of records actually stored. */
   char buffer[100];    /* Stores index if we can't look it up in table. */
   uInt4 index;         /* Current index into Wx table. */
   double vis;          /* Used to determine if we have "missing" vis. */
   char *ptr;           /* Help print english version of Wx elements. */
   uChar wx_inten;      /* Help print Wx code. */
   sInt4 HazCode;       /* Help print hazard codes. */
   uChar cover;         /* Help print coverage code. */
   int i;               /* Helps init ptr, wx_inten, hazard, and cover. */
   double lat;          /* Latitude of the current point for f_verbose */
   double lon;          /* Longitude of the current point for f_verbose */

   strncpy (filename + strlen (filename) - 3, "dbf", 3);
   if ((fp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }
/*
   codeLen = WxType->dataLen;
   CodeTable = (int *) malloc (codeLen * sizeof (int));
   for (i = 0; i < codeLen; i++) {
      CodeTable[i] = NDFD_WxTable (&(WxType->ugly[i]));
   }
*/

   /* Figure out the field lengths, and the fprintf format. */
   if (f_verbose) {
      fldLenV[5] = WxType->maxLen;
      sprintf (formBuf1, " %%08ld%%0%dd%%0%dd%%%d.%df%%%d.%df%%%ds%%08ld"
               "%%8.2f%%04ld", fldLenV[1], fldLenV[2], fldLenV[3],
               fldDecV[3], fldLenV[4], fldDecV[4], fldLenV[5]);
   } else {
      fldLenN[1] = WxType->maxLen;
      sprintf (formBuf1, " %%08ld%%%ds%%08ld%%8.2f%%04ld", fldLenN[1]);
   }

   /* Reserve enough space in each english column for "unknown" */
   for (i = 0; i < NUM_UGLY_WORD; i++) {
      if (WxType->maxEng[i] != 0) {
         sprintf (formBuf[i], "%%%ds%%03d%%02d%%010ld", WxType->maxEng[i]);
         if (f_verbose) {
            fldLenV[9 + i * 4] = WxType->maxEng[i];
         } else {
            fldLenN[5 + i * 4] = WxType->maxEng[i];
         }
      } else {
         sprintf (formBuf[i], "%%%ds%%01d%%01d%%01ld", 7);
         if (f_verbose) {
            fldLenV[9 + i * 4] = 7;
            fldLenV[10 + i * 4] = 1;
            fldLenV[11 + i * 4] = 1;
            fldLenV[12 + i * 4] = 1;
         } else {
            fldLenN[5 + i * 4] = 7;
            fldLenN[6 + i * 4] = 1;
            fldLenN[7 + i * 4] = 1;
            fldLenN[8 + i * 4] = 1;
         }
      }
   }

   if (f_verbose) {
      strncpy (namesV[5], element, 11);
      recLen = (short int) 1;
      for (cnt = 0; cnt < numColV; cnt++) {
         recLen += (short int) (fldLenV[cnt]);
      }
   } else {
      strncpy (namesN[1], element, 11);
      recLen = (short int) 1;
      for (cnt = 0; cnt < numColN; cnt++) {
         recLen += (short int) (fldLenN[cnt]);
      }
   }

   /* Start writing the header. */
   /* Write the ID and date of last update. */
   fwrite (header, sizeof (char), 4, fp);
   /* Write number of records, size of header, then length of record. */
   FWRITE_LIT (&NxNy, sizeof (sInt4), 1, fp);
   if (f_verbose) {
      si_temp = (short int) (32 + 32 * numColV + 1);
   } else {
      si_temp = (short int) (32 + 32 * numColN + 1);
   }
   FWRITE_LIT (&si_temp, sizeof (short int), 1, fp);
   FWRITE_LIT (&recLen, sizeof (short int), 1, fp);
   /* Write Reserved bytes.. 20 bytes of them.. all 0. */
   fwrite (reserved, sizeof (char), 20, fp);

   /* Write field descriptor array. */
   if (f_verbose) {
      for (x = 0; x < numColV; x++) {
         /* Already taken care of padding with 0's since we did a strncpy */
         fwrite (namesV[x], sizeof (char), 11, fp);
         fputc (typeV[x], fp);
         /* with address 0 */
         FWRITE_LIT (&(address), sizeof (sInt4), 1, fp);
         fputc (fldLenV[x], fp);
         fputc (fldDecV[x], fp);
         /* reserved already has 0's */
         fwrite (reserved, sizeof (char), 14, fp);
      }
   } else {
      for (x = 0; x < numColN; x++) {
         /* Already taken care of padding with 0's since we did a strncpy */
         fwrite (namesN[x], sizeof (char), 11, fp);
         fputc (typeN[x], fp);
         /* with address 0 */
         FWRITE_LIT (&(address), sizeof (sInt4), 1, fp);
         fputc (fldLenN[x], fp);
         fputc (fldDecN[x], fp);
         /* reserved already has 0's */
         fwrite (reserved, sizeof (char), 14, fp);
      }
   }
   /* Write trailing header character. */
   uc_temp = 13;
   fputc (uc_temp, fp);

   /* Start writing the body. */
   numRec = 0;
   for (y = 0; y < Ny; y++) {
      curData = grib_Data + y * Nx;
      /* Ny is 1 if we are creating BigPolyDbf files. */
      if (Ny == 1) {
         id = y + 1;
      } else {
         id = 10000 + y + 1;
      }
      for (x = 0; x < Nx; x++) {
         if ((!f_nMissing) || (attrib->f_miss == 0) ||
             ((attrib->f_miss == 1) && (*curData != attrib->missPri)) ||
             ((attrib->f_miss == 2) && (*curData != attrib->missPri)
              && (*curData != attrib->missSec))) {
            index = (uInt4) *curData;
            vis = 9999;
            if (index < WxType->dataLen) {
               if (WxType->f_valid[index]) {
                  if (WxType->ugly[index].minVis != 255) {
                     vis = WxType->ugly[index].minVis / 32.;
                  }
                  if (f_verbose) {
                     myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                     lat = myRound (lat, 5);
                     lon = myRound (lon, 5);
                     fprintf (fp, formBuf1, id, x + 1, y + 1, lon, lat,
                              WxType->data[index],
                              WxType->ugly[index].validIndex, vis,
                              (sInt4) WxType->ugly[index].SimpleCode);
                  } else {
                     fprintf (fp, formBuf1, id, WxType->data[index],
                              WxType->ugly[index].validIndex, vis,
                              (sInt4) WxType->ugly[index].SimpleCode);
                  }
                  for (i = 0; i < NUM_UGLY_WORD; i++) {
                     ptr = WxType->ugly[index].english[i];
                     if (ptr != NULL) {
                        wx_inten = WxType->ugly[index].wx_inten[i];
                        HazCode = WxType->ugly[index].HazCode[i];
                        cover = WxType->ugly[index].cover[i];
                        fprintf (fp, formBuf[i], ptr, wx_inten,
                                 cover, HazCode);
                     } else {
                        fprintf (fp, formBuf[i], "", 0, 0, 0);
                     }
                  }
               } else {
                  /* Handles a string of <Invalid><Invalid>... A rare event
                   * to begin with, typically this will be turned into a
                   * missing value in metaparse.c, but if we don't have any
                   * missing management, it sets f_valid to false. so we
                   * handle it here. */
                  if (f_verbose) {
                     myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                     lat = myRound (lat, 5);
                     lon = myRound (lon, 5);
                     fprintf (fp, formBuf1, id, x + 1, y + 1, lon, lat,
                              WxType->data[index],
                              WxType->ugly[index].validIndex, vis,
                              (sInt4) 9999);
                  } else {
                     fprintf (fp, formBuf1, id, WxType->data[index],
                              WxType->ugly[index].validIndex, vis,
                              (sInt4) 9999);
                  }
                  fprintf (fp, formBuf[0], "Unkown", 0, 0, 0);
                  fprintf (fp, formBuf[1], "Unkown", 0, 0, 0);
                  fprintf (fp, formBuf[2], "Unkown", 0, 0, 0);
                  fprintf (fp, formBuf[3], "Unkown", 0, 0, 0);
                  fprintf (fp, formBuf[4], "Unkown", 0, 0, 0);
               }
            } else {
               sprintf (buffer, "%ld", (long int) index);
/*               fprintf (fp, formBuf1, id, buffer, index, vis, index); */
               if (f_verbose) {
                  myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                  lat = myRound (lat, 5);
                  lon = myRound (lon, 5);
                  fprintf (fp, formBuf1, id, x + 1, y + 1, lon, lat, buffer,
                           0, vis, 0);
               } else {
                  fprintf (fp, formBuf1, id, buffer, 0, vis, 0);
               }
               fprintf (fp, formBuf[0], "Unkown", 0, 0, 0);
               fprintf (fp, formBuf[1], "Unkown", 0, 0, 0);
               fprintf (fp, formBuf[2], "Unkown", 0, 0, 0);
               fprintf (fp, formBuf[3], "Unkown", 0, 0, 0);
               fprintf (fp, formBuf[4], "Unkown", 0, 0, 0);
            }
            numRec++;
         }
         curData++;
         /* Ny is 1 if we are creating BigPolyDbf files. */
         if (Ny == 1) {
            id++;
         } else {
            id += 10000;
         }
      }
   }
   /* Update file total # of records. */
   fseek (fp, 4, SEEK_SET);
   FWRITE_LIT (&numRec, sizeof (sInt4), 1, fp);
   if (f_verbose) {
      totSize = 1 + 32 + 32 * numColV + recLen * numRec;
   } else {
      totSize = 1 + 32 + 32 * numColN + recLen * numRec;
   }
   fclose (fp);

   /* Check that .dbf is now the correct file size. */
   return checkFileSize (filename, totSize);
}

static int CreateWWADbf (char *filename, sInt4 Nx, sInt4 Ny,
                        double *grib_Data, const char *element,
                        sChar f_nMissing, gridAttribType *attrib,
                        sect2_HazardType *HazType, char f_verbose, myMaparam *map)
{
   uChar header[] = { 3, 101, 4, 20 }; /* Header info for dbf. */
   int numColN = 9;    /* The number of columns we use. */
   int numColV = 13;    /* The number of columns we use. */
   char namesN[][11] = {
      "POINTID", "element", "WWA-INDEX",
      "WWA_Code",
      "WWA_1",
      "WWA_2",
      "WWA_3",
      "WWA_4",
      "WWA_5",
   };                   /* Field (column) names. */
   char namesV[][11] = {
      "POINTID", "X", "Y", "LON", "LAT",
      "element", "WWA-INDEX",
      "WWA_Code",
      "WWA_1",
      "WWA_2",
      "WWA_3",
      "WWA_4",
      "WWA_5",
   };                   /* Field (column) names. */
   char typeN[] = {
      'N', 'C', 'N',
      'N',
      'C',
      'C',
      'C',
      'C',
      'C',
   };                   /* Type of data for columns. */
   char typeV[] = {
      'N', 'N', 'N', 'N', 'N',
      'C', 'N', 'N',
      'C',
      'C',
      'C',
      'C',
      'C',
   };                   /* Type of data for columns. */
   uChar fldLenN[] = {
      8, 0, 8,
      4,
      0,
      0,
      0,
      0,
      0,
   };                   /* field len for columns. */
   uChar fldLenV[] = {
      8, 4, 4, 10, 9,
      0, 8, 4,
      0,
      0,
      0,
      0,
      0,
   };                   /* field len for columns. */
   uChar fldDecN[] = {
      0, 0, 0,
      0,
      0,
      0,
      0,
      0,
      0,
   };                   /* field decimal lens for columns. */
   uChar fldDecV[] = {
      0, 0, 0, 5, 5,
      0, 0, 0,
      0,
      0,
      0,
      0,
      0,
   };                   /* field decimal lens for columns. */
   int cnt;             /* Used to loop over fldLen[] to compute recLen */
   sInt4 reserved[] = { 0, 0, 0, 0, 0 }; /* need 20 bytes of 0. */
   sInt4 address = 0;   /* Address to find data? */
   sInt4 NxNy = Nx * Ny; /* The total number of values. */
   FILE *fp;            /* The open .dbf file pointer. */
   int x, y;            /* Current grid cell location. */
   short int recLen;    /* Size in bytes of one cell of data */
   sInt4 id;            /* A unique decimal id for the current cell. */
   double *curData;     /* A pointer to Grib data for current cell. */
   sInt4 totSize;       /* Total size of the .dbf file. */
   uChar uc_temp;       /* Temp. storage of type unsigned char. */
   short int si_temp;   /* Temp. storage of type short int. */
   char formBuf1[100];  /* Used to format float output sent to .dbf */
   char formBuf[NUM_HAZARD_WORD][100]; /* Used to format hazard columns */
   sInt4 numRec;        /* The total number of records actually stored. */
   char buffer[100];    /* Stores index if we can't look it up in table. */
   uInt4 index;         /* Current index into WWA table. */
   char *ptr;           /* Help print english version of WWA elements. */
   int i;               /* Helps init ptr, wwa_inten. */
   double lat;          /* Latitude of the current point for f_verbose */
   double lon;          /* Longitude of the current point for f_verbose */

   strncpy (filename + strlen (filename) - 3, "dbf", 3);
   if ((fp = fopen (filename, "wb")) == NULL) {
      errSprintf ("ERROR: Problems opening %s for write.", filename);
      return -1;
   }

   /* Figure out the field lengths, and the fprintf format. */
   if (f_verbose) {
      fldLenV[5] = HazType->maxLen;
      sprintf (formBuf1, " %%08ld%%0%dd%%0%dd%%%d.%df%%%d.%df%%%ds%%08ld"
               "%%04ld", fldLenV[1], fldLenV[2], fldLenV[3],
               fldDecV[3], fldLenV[4], fldDecV[4], fldLenV[5]);
   } else {
      fldLenN[1] = HazType->maxLen;
      sprintf (formBuf1, " %%08ld%%%ds%%08ld%%04ld", fldLenN[1]);
   }

   /* Reserve enough space in each english column for "unknown" */
   for (i = 0; i < NUM_HAZARD_WORD; i++) {
      if (HazType->maxEng[i] != 0) {
         sprintf (formBuf[i], "%%%ds", HazType->maxEng[i]);
         if (f_verbose) {
            fldLenV[8 + i] = HazType->maxEng[i];
         } else {
            fldLenN[4 + i] = HazType->maxEng[i];
         }
      } else {
         sprintf (formBuf[i], "%%%ds", 7);
         if (f_verbose) {
            fldLenV[8 + i] = 7;
         } else {
            fldLenN[4 + i] = 7;
         }
      }
   }

   if (f_verbose) {
      strncpy (namesV[5], element, 11);
      recLen = (short int) 1;
      for (cnt = 0; cnt < numColV; cnt++) {
         recLen += (short int) (fldLenV[cnt]);
      }
   } else {
      strncpy (namesN[1], element, 11);
      recLen = (short int) 1;
      for (cnt = 0; cnt < numColN; cnt++) {
         recLen += (short int) (fldLenN[cnt]);
      }
   }

   /* Start writing the header. */
   /* Write the ID and date of last update. */
   fwrite (header, sizeof (char), 4, fp);
   /* Write number of records, size of header, then length of record. */
   FWRITE_LIT (&NxNy, sizeof (sInt4), 1, fp);
   if (f_verbose) {
      si_temp = (short int) (32 + 32 * numColV + 1);
   } else {
      si_temp = (short int) (32 + 32 * numColN + 1);
   }
   FWRITE_LIT (&si_temp, sizeof (short int), 1, fp);
   FWRITE_LIT (&recLen, sizeof (short int), 1, fp);
   /* Write Reserved bytes.. 20 bytes of them.. all 0. */
   fwrite (reserved, sizeof (char), 20, fp);

   /* Write field descriptor array. */
   if (f_verbose) {
      for (x = 0; x < numColV; x++) {
         /* Already taken care of padding with 0's since we did a strncpy */
         fwrite (namesV[x], sizeof (char), 11, fp);
         fputc (typeV[x], fp);
         /* with address 0 */
         FWRITE_LIT (&(address), sizeof (sInt4), 1, fp);
         fputc (fldLenV[x], fp);
         fputc (fldDecV[x], fp);
         /* reserved already has 0's */
         fwrite (reserved, sizeof (char), 14, fp);
      }
   } else {
      for (x = 0; x < numColN; x++) {
         /* Already taken care of padding with 0's since we did a strncpy */
         fwrite (namesN[x], sizeof (char), 11, fp);
         fputc (typeN[x], fp);
         /* with address 0 */
         FWRITE_LIT (&(address), sizeof (sInt4), 1, fp);
         fputc (fldLenN[x], fp);
         fputc (fldDecN[x], fp);
         /* reserved already has 0's */
         fwrite (reserved, sizeof (char), 14, fp);
      }
   }
   /* Write trailing header character. */
   uc_temp = 13;
   fputc (uc_temp, fp);

   /* Start writing the body. */
   numRec = 0;
   for (y = 0; y < Ny; y++) {
      curData = grib_Data + y * Nx;
      /* Ny is 1 if we are creating BigPolyDbf files. */
      if (Ny == 1) {
         id = y + 1;
      } else {
         id = 10000 + y + 1;
      }
      for (x = 0; x < Nx; x++) {
         if ((!f_nMissing) || (attrib->f_miss == 0) ||
             ((attrib->f_miss == 1) && (*curData != attrib->missPri)) ||
             ((attrib->f_miss == 2) && (*curData != attrib->missPri)
              && (*curData != attrib->missSec))) {
            index = (uInt4) *curData;
            if (index < HazType->dataLen) {
               if (HazType->f_valid[index]) {
                  if (f_verbose) {
                     myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                     lat = myRound (lat, 5);
                     lon = myRound (lon, 5);
                     fprintf (fp, formBuf1, id, x + 1, y + 1, lon, lat,
                              HazType->data[index],
                              HazType->haz[index].validIndex,
                              (sInt4) HazType->haz[index].SimpleCode);
                  } else {
                     fprintf (fp, formBuf1, id, HazType->data[index],
                              HazType->haz[index].validIndex,
                              (sInt4) HazType->haz[index].SimpleCode);
                  }
                  for (i = 0; i < NUM_HAZARD_WORD; i++) {
                     ptr = HazType->haz[index].english[i];
                     if (ptr != NULL) {
                        fprintf (fp, formBuf[i], ptr);
                     } else {
                        fprintf (fp, formBuf[i], "");
                     }
                  }
               } else {
                  /* Handles a string of <Invalid><Invalid>... A rare event
                   * to begin with, typically this will be turned into a
                   * missing value in metaparse.c, but if we don't have any
                   * missing management, it sets f_valid to false. so we
                   * handle it here. */
                  if (f_verbose) {
                     myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                     lat = myRound (lat, 5);
                     lon = myRound (lon, 5);
                     fprintf (fp, formBuf1, id, x + 1, y + 1, lon, lat,
                              HazType->data[index],
                              HazType->haz[index].validIndex,
                              (sInt4) 9999);
                  } else {
                     fprintf (fp, formBuf1, id, HazType->data[index],
                              HazType->haz[index].validIndex,
                              (sInt4) 9999);
                  }
                  fprintf (fp, formBuf[0], "Unkown");
                  fprintf (fp, formBuf[1], "Unkown");
                  fprintf (fp, formBuf[2], "Unkown");
                  fprintf (fp, formBuf[3], "Unkown");
                  fprintf (fp, formBuf[4], "Unkown");
               }
            } else {
               sprintf (buffer, "%ld", (long int) index);
/*               fprintf (fp, formBuf1, id, buffer, index, vis, index); */
               if (f_verbose) {
                  myCxy2ll (map, x + 1, y + 1, &lat, &lon);
                  lat = myRound (lat, 5);
                  lon = myRound (lon, 5);
                  fprintf (fp, formBuf1, id, x + 1, y + 1, lon, lat, buffer,
                           0, 0);
               } else {
                  fprintf (fp, formBuf1, id, buffer, 0, 0);
               }
               fprintf (fp, formBuf[0], "Unkown");
               fprintf (fp, formBuf[1], "Unkown");
               fprintf (fp, formBuf[2], "Unkown");
               fprintf (fp, formBuf[3], "Unkown");
               fprintf (fp, formBuf[4], "Unkown");
            }
            numRec++;
         }
         curData++;
         /* Ny is 1 if we are creating BigPolyDbf files. */
         if (Ny == 1) {
            id++;
         } else {
            id += 10000;
         }
      }
   }
   /* Update file total # of records. */
   fseek (fp, 4, SEEK_SET);
   FWRITE_LIT (&numRec, sizeof (sInt4), 1, fp);
   if (f_verbose) {
      totSize = 1 + 32 + 32 * numColV + recLen * numRec;
   } else {
      totSize = 1 + 32 + 32 * numColN + recLen * numRec;
   }
   fclose (fp);

   /* Check that .dbf is now the correct file size. */
   return checkFileSize (filename, totSize);
}

/*****************************************************************************
 * GridCompute() -- Review 12/2002
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This computes all the lat/lon values of the grid, and returns the
 * bounding lat/lon region of the grid.
 *   If f_corner, then it computes the 4 corners of the grid cells, instead
 * of the center points.
 *
 * ARGUMENTS
 *      map = Holds the current map projection info to compute from. (Input)
 *      gds = Grid Definiton from the parsed GRIB msg to write. (Input)
 *       dp = A lat/lon array large enough for gds->Nx*Ny points. (Output)
 * f_corner = 1 compute corners of grid cells, 0 compute centers. (Input)
 *
 * FILES/DATABASES: None
 *
 * RETURNS: void
 *
 * HISTORY
 *  11/2002 Arthur Taylor (MDL/RSIS): Created.
 *  12/2002 (RY,FC,MA,&TB): Code Review.
 *   4/2003 AAT: Added f_corner option.
 *   5/2003 AAT: Removed requirement for computation of min/max lat/lon.
 *   5/2003 AAT: Decided to have 1,1 be lower left corner in .shp files.
 *   3/2004 AAT: Switched to 6 decimal accuracy on lat/lon
 *               (because that is native of GRIB2 messages)
 *
 * NOTES
 * 1) Round to 6 decimals (because that is the accuracy of GRIB2 messages)
 *    (6 is too much for HP/Linux match)
 *    (7 is too much for Linux/Cygwin match)
 * 2) The grid for spatial analyst starts in upper left corner.
 *    So for consistancy sake we also start in upper left corner.
 *****************************************************************************
 */
static void GridCompute (myMaparam *map, gdsType *gds, LatLon *dp,
                         sChar f_corner, sChar LatLon_Decimal)
{
   LatLon *cur_dp;      /* Current lat/lon point. */
   uInt4 x, y;          /* loop index while computing the grid. */

   cur_dp = dp;
   for (y = 1; y <= gds->Ny + f_corner; y++) {
      for (x = 1; x <= gds->Nx + f_corner; x++) {
         myCxy2ll (map, x - f_corner * .5, y - f_corner * .5, &cur_dp->lat,
                   &cur_dp->lon);
         cur_dp->lat = myRound (cur_dp->lat, LatLon_Decimal);
         cur_dp->lon = myRound (cur_dp->lon, LatLon_Decimal);
         cur_dp++;
      }
   }
}

/*****************************************************************************
 * gribWriteShp() -- Review 12/2002
 *
 * Arthur Taylor / MDL
 *
 * PURPOSE
 *   This controls the creation of a .shp/.shx/.dbf file set.  These are
 * Esri ArcView specific file formats, but it was felt that we should provide
 * the capability of writing to them as it is a common GIS package.  After
 * it generates the .shp file set, it then creates a .ave script which can
 * be run to make sure that the projection is set in ArcView the way the Grib
 * data wanted it.  Main reason is to get the correct radius of the earth.
 *
 * ARGUMENTS
 *   Filename = Name of file to save to. (Input)
 *  grib_Data = The grib2 data to write. (Input)
 *        gds = Grid Definition from the parsed GRIB msg to write. (Input)
 *     f_poly = 2 => BigPoly, 1 => polygon, 0 => point .shp file, (Input)
 * f_nMissing = True if we do not want missing values in the .shp file. (In)
 *    decimal = How many decimals to round to. (Input)
 *  f_verbose = True if we want the verbose output. (Input)
 *
 * FILES/DATABASES:
 *   Either Calls CreateShpPnt to create the Esri .shp and .shx files.
 *   or    Calls CreateShpPoly to create the ESRI .shp and .shx files.
 *   Calls CreateDbf to create the Esri .dbf file.
 *   Calls gribWriteEsriAve to create the Esri .ave script.
 *
 * RETURNS: int (could use errSprintf())
 *  0 = OK
 * -1 = illegal declaration of the grid size.
 * -2 = Problems opening the files.
 * -3 = un-supported map projection.
 * -4 = Error in CreateShpPnt
 * -5 = Error in CreateDbf
 * -6 = Error in CreatePrj
 *  1 = invalid parameters for gribWriteEsriAve.
 *
 * HISTORY
 *   9/2002 Arthur Taylor (MDL/RSIS): Created.
 *  12/2002 (RY,FC,MA,&TB): Code Review.
 *   4/2003 AAT: Added calls to CreateShpPoly.
 *   5/2003 AAT: Switched from gdsType to grib_MetaData type so we can use
 *          the info record to determine good field lengths for the .dbf file.
 *   5/2003 AAT: Added rounding to decimal.
 *   5/2003 AAT: Enabled other spherical earths.
 *   6/2003 AAT: Switched to a Warning and then averaging if Dx != Dy.
 *   7/2003 AAT: Bug dealing with min = max = 0 for dbf file.
 *   7/2003 AAT: Proper handling of Dx != Dy.
 *   7/2003 AAT: switched to checking against element name for Wx instead
 *          of pds2.sect2.ptrType == GS2_WXTYPE
 *  10/2003 AAT: Added Calls to CreateBigPolyShp()
 *   1/2005 AAT: Added Call to CreatePrj()
 *   1/2005 AAT: Modified for verbose output
 *
 * NOTES
 * 1) Order is .shp/.shx then .dbf, then .ave.  If .ave doesn't work they have
 * the important stuff.
 *****************************************************************************
 */
int gribWriteShp (const char *Filename, double *grib_Data,
                  grib_MetaData *meta, sChar f_poly, sChar f_nMissing,
                  sChar decimal, sChar LatLon_Decimal, char f_verbose)
{
   myMaparam map;       /* Used to compute the grid lat/lon points. */
   char *filename;      /* local copy of the filename. */
   int nameLen;         /* length of filename. */
   LatLon *dp;          /* Array of lat/lon points. */
   double orient;       /* Orientation longitude of projection (where N is
                         * up.) (between -180 and 180) */
   gdsType *gds = &(meta->gds); /* Simplifies references to the gds data. */
   double max, min;     /* Used to figure out the largest number we have to
                         * store in the .dbf file. */
   uChar fieldLen1;     /* How many char to store the max number in. */
   uChar fieldLen2;     /* How many char to store the min number in. */
   polyType *poly;      /* list of chains that represent large polygons. */
   int numPoly;         /* number of element in poly. */
   double *polyData;    /* Data values for each poly (no missing values) */
   char EsriName[11];   /* element name shortened to 11 characters */
   int len;             /* Length of desired element name. */

   nameLen = strlen (Filename);
   if (nameLen < 4) {
      errSprintf ("ERROR: File %s is too short in length (it may need an "
                  "extension?)", Filename);
      return -2;
   }
   filename = (char *) malloc (nameLen + 1 * sizeof (char));
   strncpy (filename, Filename, nameLen);
   filename[nameLen] = '\0';

   /* Check that gds is valid before setting up map projection. */
   if (GDSValid (&meta->gds) != 0) {
      preErrSprintf ("ERROR: Sect3 was not Valid.\n");
      free (filename);
      return -3;
   }
   /* Set up the map projection. */
   SetMapParamGDS (&map, gds);
   /* Create the .shp/.shx files */
   if (f_poly == 1) {   /* Small poly */
      dp = (LatLon *) malloc ((gds->Nx + 1) * (gds->Ny + 1) *
                              sizeof (LatLon));
      /* Compute the lat/lon grid. */
      GridCompute (&map, gds, dp, 1, LatLon_Decimal);
      if (CreateShpPoly (filename, dp, gds->Nx, gds->Ny, f_nMissing,
                         grib_Data, &(meta->gridAttrib), LatLon_Decimal)
          != 0) {
         free (dp);
         free (filename);
         return -4;
      }
      free (dp);
   } else if (f_poly == 2) { /* Big poly */
      NewPolys (&poly, &numPoly);

      /* Convert the grid to a list of polygon chains. */
      /* Assumes data came from ParseGrid() so scan flag == "0100". */
      Grid2BigPoly (&poly, &numPoly, gds->Nx, gds->Ny, grib_Data);

      /* Following removes all "missing" values from the chains. */
      gribCompactPolys (poly, &numPoly, f_nMissing, &(meta->gridAttrib),
                        &polyData);

      /* Following is for experimenting with dateline issue. */
      ConvertChain2LtLn (poly, numPoly, &map, LatLon_Decimal);

      /* Following saves the chain of lat/lons to a shp/shx file. */
      if (CreateBigPolyShp (filename, poly, numPoly) != 0) {
         FreePolys (poly, numPoly);
         free (filename);
         return -4;
      }

      FreePolys (poly, numPoly);
   } else {             /* point */
#ifdef TEST
      dp = (LatLon *) malloc (gds->numPts * (sizeof (LatLon)));
      /* Compute the lat/lon grid. */
      GridCompute (&map, gds, dp, 0, LatLon_Decimal);
#endif
      if (CreateShpPnt (filename, &map, gds, LatLon_Decimal, f_nMissing,
                        grib_Data, &(meta->gridAttrib)) != 0) {
#ifdef TEST
         free (dp);
#endif
         free (filename);
         return -4;
      }
#ifdef TEST
      free (dp);
#endif
   }

   /* Create the .prj file */
   if (CreatePrj (filename, gds) != 0) {
      free (filename);
      return -6;
   }

   /* Since ESRI was lazy, we null terminate the .dbf column name.  Excel
    * doesn't need this, but ESRI does. */
   strncpy (EsriName, meta->element, 10);
   len = strlen (meta->element);
   if (len > 10) {
      if (strncmp (meta->element, "Prob", 4) == 0) {
         EsriName[0] = 'P';
         if ((len > 13) && (strncmp (meta->element + 4, "Wind", 4) == 0)) {
            EsriName[1] = 'w';
            strncpy (EsriName + 2, meta->element + 6, 10);
         } else {
            strncpy (EsriName + 1, meta->element + 4, 10);
         }
      }
   }
   EsriName[10] = '\0';


   /* Create the .dbf files */
   if (strcmp (EsriName, "Wx") == 0) {
      if (f_poly == 2) {
         if (CreateWxDbf (filename, numPoly, 1, polyData, EsriName,
                          f_nMissing, &(meta->gridAttrib),
                          (sect2_WxType *) &(meta->pds2.sect2.wx), 0,
                          &map) != 0) {
            free (filename);
            free (polyData);
            return -5;
         }
         free (polyData);
      } else {
         if (CreateWxDbf (filename, gds->Nx, gds->Ny, grib_Data,
                          EsriName, f_nMissing, &(meta->gridAttrib),
                          (sect2_WxType *) &(meta->pds2.sect2.wx),
                          f_verbose, &map) != 0) {
            free (filename);
            return -5;
         }
      }
   } else if (strcmp (EsriName, "WWA") == 0) {
      if (f_poly == 2) {
         if (CreateWWADbf (filename, numPoly, 1, polyData, EsriName,
                          f_nMissing, &(meta->gridAttrib),
                          (sect2_HazardType *) &(meta->pds2.sect2.hazard), 0,
                          &map) != 0) {
            free (filename);
            free (polyData);
            return -5;
         }
         free (polyData);
      } else {
         if (CreateWWADbf (filename, gds->Nx, gds->Ny, grib_Data,
                          EsriName, f_nMissing, &(meta->gridAttrib),
                          (sect2_HazardType *) &(meta->pds2.sect2.hazard),
                          f_verbose, &map) != 0) {
            free (filename);
            return -5;
         }
      }
   } else {
      max = meta->gridAttrib.max;
      min = meta->gridAttrib.min;
      switch (meta->gridAttrib.f_miss) {
         case 2:
            if (max < meta->gridAttrib.missSec) {
               max = meta->gridAttrib.missSec;
            } else if (min > meta->gridAttrib.missSec) {
               min = meta->gridAttrib.missSec;
            }
            /* Intentionally fall through. */
         case 1:
            if (max < meta->gridAttrib.missPri) {
               max = meta->gridAttrib.missPri;
            } else if (min > meta->gridAttrib.missPri) {
               min = meta->gridAttrib.missPri;
            }
      }
      if (max >= 1) {
         fieldLen1 = decimal + 1 + (int) (1 + log10 (max));
      } else if (max >= 0) {
         fieldLen1 = decimal + 1 + 1;
      } else if (max > -1) {
         fieldLen1 = decimal + 1 + 1 + 1;
      } else {
         fieldLen1 = decimal + 1 + (int) (1 + log10 (-1 * max)) + 1;
      }
      if (min >= 1) {
         fieldLen2 = decimal + 1 + (int) (1 + log10 (min));
      } else if (min >= 0) {
         fieldLen2 = decimal + 1 + 1;
      } else if (min > -1) {
         fieldLen2 = decimal + 1 + 1 + 1;
      } else {
         fieldLen2 = decimal + 1 + (int) (1 + log10 (-1 * min)) + 1;
      }
      if (fieldLen1 > fieldLen2) {
         if (f_poly == 2) {
            if (CreateDbf (filename, numPoly, 1, polyData, EsriName,
                           fieldLen1, decimal, f_nMissing,
                           &(meta->gridAttrib), 0, &map) != 0) {
               free (filename);
               free (polyData);
               return -5;
            }
            free (polyData);
         } else {
            if (CreateDbf (filename, gds->Nx, gds->Ny, grib_Data,
                           EsriName, fieldLen1, decimal, f_nMissing,
                           &(meta->gridAttrib), f_verbose, &map) != 0) {
               free (filename);
               return -5;
            }
         }
      } else {
         if (f_poly == 2) {
            if (CreateDbf (filename, numPoly, 1, polyData, EsriName,
                           fieldLen2, decimal, f_nMissing,
                           &(meta->gridAttrib), 0, &map) != 0) {
               free (filename);
               free (polyData);
               return -5;
            }
            free (polyData);
         } else {
            if (CreateDbf (filename, gds->Nx, gds->Ny, grib_Data,
                           EsriName, fieldLen2, decimal, f_nMissing,
                           &(meta->gridAttrib), f_verbose, &map) != 0) {
               free (filename);
               return -5;
            }
         }
      }
   }

   /* Create the .ave file */
   strncpy (filename + nameLen - 3, "ave", 3);
   orient = gds->orientLon;
   while (orient > 180) {
      orient -= 360;
   }
   while (orient < -180) {
      orient += 360;
   }
   if (gribWriteEsriAve (filename, gds, orient) != 0) {
      preErrSprintf ("gribWriteEsriAve had the following error\n");
      free (filename);
      return 1;
   }
   free (filename);
   return 0;
}