g2_addfield.c

Go to the documentation of this file.

 
#include "grib2_int.h"
#include <stdio.h>
#include <stdlib.h>
 
g2int
g2_addfield(unsigned char *cgrib, g2int ipdsnum, g2int *ipdstmpl,
            float *coordlist, g2int numcoord, g2int idrsnum, g2int *idrstmpl,
            float *fld, g2int ngrdpts, g2int ibmap, g2int *bmap)
{
    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;
    float *pfld;
    gtemplate *mappds, *mapdrs;
#if defined USE_PNG || defined USE_JPEG2000 || defined USE_OPENJPEG || defined USE_AEC
    unsigned int allones = 4294967295u;
    g2int width, height, iscan, itemp;
#endif
    int ret;
 
    /* Check for GRIB header and terminator. Translate the error codes
     * to the legacy G2 error codes. */
    if ((ret = g2c_check_msg(cgrib, &lencurr, 1)))
    {
        if (ret == G2C_ENOTGRIB)
            return G2_ADD_MSG_INIT;
        if (ret == G2C_EMSGCOMPLETE)
            return G2_ADD_MSG_COMPLETE;
    }
 
    /* 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);
            return G2_BAD_SEC_COUNTS;
        }
    }
 
    /* 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 %ld was the last found in given GRIB message.\n",
               isecnum);
        return G2_BAD_SEC;
    }
    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");
        return G2_ADDFIELD_BAD_GDS;
    }
 
    /* 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)))
        return G2_ADDFIELD_BAD_PDT;
 
    /* 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++;
        }
    }
    if (mappds->