Skip to content

Commit 81d4e56

Browse files
authored
Improve header checker for netCDF grids (#4251)
* Clarify and improve the netcdf header checker I added more temporary variables to the neader checker so it is easies to determine what the tests are based on. Also, a netCDF file with a y-array that went from 0 to -50 caused trouble in a check for no particular reason. This version should yield less spurious warnings. * Update gmt_nc.c * Update gmt_nc.c
1 parent c3eebc7 commit 81d4e56

File tree

1 file changed

+20
-10
lines changed

1 file changed

+20
-10
lines changed

src/gmt_nc.c

+20-10
Original file line numberDiff line numberDiff line change
@@ -379,14 +379,14 @@ GMT_LOCAL void gmtnc_set_optimal_chunksize (struct GMT_CTRL *GMT, struct GMT_GRI
379379
}
380380

381381
GMT_LOCAL bool gmtnc_not_obviously_global (double *we) {
382-
/* If range is not 360 and boundaries are not 0, 180, 360 then we pass */
382+
/* If range is not 360 and boundaries are not 0, 180, 360 then we pass true */
383383
if (!gmt_M_360_range (we[0], we[1])) return true;
384384
if (!(doubleAlmostEqualZero (we[0], 0.0) || doubleAlmostEqual (we[0], -180.0))) return true;
385385
return false;
386386
}
387387

388388
GMT_LOCAL bool gmtnc_not_obviously_polar (double *se) {
389-
/* If range is not 180 and boundaries are not -90/90 then we pass */
389+
/* If range is not 180 and boundaries are not -90/90 then we pass true */
390390
if (!gmt_M_180_range (se[0], se[1])) return true;
391391
if (!(doubleAlmostEqual (se[0], -90.0) && doubleAlmostEqual (se[1], 90.0))) return true;
392392
return false;
@@ -395,7 +395,7 @@ GMT_LOCAL bool gmtnc_not_obviously_polar (double *se) {
395395
GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, char job) {
396396
int j, err, has_vector, has_range, registration, var_spacing = 0;
397397
int old_fill_mode, status;
398-
double dummy[2], *xy = NULL;
398+
double dummy[2], min, max, *xy = NULL;
399399
char dimname[GMT_GRID_UNIT_LEN80], coord[GMT_LEN8];
400400
nc_type z_type;
401401
bool save_xy_array = !strncmp (GMT->init.module_name, "grd2xyz", 7U);
@@ -572,7 +572,7 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
572572

573573
if (job == 'r') {
574574
bool set_reg = true;
575-
double dx = 0, dy = 0;
575+
double dx = 0, dy = 0, threshold = 0.0;
576576
/* Get global information */
577577
if (gmtlib_nc_get_att_text (GMT, ncid, NC_GLOBAL, "title", header->title, GMT_GRID_TITLE_LEN80))
578578
gmtlib_nc_get_att_text (GMT, ncid, z_id, "long_name", header->title, GMT_GRID_TITLE_LEN80);
@@ -636,14 +636,18 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
636636
nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_max", &dummy[1])));
637637

638638
if (has_vector && has_range) { /* Has both so we can do a basic sanity check */
639-
if (fabs (dummy[0] - xy[0]) > ((0.5+GMT_CONV5_LIMIT) * dx) || fabs (dummy[1] - xy[header->n_columns-1]) > ((0.5+GMT_CONV5_LIMIT) * dx)) {
639+
threshold = (0.5+GMT_CONV5_LIMIT) * dx;
640+
min = xy[0]; max = xy[header->n_columns-1];
641+
if (min > max) gmt_M_double_swap (min, max); /* Got it backwards in the array */
642+
if (fabs (dummy[0] - min) > threshold || fabs (dummy[1] - max) > threshold) {
640643
if (gmtnc_not_obviously_global (dummy)) {
641644
GMT_Report (GMT->parent, GMT_MSG_WARNING, "The x-coordinates and range attribute are in conflict; must rely on coordinates only\n");
642645
dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1];
643646
has_range = false; /* Since useless information */
644647
/* For registration, we have to assume that the actual range is an integer multiple of increment.
645648
* If so, then if the coordinates are off by 0.5*dx then we assume we have pixel registration */
646-
if (set_reg && fabs (fmod (dummy[0], dx)) > ((0.5-GMT_CONV5_LIMIT) * dx)) { /* Pixel registration */
649+
threshold = (0.5-GMT_CONV5_LIMIT) * dx;
650+
if (set_reg && fabs (fmod (dummy[0], dx)) > threshold) { /* Pixel registration */
647651
registration = header->registration = GMT_GRID_PIXEL_REG;
648652
dummy[0] -= 0.5 * dx; dummy[1] += 0.5 * dx;
649653
if (gmt_M_360_range (dummy[0], dummy[1]))
@@ -666,7 +670,8 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
666670
}
667671
else if (has_vector) { /* No attribute for range, use coordinates */
668672
dummy[0] = xy[0], dummy[1] = xy[header->n_columns-1];
669-
if (set_reg && fabs (fmod (dummy[0], dx)) > ((0.5-GMT_CONV5_LIMIT) * dx)) { /* Most likely pixel registration since off by dx/2 */
673+
threshold = (0.5-GMT_CONV5_LIMIT) * dx;
674+
if (set_reg && fabs (fmod (dummy[0], dx)) > threshold) { /* Most likely pixel registration since off by dx/2 */
670675
registration = header->registration = GMT_GRID_PIXEL_REG;
671676
dummy[0] -= 0.5 * dx; dummy[1] += 0.5 * dx;
672677
if (gmt_M_360_range (dummy[0], dummy[1]))
@@ -725,14 +730,18 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
725730
nc_get_att_double (ncid, ids[HH->xy_dim[1]], "valid_max", &dummy[1])));
726731

727732
if (has_vector && has_range) { /* Has both so we can do a basic sanity check */
728-
if (fabs (dummy[0] - xy[0]) > ((0.5+GMT_CONV5_LIMIT) * dy) || fabs (dummy[1] - xy[header->n_rows-1]) > ((0.5+GMT_CONV5_LIMIT) * dy)) {
733+
threshold = (0.5+GMT_CONV5_LIMIT) * dy;
734+
min = xy[0]; max = xy[header->n_rows-1];
735+
if (min > max) gmt_M_double_swap (min, max); /* Got it backwards in the array */
736+
if (fabs (dummy[0] - min) > threshold || fabs (dummy[1] - max) > threshold) {
729737
if (gmtnc_not_obviously_polar (dummy)) {
730738
GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict; must rely on coordinates only\n");
731739
dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
732740
has_range = false; /* Since useless information */
733741
/* Registration was set using x values, so here we just check that we get the same result.
734742
* If the coordinates are off by 0.5*dy then we assume we have pixel registration */
735-
if (fabs (fmod (dummy[0], dy)) > ((0.5-GMT_CONV5_LIMIT) * dy)) { /* Pixel registration? */
743+
threshold = (0.5-GMT_CONV5_LIMIT) * dy;
744+
if (fabs (fmod (dummy[0], dy)) > threshold) { /* Pixel registration? */
736745
if (header->registration == GMT_GRID_NODE_REG) /* No, somehow messed up now */
737746
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
738747
else { /* Pixel registration confirmed */
@@ -758,8 +767,9 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
758767
}
759768
}
760769
else if (has_vector) { /* No attribute for range, use coordinates */
770+
threshold = (0.5-GMT_CONV5_LIMIT) * dy;
761771
dummy[0] = xy[0], dummy[1] = xy[header->n_rows-1];
762-
if ((fabs (fmod (dummy[0], dy)) > ((0.5-GMT_CONV5_LIMIT) * dy))) { /* Most likely pixel registration since off by dy/2 */
772+
if (fabs (fmod (dummy[0], dy)) > threshold) { /* Most likely pixel registration since off by dy/2 */
763773
if (header->registration == GMT_GRID_NODE_REG) /* No, somehow messed up now */
764774
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
765775
else { /* Pixel registration confirmed */

0 commit comments

Comments
 (0)