g2_addfield.c

Go to the documentation of this file.

 
#include <stdio.h>
#include <stdlib.h>
#include "grib2.h"
 
g2int getdim(unsigned char *, g2int *, g2int *, g2int *);
g2int getpoly(unsigned char *, g2int *, g2int *, g2int *);
void simpack(g2float *,  g2int,  g2int *,  unsigned char *,  g2int *);
void cmplxpack(g2float *,  g2int,  g2int,  g2int *,  unsigned char *,  g2int *);
void specpack(g2float *, g2int, g2int, g2int, g2int, g2int *, unsigned char *, 
              g2int *);
#ifdef USE_PNG
void pngpack(g2float *, g2int, g2int, g2int *, unsigned char *, g2int *);
#endif  /* USE_PNG */
#if defined USE_JPEG2000 || defined USE_OPENJPEG
void jpcpack(g2float *, g2int, g2int, g2int *, unsigned char *, g2int *);
#endif  /* USE_JPEG2000 */
 
 
g2int
g2_addfield(unsigned char *cgrib, g2int ipdsnum, g2int *ipdstmpl,
            g2float *coordlist, g2int numcoord, g2int idrsnum, g2int *idrstmpl,
            g2float *fld, g2int ngrdpts, g2int ibmap, g2int *bmap)
{
    g2int ierr = 0;
    static unsigned char G = 0x47; /* 'G' */
    static unsigned char R = 0x52; /* 'R' */
    static unsigned char I = 0x49; /* 'I' */
    static unsigned char B = 0x42; /* 'B' */
    static unsigned char s7 = 0x37; /* '7' */
 
    unsigned char *cpack;
    static g2int zero = 0, one = 1, four = 4, five = 5, six = 6, seven = 7;
    const g2int minsize = 50000;
    g2int iofst, ibeg, lencurr, len, nsize;
    g2int ilen, isecnum, i, nbits, temp, left;
    g2int ibmprev, j, lcpack, ioctet, newlen, ndpts;
    g2int lensec4, lensec5, lensec6, lensec7;
    g2int issec3 = 0, isprevbmap = 0, lpos3 = 0, JJ, KK, MM;
    g2int *coordieee;
    g2int width, height, iscan, itemp;
    g2float *pfld;
    gtemplate *mappds, *mapdrs;
    unsigned int allones = 4294967295u;
 
    /* Check to see if beginning of GRIB message exists. */
    if (cgrib[0] != G || cgrib[1] != R || cgrib[2] != I || cgrib[3] != B)
    {
        printf("g2_addfield: GRIB not found in given message.\n");
        printf("g2_addfield: Call to routine g2_create required to initialize GRIB messge.\n");
        ierr = -1;
        return ierr;
    }
 
    /* Get current length of GRIB message. */
    gbit(cgrib, &lencurr, 96, 32);
 
    /* Check to see if GRIB message is already complete. */
    if (cgrib[lencurr-4] == s7 && cgrib[lencurr-3] == s7 &&
         cgrib[lencurr-2] == s7 && cgrib[lencurr-1] == s7)
    {
        printf("g2_addfield: GRIB message already complete.  Cannot add new section.\n");
        ierr = -2;
        return ierr;
    }
 
    /* Loop through all current sections of the GRIB message to find
     * the last section number. */
    len = 16;    /* length of Section 0 */
    for (;;)
    {
        /* Get number and length of next section. */
        iofst = len * 8;
        gbit(cgrib, &ilen, iofst, 32);
        iofst = iofst + 32;
        gbit(cgrib, &isecnum, iofst, 8);
        iofst = iofst + 8;
 
        /* Check if previous Section 3 exists. */
        if (isecnum == 3)
        {
            issec3 = 1;
            lpos3 = len;
        }
        /* Check if a previous defined bitmap exists. */
        if (isecnum == 6)
        {
            gbit(cgrib, &ibmprev, iofst, 8);
            iofst = iofst + 8;
            if (ibmprev >= 0 && ibmprev <= 253)
                isprevbmap = 1;
        }
        len = len + ilen;
 
        /* Exit loop if last section reached. */
        if (len == lencurr)
            break;
 
        /* If byte count for each section doesn't match current */
        /* total length, then there is a problem. */
        if (len > lencurr)
        {
            printf("g2_addfield: Section byte counts don''t add to total.\n");
            printf("g2_addfield: Sum of section byte counts = %ld\n", len);
            printf("g2_addfield: Total byte count in Section 0 = %ld\n", lencurr);
            ierr = -3;
            return ierr;
        }
    }
 
    /* Sections 4 through 7 can only be added after section 3 or 7. */
    if ((isecnum != 3) && (isecnum != 7))
    {
        printf("g2_addfield: Sections 4-7 can only be added after Section 3 or 7.\n");
        printf("g2_addfield: Section ',isecnum,' was the last found in given GRIB message.\n");
        ierr = -4;
        return ierr;
    }
    else if (!issec3)
    {
        /* Sections 4 through 7 can only be added if section 3 was previously defined. */
        printf("g2_addfield: Sections 4-7 can only be added if Section 3 was previously included.\n");
        printf("g2_addfield: Section 3 was not found in given GRIB message.\n");
        printf("g2_addfield: Call to routine addgrid required to specify Grid definition.\n");
        ierr = -6;
        return ierr;
    }
 
    /* Add Section 4  - Product Definition Section. */
    ibeg = lencurr * 8;        /* Calculate offset for beginning of section 4 */
    iofst = ibeg + 32;         /* leave space for length of section */
    sbit(cgrib, &four, iofst, 8);     /* Store section number (4) */
    iofst = iofst + 8;
    sbit(cgrib, &numcoord, iofst, 16);   /* Store num of coordinate values */
    iofst = iofst + 16;
    sbit(cgrib, &ipdsnum, iofst, 16);    /* Store Prod Def Template num. */
    iofst = iofst + 16;
 
    /* Get Product Definition Template. */
    if (!(mappds = getpdstemplate(ipdsnum)))
    {          /* undefined template */
        ierr = -5;
        return ierr;
    }
 
    /* Extend the Product Definition Template, if necessary.  The */
    /* number of values in a specific template may vary depending on */
    /* data specified in the "static" part of the template. */
    if (mappds->needext)
    {
        free(mappds);
        mappds = extpdstemplate(ipdsnum, ipdstmpl);
    }
 
    /* Pack up each input value in array ipdstmpl into the the */
    /* appropriate number of octets, which are specified in */
    /* corresponding entries in array mappds. */
    for (i = 0; i < mappds->maplen; i++)
    {
        nbits = abs(mappds->map[i]) * 8;
        if ((mappds->map[i] >= 0) || (ipdstmpl[i] >= 0))
            sbit(cgrib, ipdstmpl + i, iofst, nbits);
        else
        {
            sbit(cgrib, &one, iofst, 1);
            temp = abs(ipdstmpl[i]);
            sbit(cgrib, &temp, iofst + 1, nbits-1);
        }
        iofst = iofst + nbits;
    }
 
    /* Pack template extension, if appropriate. */
    j = mappds->maplen;
    if (mappds->needext && (mappds->extlen > 0))
    {
        for (i = 0; i < mappds->extlen; i++)
        {
            nbits = abs(mappds->ext[i]) * 8;
            if ((mappds->ext[i] >= 0) || (ipdstmpl[j] >= 0))
                sbit(cgrib, ipdstmpl + j, iofst, nbits);
            else
            {
                sbit(cgrib, &one, iofst, 1);
                temp = abs(ipdstmpl[j]);
                sbit(cgrib, &temp, iofst + 1, nbits-1);
            }
            iofst = iofst + nbits;
            j++;
        }
    }
    free(mappds);
 
    /* Add Optional list of vertical coordinate values after the */
    /* Product Definition Template, if necessary. */
    if (numcoord != 0)
    {
        coordieee = calloc(numcoord, sizeof(g2int));
        mkieee(coordlist, coordieee, numcoord);
        sbits(cgrib, coordieee, iofst, 32, 0, numcoord);
        iofst = iofst + (32*numcoord);
        free(coordieee);
    }
 
    /* Calculate length of section 4 and store it in octets 1-4 of */
    /* section 4. */
    lensec4 = (iofst - ibeg) / 8;
    sbit(cgrib, &lensec4, ibeg, 32);
 
    /* Pack Data using appropriate algorithm Get Data Representation */
    /* Template */
    if (!(mapdrs = getdrstemplate(idrsnum)))
    {
        ierr = -5;
        return ierr;
    }
 
    /* Contract data field, removing data at invalid grid points, if */
    /* bit-map is provided with field. */
    if (ibmap == 0 || ibmap == 254)
    {
        pfld = malloc(ngrdpts * sizeof(g2float));
        ndpts = 0;
        for (j = 0; j < ngrdpts; j++)
        {
            if (bmap[j] == 1)
                pfld[ndpts++] = fld[j];
        }
    }
    else
    {
        ndpts = ngrdpts;
        pfld = fld;
    }
 
    /* Allocate storage for the packed data. */
    nsize = ndpts * 4;
    if (nsize < minsize)
        nsize = minsize;
    cpack = malloc(nsize);
 
    /* Call packing function based on idrsnum. */
    if (idrsnum == 0)           /*  Simple Packing */
        simpack(pfld, ndpts, idrstmpl, cpack, &lcpack);
    else if (idrsnum == 2 || idrsnum == 3)           /*  Complex Packing */
        cmplxpack(pfld, ndpts, idrsnum, idrstmpl, cpack, &lcpack);
    else if (idrsnum == 50)  /* Sperical Harmonic Simple Packing */
    { 
        simpack(pfld + 1, ndpts - 1, idrstmpl, cpack, &lcpack);
        mkieee(pfld, idrstmpl + 4, 1);  /* ensure RE(0, 0) value is IEEE format */
    }
    else if (idrsnum == 51)      /* Sperical Harmonic Complex Packing */
    {   
        getpoly(cgrib + lpos3, &JJ, &KK, &MM);
        if (JJ != 0 && KK != 0 && MM != 0)
            specpack(pfld, ndpts, JJ, KK, MM, idrstmpl, cpack, &lcpack);
        else
        {
            printf("g2_addfield: Cannot pack DRT 5.51.\n");
            return -9;
        }
    }
#if defined USE_JPEG2000 || defined USE_OPENJPEG
    else if (idrsnum == 40 || idrsnum == 40000)
    {    /*  JPEG2000 encoding  */
        if (ibmap == 255)
        {
            getdim(cgrib + lpos3, &width, &height, &iscan);
            if (width == 0 || height == 0)
            {
                width = ndpts;
                height = 1;
            }
            else if (width == allones || height == allones)
            {
                width = ndpts;
                height = 1;
            }
            else if ((iscan&32) == 32)
            {   /* Scanning mode: bit 3  */
                itemp = width;
                width = height;
                height = itemp;
            }
        }
        else
        {
            width = ndpts;
            height = 1;
        }
        lcpack = nsize;
        jpcpack(pfld, width, height, idrstmpl, cpack, &lcpack);
    }
#endif  /* USE_JPEG2000 */
#ifdef USE_PNG
    else if (idrsnum == 41 || idrsnum == 40010)
    {      /*  PNG encoding   */
        if (ibmap == 255)
        {
            getdim(cgrib + lpos3, &width, &height, &iscan);
            if (width==0 || height==0)
            {
                width = ndpts;
                height = 1;
            }
            else if (width == allones || height == allones)
            {
                width = ndpts;
                height = 1;
            }
            else if ((iscan & 32) == 32)
            {   /* Scanning mode: bit 3  */
                itemp = width;
                width = height;
                height = itemp;
            }
        }
        else
        {
            width = ndpts;
            height = 1;
        }
        pngpack(pfld, width, height, idrstmpl, cpack, &lcpack);
    }
#endif  /* USE_PNG */
    else
    {
        printf("g2_addfield: Data Representation Template 5.%ld not yet implemented.\n", idrsnum);
        ierr = -7;
        return ierr;
    }
    if (ibmap == 0 || ibmap == 254)
    {      /* free temp space */
        if (fld != pfld)
            free(pfld);
    }
    if (lcpack < 0)
    {
        if (cpack)
            free(cpack);
        ierr = -10;
        return ierr;
    }
 
    /*  Add Section 5  - Data Representation Section */
    ibeg = iofst;            /*   Calculate offset for beginning of section 5 */
    iofst = ibeg + 32;         /*   leave space for length of section */
    sbit(cgrib, &five, iofst, 8);     /* Store section number (5) */
    iofst = iofst + 8;
    sbit(cgrib, &ndpts, iofst, 32);    /* Store num of actual data points */
    iofst = iofst + 32;
    sbit(cgrib, &idrsnum, iofst, 16);    /* Store Data Repr. Template num. */
    iofst = iofst + 16;
 
    /*   Pack up each input value in array idrstmpl into the */
    /*   the appropriate number of octets, which are specified in */
    /*   corresponding entries in array mapdrs. */
    for (i = 0; i < mapdrs->maplen; i++)
    {
        nbits = abs(mapdrs->map[i]) * 8;
        if ((mapdrs->map[i] >= 0) || (idrstmpl[i] >= 0))
            sbit(cgrib, idrstmpl + i, iofst, nbits);
        else
        {
            sbit(cgrib, &one, iofst, 1);
            temp = abs(idrstmpl[i]);
            sbit(cgrib, &temp, iofst + 1, nbits - 1);
        }
        iofst = iofst + nbits;
    }
    free(mapdrs);
 
    /* Calculate length of section 5 and store it in octets */
    /* 1-4 of section 5. */
    lensec5 = (iofst - ibeg) / 8;
    sbit(cgrib, &lensec5, ibeg, 32);
 
    /* Add Section 6  - Bit-Map Section */
    ibeg = iofst;            /*   Calculate offset for beginning of section 6 */
    iofst = ibeg + 32;         /*   leave space for length of section */
    sbit(cgrib, &six, iofst, 8);     /* Store section number (6) */
    iofst = iofst + 8;
    sbit(cgrib, &ibmap, iofst, 8);    /* Store Bit Map indicator */
    iofst = iofst + 8;
 
    /*  Store bitmap, if supplied */
    if (ibmap == 0)
    {
        sbits(cgrib, bmap, iofst, 1, 0, ngrdpts);    /* Store BitMap */
        iofst = iofst + ngrdpts;
    }
 
    /*  If specifying a previously defined bit-map, make sure */
    /*  one already exists in the current GRIB message. */
    if (ibmap == 254 && !isprevbmap)
    {
        printf("g2_addfield: Requested previously defined bitmap,");
        printf(" but one does not exist in the current GRIB message.\n");
        ierr = -8;
        return ierr;
    }
 
    /* Calculate length of section 6 and store it in octets */
    /* 1-4 of section 6. Pad to end of octect, if necessary. */
    left = 8 - (iofst % 8);
    if (left != 8)
    {
        sbit(cgrib, &zero, iofst, left);     /* Pad with zeros to fill Octet */
        iofst = iofst + left;
    }
    lensec6 = (iofst-ibeg) / 8;
    sbit(cgrib, &lensec6, ibeg, 32);
 
    /* Add Section 7  - Data Section */
    ibeg = iofst;            /*   Calculate offset for beginning of section 7 */
    iofst = ibeg + 32;        /*   leave space for length of section */
    sbit(cgrib, &seven, iofst, 8);    /* Store section number (7) */
    iofst = iofst + 8;
 
    /* Store Packed Binary Data values, if non-constant field. */
    if (lcpack != 0)
    {
        ioctet = iofst / 8;
        /*cgrib(ioctet + 1:ioctet + lcpack)=cpack(1:lcpack) */
        for (j = 0; j < lcpack; j++)
            cgrib[ioctet + j] = cpack[j];
        iofst = iofst + (8*lcpack);
    }
 
    /* Calculate length of section 7 and store it in octets */
    /* 1-4 of section 7. */
    lensec7 = (iofst - ibeg) / 8;
    sbit(cgrib, &lensec7, ibeg, 32);
 
    if (cpack != 0)
        free(cpack);
 
    /*  Update current byte total of message in Section 0 */
    newlen = lencurr + lensec4 + lensec5 + lensec6 + lensec7;
    sbit(cgrib, &newlen, 96, 32);
 
    return newlen;
 
}