Skip to content

Commit

Permalink
Patch for errors in NASA OceanColor L3 nc files. Specially fix the UD…
Browse files Browse the repository at this point in the history
… flip issue. (#8380)

NASA MODIS L3 grids are a mess, with wrong coordinates and the valid_min/max attributes in opposite order of that of the lat vector. That is, the lat vector is top-down but valid_min = -90 and valid_max = 90. Furthermore, the grid is grid-registered and the start lat should be -90 + 180/4321/2 = -89.97917148808146 but it wrongly starts at -89.979171752929688 (well ends at this number because the vec is top-down). Same mess with the lon vector.
  • Loading branch information
joa-quim committed Feb 24, 2024
1 parent 1a12313 commit de507ff
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 8 deletions.
3 changes: 3 additions & 0 deletions src/gmt_gdalread.c
Original file line number Diff line number Diff line change
Expand Up @@ -446,6 +446,9 @@ GMT_LOCAL int gmtgdalread_populate_metadata (struct GMT_CTRL *GMT, struct GMT_GD

Ctrl->RasterCount = raster_count = GDALGetRasterCount(hDataset);

if (raster_count == 0)
GMT_Report(GMT->parent, GMT_MSG_ERROR, "No layers found in this file. Likely data is stored in SUBDATASETS\n.Must provide the SUBDATASET name.\n");

/* ------------------------------------------------------------------------- */
/* Get some metadata for each band. */
/* ------------------------------------------------------------------------- */
Expand Down
42 changes: 34 additions & 8 deletions src/gmt_nc.c
Original file line number Diff line number Diff line change
Expand Up @@ -688,6 +688,19 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
nc_get_att_double (ncid, ids[HH->xy_dim[0]], "valid_max", &dummy[1])));

if (has_vector && has_range) { /* Has both so we can do a basic sanity check */

if ((strstr(HH->name, ".L3m.") != NULL) && ((strstr(HH->name, "_MODIS") != NULL) || (strstr(HH->name, "_VIIRS") != NULL))) {
/* Patch for the OceanColor level 3 gids that have WRONG coordinates, and advertize a PIXEL registration
valid_min|max whilst their lon vector is GRID registered.
*/
double inc;
dummy[0] = -180.0 + 180.0 / (header->n_columns + 1); xy[0] = dummy[0];
dummy[1] = 180.0 - 180.0 / (header->n_columns + 1);
inc = gmt_M_get_inc(GMT, dummy[0], dummy[1], header->n_columns, registration);
for (unsigned k = 1; k < header->n_columns; k++)
xy[k] = xy[0] + k * inc;
}

threshold = (0.5+GMT_CONV5_LIMIT) * dx;
min = xy[0]; max = xy[header->n_columns-1];
if (min > max) gmt_M_double_swap (min, max); /* Got it backwards in the array */
Expand Down Expand Up @@ -743,13 +756,6 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
if (gmt_M_is_dnan(header->inc[GMT_X]) || gmt_M_is_zero (header->inc[GMT_X])) header->inc[GMT_X] = 1.0;
if (header->n_columns == 1) registration = GMT_GRID_PIXEL_REG; /* The only way to have a grid like that */

#ifdef NC4_DEBUG
GMT_Report (GMT->parent, GMT_MSG_WARNING, "x registration: %u\n", header->registration);
GMT_Report (GMT->parent, GMT_MSG_WARNING, "x dummy: %g %g\n", dummy[0], dummy[1]);
GMT_Report (GMT->parent, GMT_MSG_WARNING, "x[0] x[nx-1]: %g %g\n", xy[0], xy[header->n_columns-1]);
GMT_Report (GMT->parent, GMT_MSG_WARNING, "xinc: %g %g\n", header->inc[GMT_X]);
#endif

/* Extend x boundaries by half if we found pixel registration */
if (registration == GMT_GRID_NODE_REG && header->registration == GMT_GRID_PIXEL_REG)
header->wesn[XLO] = dummy[0] - header->inc[GMT_X] / 2.0,
Expand Down Expand Up @@ -786,6 +792,14 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
if (has_vector && has_range) { /* Has both so we can do a basic sanity check */
threshold = (0.5+GMT_CONV5_LIMIT) * dy;
min = xy[0]; max = xy[header->n_rows-1];

/* Need to do this here because NASA MODIS L3 grids are a mess, with wrong coordinates and the valid_min/max attributes
in opposite order of that of the lat vector. That is, the lat vector is top-down but valid_min = -90 and valid_max = 90.
Furthermore, the grid is grid-registered and the start lat should be -90 + 180/4321/2 = -89.97917148808146 but it wrongly
starts at -89.979171752929688 (well ends at this number because the vec is top-down). Same mess with the lon vector.
*/
if (min > max) HH->row_order = k_nc_start_north;

if (min > max) gmt_M_double_swap (min, max); /* Got it backwards in the array */
if (fabs (dummy[0] - min) > threshold || fabs (dummy[1] - max) > threshold) {
if (gmtnc_not_obviously_polar (dummy)) {
Expand All @@ -807,7 +821,19 @@ GMT_LOCAL int gmtnc_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *head
}
}
else {
GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict but range is exactly 180; we rely on this range\n");
if ((strstr(HH->name, ".L3m.") != NULL) && ((strstr(HH->name, "_MODIS") != NULL) || (strstr(HH->name, "_VIIRS") != NULL))) {
/* Patch for the OceanColor level 3 gids that have WRONG coordinates, and advertize a PIXEL registration
valid_min|max whilst their lat vector is GRID registered.
*/
double inc;
dummy[0] = -90.0 + 180.0 / (header->n_rows + 1); xy[0] = dummy[0];
dummy[1] = 90.0 - 180.0 / (header->n_rows + 1);
inc = gmt_M_get_inc(GMT, dummy[0], dummy[1], header->n_rows, registration);
for (unsigned int k = 1; k < header->n_rows; k++)
xy[k] = xy[0] + k * inc;
}
else
GMT_Report (GMT->parent, GMT_MSG_WARNING, "The y-coordinates and range attribute are in conflict but range is exactly 180; we rely on this range\n");
if ((header->n_rows%2) == 1 && header->registration == GMT_GRID_NODE_REG) /* Pixel registration? */
GMT_Report (GMT->parent, GMT_MSG_WARNING, "Guessing of registration in conflict between x and y, using %s\n", regtype[header->registration]);
else
Expand Down

0 comments on commit de507ff

Please sign in to comment.