Skip to content

Commit bc4eed0

Browse files
PaulWesselgithub-actions[bot]
authored andcommitted
Add special check for non-rotated global grids (#3849)
* Add special check for non-rotated global grids We usually read global grids from file and they get rotated depending on the wesn selection. However, when a global grid is passed from another environment, as GMT_IS_REFERENCE, we may be projeting that image in gmt_grd_project and if the repeated meridians fall inside the map (since no rotation has happened) they are counted twice, yielding a slightly different result. Only applies to gridline-registrered grids. This fix detects this special situation and skips one of the repeated meridians. * Handle meridian to be duplicated * Update gmt_map.c * Update gmt_project.h
1 parent 97bfe88 commit bc4eed0

File tree

4 files changed

+49
-3
lines changed

4 files changed

+49
-3
lines changed

src/gmt_map.c

Lines changed: 44 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7956,13 +7956,24 @@ int gmt_grd_project (struct GMT_CTRL *GMT, struct GMT_GRID *I, struct GMT_GRID *
79567956

79577957
O->header->z_min = FLT_MAX; O->header->z_max = -FLT_MAX; /* Min/max for out */
79587958
if (GMT->common.n.antialias) { /* Blockaverage repeat pixels, at least the first ~32767 of them... */
7959-
int n_columns = O->header->n_columns, n_rows = O->header->n_rows;
7959+
bool skip_repeat = false, duplicate_east = false;
7960+
int n_columns = O->header->n_columns, n_rows = O->header->n_rows, kase, duplicate_col;
79607961
nz = gmt_M_memory (GMT, NULL, O->header->size, short int);
79617962
/* Cannot do OPENMP yet here since it would require a reduction into an output array (nz) */
7963+
kase = gmt_whole_earth (GMT, I->header->wesn, GMT->common.R.wesn);
7964+
if (kase == 1 && I->header->registration == GMT_GRID_NODE_REG) {
7965+
/* Need to avoid giving the repeated east and west meridian values double weight when they plot inside the image.
7966+
* This is only likely to happen when external global grids are passed in via GMT_IS_REFERENCE. */
7967+
skip_repeat = true; /* Since both -180 and +180 fall inside the image, we only want to use one of then (-180) */
7968+
duplicate_east = gmt_M_is_periodic (GMT); /* Since the meridian corresponding to map west only appears once we may need to duplicate on east */
7969+
if (duplicate_east) /* Find the column in I->data that corresponds to the longitude of the left boundary of the map */
7970+
duplicate_col = gmt_M_grd_x_to_col (GMT, GMT->common.R.wesn[XLO], I->header);
7971+
}
79627972

79637973
gmt_M_row_loop (GMT, I, row_in) { /* Loop over the input grid row coordinates */
79647974
if (gmt_M_is_rect_graticule (GMT)) y_proj = y_in_proj[row_in];
79657975
gmt_M_col_loop (GMT, I, row_in, col_in, ij_in) { /* Loop over the input grid col coordinates */
7976+
if (skip_repeat && col_in == 0) continue;
79667977
if (gmt_M_is_rect_graticule (GMT))
79677978
x_proj = x_in_proj[col_in];
79687979
else if (inverse)
@@ -7991,6 +8002,38 @@ int gmt_grd_project (struct GMT_CTRL *GMT, struct GMT_GRID *I, struct GMT_GRID *
79918002
if (O->data[ij_out] < O->header->z_min) O->header->z_min = O->data[ij_out];
79928003
else if (O->data[ij_out] > O->header->z_max) O->header->z_max = O->data[ij_out];
79938004
}
8005+
if (duplicate_east) {
8006+
ij_in = ij_in - I->header->n_columns + duplicate_col; /* Rewind to the col to be duplicated */
8007+
col_in = duplicate_col;
8008+
if (gmt_M_is_rect_graticule (GMT))
8009+
x_proj = GMT->current.proj.rect[XHI];
8010+
else if (inverse)
8011+
gmt_xy_to_geo (GMT, &x_proj, &y_proj, x_in[col_in], y_in[row_in]);
8012+
else {
8013+
if (GMT->current.map.outside (GMT, x_in[col_in], y_in[row_in])) continue; /* Quite possible we are beyond the horizon */
8014+
gmt_geo_to_xy (GMT, GMT->common.R.wesn[XHI], y_in[row_in], &x_proj, &y_proj);
8015+
}
8016+
8017+
/* Here, (x_proj, y_proj) is the projected grid point. Now find nearest node on the output grid */
8018+
8019+
row_out = gmt_M_grd_y_to_row (GMT, y_proj, O->header);
8020+
if (row_out < 0 || row_out >= n_rows) continue; /* Outside our grid region */
8021+
col_out = gmt_M_grd_x_to_col (GMT, x_proj, O->header);
8022+
if (col_out < 0 || col_out >= n_columns) continue; /* Outside our grid region */
8023+
8024+
/* OK, this projected point falls inside the projected grid's rectangular domain */
8025+
8026+
ij_out = gmt_M_ijp (O->header, row_out, col_out); /* The output node */
8027+
if (nz[ij_out] == 0) O->data[ij_out] = 0.0f; /* First time, override the initial value */
8028+
if (nz[ij_out] < SHRT_MAX) { /* Avoid overflow */
8029+
O->data[ij_out] += I->data[ij_in]; /* Add up the z-sum inside this rect... */
8030+
nz[ij_out]++; /* ..and how many points there were */
8031+
}
8032+
if (gmt_M_is_fnan (O->data[ij_out])) continue;
8033+
if (O->data[ij_out] < O->header->z_min) O->header->z_min = O->data[ij_out];
8034+
else if (O->data[ij_out] > O->header->z_max) O->header->z_max = O->data[ij_out];
8035+
8036+
}
79948037
}
79958038
}
79968039

src/gmt_project.h

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,9 @@ enum GMT_enum_zdown {GMT_ZDOWN_R = 0, /* Default: Annotating radius */
140140
GMT_ZDOWN_ZP = 2, /* Annotating planetary radius - r */
141141
GMT_ZDOWN_ZR = 3}; /* Annotating given radius - r */
142142

143+
/* gmt_M_is_periodic means the east and west meridians of a global map are separated */
144+
#define gmt_M_is_periodic(C) (gmt_M_is_cylindrical (C) || gmt_M_is_misc (C))
145+
143146
/* gmt_M_is_rect_graticule means parallels and meridians are orthogonal, but does not imply linear spacing */
144147
#define gmt_M_is_rect_graticule(C) (C->current.proj.projection <= GMT_MILLER)
145148

src/psxy.c

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1215,7 +1215,7 @@ EXTERN_MSC int GMT_psxy (void *V_API, int mode, void *args) {
12151215
/* Determine if we need to worry about repeating periodic symbols */
12161216
if ((Ctrl->N.mode == PSXY_CLIP_REPEAT || Ctrl->N.mode == PSXY_NO_CLIP_REPEAT) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && gmt_M_is_geographic (GMT, GMT_IN)) {
12171217
/* Only do this for projection where west and east are split into two separate repeating boundaries */
1218-
periodic = (gmt_M_is_cylindrical (GMT) || gmt_M_is_misc (GMT));
1218+
periodic = gmt_M_is_periodic (GMT);
12191219
}
12201220
n_times = (periodic) ? 2 : 1; /* For periodic boundaries we plot each symbol twice to allow for periodic clipping */
12211221

src/psxyz.c

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -858,7 +858,7 @@ EXTERN_MSC int GMT_psxyz (void *V_API, int mode, void *args) {
858858
/* Determine if we need to worry about repeating periodic symbols */
859859
if (clip_set && (Ctrl->N.mode == PSXYZ_CLIP_REPEAT || Ctrl->N.mode == PSXYZ_NO_CLIP_REPEAT) && gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && gmt_M_is_geographic (GMT, GMT_IN)) {
860860
/* Only do this for projection where west and east are split into two separate repeating boundaries */
861-
periodic = (gmt_M_is_cylindrical (GMT) || gmt_M_is_misc (GMT));
861+
periodic = gmt_M_is_periodic (GMT);
862862
if (S.symbol == GMT_SYMBOL_GEOVECTOR) periodic = false;
863863
}
864864
n_times = (periodic) ? 2 : 1; /* For periodic boundaries we plot each symbol twice to allow for periodic clipping */

0 commit comments

Comments
 (0)