better handling for lat/lon (#425)

* handling dimension metadata * development of dimension code * trying dims * adding test file * adding test file * fixing doxygen * fixed warning * fixed warning * fixed warning * fixed leak * added extra testing to CMakeFiles * added extra testing to CMakeFiles * changes * fixing? * more testing * more testing * fixed warning --------- Co-authored-by: Ed <edwardjameshartnett@gmail.com>
NOAA-EMC · Aug 3, 2023 · 67a2500 · 67a2500
1 parent 7abc8f5
commit 67a2500
Show file tree

Hide file tree

Showing 12 changed files with 541 additions and 43 deletions.
diff --git a/.github/workflows/developer.yml b/.github/workflows/developer.yml
@@ -57,7 +57,7 @@ jobs:
       uses: actions/cache@v3
       with:
         path: ~/data
-        key: data-1
+        key: data-2
 
     - name: build
       run: |

diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -23,7 +23,8 @@ option(USE_OpenJPEG "Use OpenJPEG library" OFF)
 option(BUILD_SHARED_LIBS "Build shared libraries" ON)
 option(BUILD_STATIC_LIBS "Build static libraries" ON)
 option(FTP_TEST_FILES "Fetch and test with files on FTP site." OFF)
-option(LARGE_FTP_TESTS "When testing with FTP files, include very large (> 2GB) test files - this is slow!" OFF)
+option(FTP_LARGE_TEST_FILES "Fetch and test with very large files on FTP site." OFF)
+option(FTP_EXTRA_TEST_FILES "Test with files not yet available via FTP." OFF)
 option(LOGGING "Turn on internal logging messages. Only useful to g2c developers." OFF)
 option(PTHREADS "Turn on thread-safty with pthreads." OFF)
 option(UTILS "Build and install some utility programs." ON)

diff --git a/src/g2cdegrib2.c b/src/g2cdegrib2.c
@@ -611,8 +611,8 @@ g2c_degrib2(int g2cid, const char *fileout)
             char abbrev[G2C_MAX_NOAA_ABBREV_LEN + 1];
             char level_desc[G2C_MAX_TYPE_OF_FIXED_SURFACE_LEN + 1];
             char date_time[100 + 1];
-            float *data;
-            float total, max, min;
+            float *data = NULL;
+            float total = 0.0, max = 0.0, min = 0.0;
             int t;
 
             fprintf(f, "\n");
@@ -726,31 +726,39 @@ g2c_degrib2(int g2cid, const char *fileout)
 
             /* Now read the data and find the min, max, and average values. */
 
-            /* Allocate storage for the data. */
-            if (!(data = malloc(sec5_info->num_data_points * sizeof(float))))
-                return G2C_ERROR;
-
+            LOG((5, "sec5_info->num_data_points %ld", sec5_info->num_data_points));
+            if (sec5_info->num_data_points)
+            {
+                /* Allocate storage for the data. */
+                if (!(data = malloc(sec5_info->num_data_points * sizeof(float))))
+                    return G2C_ERROR;
+            }
+
             /* Get the data from message 0, product 0. */
             if ((ret = g2c_get_prod(g2cid, msg->msg_num, fld, NULL, data)))
                 return ret;
 
             /* Find min/max/avg. */
-            max = data[0];
-            min = data[0];
-            total = data[0];
-            for (i = 1; i < sec5_info->num_data_points; i++)
+            if (sec5_info->num_data_points)
             {
-                total += data[i];
-                if (data[i] > max)
-                    max = data[i];
-                if (data[i] < min)
-                    min = data[i];
+                max = data[0];
+                min = data[0];
+                total = data[0];
+                for (i = 1; i < sec5_info->num_data_points; i++)
+                {
+                    total += data[i];
+                    if (data[i] > max)
+                        max = data[i];
+                    if (data[i] < min)
+                        min = data[i];
+                }
             }
             fprintf(f, "( PARM= %s ) :  MIN=%25.8f AVE=%25.8f MAX=%25.8f\n",
                     abbrev, min, total/sec5_info->num_data_points, max);
 
             /* Free the data. */
-            free(data);
+            if (sec5_info->num_data_points)
+                free(data);
 
             total_fields++;
         }

diff --git a/src/g2cfile.c b/src/g2cfile.c
@@ -282,6 +282,67 @@ find_available_g2cid(int *g2cid)
     return ret;
 }
 
+
+/**
+ * Determine the dimension information from the section 3 metadata.
+ *
+ * See (GRIB2 - SECTION 3 GRID DEFINITION
+ * SECTION)[https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect3.shtml].
+ *
+ * For a list of grid definitions see [GRIB2 - TABLE 3.1 Grid
+ * Definition Template
+ * Number](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-1.shtml).
+ *
+ * @param sec G2C_SECTION3_INFO_T struct.
+ *
+ * @return
+ * - ::G2C_NOERROR No error.
+ * - ::G2C_EINVAL Invalid input.
+ * - ::G2C_ENOMEM Out of memory.
+ *
+ * @author Ed Hartnett @date Sep 15, 2022
+ */
+static int
+determine_dims(G2C_SECTION_INFO_T *sec)
+{
+    G2C_DIM_INFO_T *d0, *d1;
+    G2C_SECTION3_INFO_T *sec3_info;
+    int d;
+
+    sec3_info = (G2C_SECTION3_INFO_T *)(sec->sec_info);
+    d0 = &(sec3_info->dim[0]);
+    d1 = &(sec3_info->dim[1]);
+
+    /* Based on the grid definition template number, and the contents
+     * of the template, decide the len, name, and values of the two
+     * dimensions. */
+    switch (sec3_info->grid_def)
+    {
+    case 0:
+        LOG((5, "determine_dim allocating storage for lat/lon values"));
+        d0->len = sec->template[8];
+        strncpy(d0->name, LATITUDE, G2C_MAX_NAME);
+        if (!(d0->value = malloc(d0->len * sizeof(float))))
+            return G2C_ENOMEM;
+        d0->value[0] = sec->template[11];
+        for (d = 1; d < d0->len; d++)
+            d0->value[d] = d0->value[d - 1] - sec->template[16];
+
+        d1->len = sec->template[7];
+        strncpy(d1->name, LONGITUDE, G2C_MAX_NAME);
+        if (!(d1->value = malloc(d1->len * sizeof(float))))
+            return G2C_ENOMEM;
+        d1->value[0] = sec->template[12];
+        for (d = 1; d < d1->len; d++)
+            d1->value[d] = d1->value[d - 1] - sec->template[17];
+        break;
+    default:
+        break;
+    }
+
+    return G2C_NOERROR;
+}
+
 /**
  * Read the metadata from section 3 (Grid Definition Section) of a
  * GRIB2 message.
@@ -293,7 +354,7 @@ find_available_g2cid(int *g2cid)
  *
  * @param f FILE pointer to open GRIB2 file.
  * @param rw_flag ::G2C_FILE_WRITE if function should write,
- * ::G2C_FILE_READ if it should read.
+ * ::G2C_FILE_READ (0) if it should read. 
  * @param sec Pointer to the G2C_SECTION_INFO_T struct.
  *
  * @return
@@ -369,6 +430,11 @@ g2c_rw_section3_metadata(FILE *f, int rw_flag, G2C_SECTION_INFO_T *sec)
     if (!rw_flag)
         sec->sec_info = sec3_info;
 
+    /* Figure out the dimensions, if we can. */
+    if (!rw_flag)
+        if ((ret = determine_dims(sec)))
+            return ret;
+
     LOG((6, "finished reading or writing section 3 at file position %ld", ftell(f)));
     return G2C_NOERROR;
 }
@@ -1139,6 +1205,17 @@ free_metadata(int g2cid)
             stmp = sec->next;
             if (sec->template)
                 free(sec->template);
+            /* Free dim info in section 3. */
+            if (sec->sec_num == 3)
+            {
+                LOG((5, "free_metadata freeing storage for lat/lon values"));
+                float *v0 = ((G2C_SECTION3_INFO_T *)(sec->sec_info))->dim[0].value;
+                float *v1 = ((G2C_SECTION3_INFO_T *)(sec->sec_info))->dim[1].value;
+                if (v0)
+                    free(v0);
+                if (v1)
+                    free(v1);
+            }
             if (sec->sec_info)
                 free(sec->sec_info);
             free(sec);

diff --git a/src/g2cinq.c b/src/g2cinq.c
@@ -300,3 +300,83 @@ g2c_inq_prod(int g2cid, int msg_num, int prod_num, int *pds_template_len,
 
     return ret;
 }
+
+/**
+ * Learn about the one of the dimensions of a GRIB2 product.
+ *
+ * @param g2cid ID of the opened file, as from g2c_open().
+ * @param msg_num Number of the message in the file, starting with the
+ * first message as 0.
+ * @param prod_num Product number.
+ * @param dim_num Dimension number, with the first dimension as 0.
+ * @param len Pointer that gets the length of this dimension. Ignored if NULL.
+ * @param name Pointer that gets the name of this dimension. Must have
+ * memory of size G2C_MAX_NAME. Ignored if NULL.
+ * @param val Pointer that gets array of dimension values, of length
+ * len. Ignored if NULL.
+ *
+ * @return
+ * - ::G2C_NOERROR No error.
+ * - ::G2C_EBADID File ID not found.
+ * - ::G2C_ENOMSG Message not found.
+ * - ::G2C_ENOPRODUCT Product not found.
+ * - ::G2C_ENOSECTION GDS not found.
+ *
+ * @author Ed Hartnett @date 10/21/22
+ */
+int
+g2c_inq_dim(int g2cid, int msg_num, int prod_num, int dim_num, size_t *len,
+                char *name, float *val)
+{
+    G2C_MESSAGE_INFO_T *msg;
+    G2C_SECTION_INFO_T *sec4, *sec3;
+    G2C_DIM_INFO_T *dim;
+    int d;
+    int ret = G2C_NOERROR;
+
+    /* Is this an open GRIB2 file? */
+    if (g2cid < 0 || g2cid > G2C_MAX_FILES)
+        return G2C_EBADID;
+
+    /* If using threading, lock the mutex. */
+    MUTEX_LOCK(m);
+
+    if (g2c_file[g2cid].g2cid != g2cid)
+        return G2C_EBADID;
+
+    /* Find the message. */
+    for (msg = g2c_file[g2cid].msg; msg; msg = msg->next)
+        if (msg->msg_num == msg_num)
+            break;
+    if (!msg)
+        return G2C_ENOMSG;
+
+    /* Find the product. After this, sec4 will point to the
+     * appropropriate section 4 G2C_SECTION_INFO_T. */
+    for (sec4 = msg->sec; sec4; sec4 = sec4->next)
+        if (sec4->sec_num == 4 && ((G2C_SECTION4_INFO_T *)sec4->sec_info)->field_num == prod_num)
+            break;
+    if (!sec4)
+        return G2C_ENOPRODUCT;
+    /* sec4_info = (G2C_SECTION4_INFO_T *)sec4->sec_info; */
+
+    /* Find the GDS. */
+    for (sec3 = sec4->prev; sec3; sec3 = sec3->prev)
+        if (sec3->sec_num == 3)
+            break;
+    if (!sec3)
+        return G2C_ENOSECTION;
+    dim = &((G2C_SECTION3_INFO_T *)sec3->sec_info)->dim[dim_num];
+
+    /* Give the caller the info they want. */
+    if (len)
+        *len = dim->len;
+    if (name)
+        strncpy(name, dim->name, G2C_MAX_NAME);
+    if (val)
+        for (d = 0; d < dim->len; d++)
+            val[d] = dim->value[d];
+
+
+    return ret;
+}
diff --git a/src/grib2.h.in b/src/grib2.h.in
@@ -322,18 +322,8 @@ int g2c_inq_msg_time(int g2cid, int msg_num, unsigned char *sig_ref_time, short
 int g2c_inq_prod(int g2cid, int msg_num, int prod_num, int *pds_template_len,
                  int *pds_template, int *gds_template_len, int *gds_template,
                  int *drs_template_len, int *drs_template);
-
-/* Inquiry functions. */
-int g2c_inq(int g2cid, int *num_msg);
-int g2c_inq_msg(int g2cid, int msg_num, unsigned char *discipline, int *num_fields,
-                int *num_local, short *center, short *subcenter, unsigned char *master_version,
-                unsigned char *local_version);
-int g2c_inq_msg_time(int g2cid, int msg_num, unsigned char *sig_ref_time, short *year,
-                     unsigned char *month, unsigned char *day, unsigned char *hour,
-                     unsigned char *minute, unsigned char *second);
-int g2c_inq_prod(int g2cid, int msg_num, int prod_num, int *pds_template_len,
-                 int *pds_template, int *gds_template_len, int *gds_template,
-                 int *drs_template_len, int *drs_template);
+int g2c_inq_dim(int g2cid, int msg_num, int prod_num, int dim_num, size_t *len,
+                char *name, float *val);
 
 /* Getting data. */
 int g2c_get_prod(int g2cid, int msg_num, int prod_num, int *num_data_points,

diff --git a/src/grib2_int.h b/src/grib2_int.h
@@ -2,11 +2,6 @@
  * @brief Header file with internal function prototypes NCEPLIBS-g2c
  * library.
  *
- * ### Program History Log
- * Date | Programmer | Comments
- * -----|------------|---------
- * 2021-11-08 | Ed Hartnett | Initial
- *
  * @author Ed Hartnett @date 2021-11-08
  */
 
@@ -57,6 +52,12 @@
 #define FOUR_BYTES 4 /**< Four bytes. */
 #define EIGHT_BYTES 8 /**< Eight bytes. */
 
+/** Latitude */
+#define LATITUDE "Latitude"
+
+/** Longitude */
+#define LONGITUDE "Longitude"
+
 /* For thread-safety, use these macros. */
 #ifdef PTHREADS
 
@@ -155,6 +156,15 @@ typedef struct g2c_section_info
     int template_len; /**< Number of entries in template. */
 } G2C_SECTION_INFO_T;
 
+/** Keep information about dimensions defined in section 3. */
+typedef struct g2c_dim_info
+{
+    int dimid; /**< Dimension ID. */
+    size_t len; /**< Length of dimension. */
+    char name[G2C_MAX_NAME + 1]; /**< Name of dimension. */
+    float *value; /**< Array of dimension values. */
+} G2C_DIM_INFO_T;
+
 /** Information about [Section 3 GRID DEFINITION
  * SECTION](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_sect3.shtml). */
 typedef struct g2c_section3_info
@@ -165,6 +175,7 @@ typedef struct g2c_section3_info
     unsigned char interp_list; /**< Interpetation of list of numbers defining number of points (See Table 3.11). */
     unsigned short grid_def; /**< Grid definition template number (= N) (See Table 3.1). */
     int *optional; /**< Optional list of numbers defining number of points. */
+    G2C_DIM_INFO_T dim[2]; /**< Dimension information. */    
 } G2C_SECTION3_INFO_T;
 
 /** Information about [Section 4 PRODUCT DEFINITION