Skip to content

Make the Spilhaus projection also usable from grdimage. #8746

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
joa-quim merged 2 commits into from
Jun 11, 2025
Merged

Make the Spilhaus projection also usable from grdimage. #8746

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6a9e087
Make the Spilhaus projection also usable from grdimage.
joa-quim Jun 11, 2025
2e973c5
Test that GDAL has the minimum PROJ v9.6 to run the Spilhaus proj.
joa-quim Jun 11, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions src/gmt_init.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16456,6 +16456,7 @@ void gmt_detect_oblique_region (struct GMT_CTRL *GMT, char *file) {
if (!GMT->common.R.active[RSET]) return; /* No -R given, presumably use whole grid or image */
if (!(gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && gmt_M_180_range (GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]))) return; /* Gave -Rd or -Rg so need to probe more*/
if (gmt_M_is_azimuthal (GMT) && doubleAlmostEqual (fabs (GMT->current.proj.lat0), 90.0) && !GMT->common.R.oblique) return; /* Nothing to do */
if (GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) return; /* This is one is square */
gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Save the region we were given */

if (gmt_map_setup (GMT, GMT->common.R.wesn)) /* Set up projection */
Expand Down Expand Up @@ -18544,12 +18545,14 @@ int gmt_parse_common_options (struct GMT_CTRL *GMT, char *list, char option, cha
}
if (GMT->current.gdal_read_in.hCT_fwd) OCTDestroyCoordinateTransformation(GMT->current.gdal_read_in.hCT_fwd);
if (GMT->current.gdal_read_in.hCT_inv) OCTDestroyCoordinateTransformation(GMT->current.gdal_read_in.hCT_inv);
GMT->current.gdal_read_in.hCT_fwd = gmt_OGRCoordinateTransformation (GMT, source, dest);
GMT->current.gdal_read_in.hCT_inv = gmt_OGRCoordinateTransformation (GMT, dest, source);
GMT->current.gdal_read_in.hCT_fwd = gmt_OGRCoordinateTransformation(GMT, source, dest);
GMT->current.gdal_read_in.hCT_inv = gmt_OGRCoordinateTransformation(GMT, dest, source);
GMT->current.proj.projection = strstr(dest, "spilhaus") ? GMT_PROJ4_SPILHAUS : GMT_PROJ4_PROJS; /* Special case for spilhaus */
GMT->common.J.active = true;
if (GMT->current.gdal_read_in.hCT_fwd == NULL || GMT->current.gdal_read_in.hCT_inv == NULL)
error = 1;
else if (error && strstr(dest, "+proj=")) /* In this case it arrived here with error = 1 */
error = 0;
}
else { /* Horizontal map projection */
error += (gmt_M_check_condition (GMT, GMT->common.J.active, "Option -J given more than once\n") ||
Expand Down
24 changes: 18 additions & 6 deletions src/gmt_map.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5212,7 +5212,7 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {

*search = false;
switch (GMT->current.proj.projection_GMT) {
case GMT_LINEAR: error = gmtmap_init_linear (GMT, search); break; /* Linear transformations */
//case GMT_LINEAR: error = gmtmap_init_linear (GMT, search); break; /* Linear transformations */
case GMT_POLAR: error = gmtmap_init_polar (GMT, search); break; /* Both lon/lat are actually theta, radius */
case GMT_MERCATOR: error = gmtmap_init_merc (GMT, search); break; /* Standard Mercator projection */
case GMT_STEREO: error = gmtmap_init_stereo (GMT, search); break; /* Stereographic projection */
Expand Down Expand Up @@ -5247,8 +5247,8 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {
GMT->current.proj.scale[GMT_X] = GMT->current.proj.scale[GMT_Y] = GMT->current.proj.proj4_scl;
GMT->current.map.n_lon_nodes = 360; GMT->current.map.n_lat_nodes = 180;
if (GMT->common.R.oblique) {
gmt_proj4_fwd (GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
gmt_proj4_fwd (GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &xmax, &ymax);
gmt_proj4_fwd(GMT, GMT->common.R.wesn[XLO], GMT->common.R.wesn[YLO], &xmin, &ymin);
gmt_proj4_fwd(GMT, GMT->common.R.wesn[XHI], GMT->common.R.wesn[YHI], &xmax, &ymax);
GMT->current.map.outside = &gmtmap_rect_outside;
GMT->current.map.crossing = &gmtmap_rect_crossing;
GMT->current.map.overlap = &gmtmap_rect_overlap;
Expand All @@ -5258,7 +5258,7 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {
GMT->current.map.frame.check_side = true;
}
else {
gmtmap_xy_search (GMT, &xmin, &xmax, &ymin, &ymax, GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI], GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]);
gmtmap_xy_search(GMT, &xmin, &xmax, &ymin, &ymax, GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI], GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]);
GMT->current.map.outside = &gmtmap_wesn_outside;
GMT->current.map.crossing = &gmtmap_wesn_crossing;
GMT->current.map.overlap = &gmtmap_wesn_overlap;
Expand All @@ -5267,8 +5267,20 @@ GMT_LOCAL int gmtmap_init_polyconic (struct GMT_CTRL *GMT, bool *search) {
GMT->current.map.right_edge = &gmtmap_right_rect;
GMT->current.map.frame.horizontal = 2;
}
gmtmap_setinfo (GMT, xmin, xmax, ymin, ymax, GMT->current.proj.proj4_scl);
gmtmap_setinfo(GMT, xmin, xmax, ymin, ymax, GMT->current.proj.proj4_scl);
if (GMT->current.setting.map_frame_type & GMT_IS_FANCY) GMT->current.setting.map_frame_type = GMT_IS_PLAIN;
#if (GDAL_VERSION_MAJOR >= 3 && GDAL_VERSION_MINOR >= 9 || GDAL_VERSION_MAJOR >= 4)
if (GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) { /* Spilhaus issues many "ERROR 1: Point outside of projection domain" */
#if ((PROJ_VERSION_MAJOR < 9) || (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR < 6))
GMT_Report(GMT->parent, GMT_MSG_ERROR, "To use the Spilhaus projection you need a GDAL library that was compiled with a PROJ version >= 9.6 and you don’t have that.\n");
#endif
#ifdef WIN32
CPLSetConfigOption("CPL_LOG", "NUL");
#else
CPLSetConfigOption("CPL_LOG", "/dev/null");
#endif
}
#endif
break;
}
/* Now we only have to replace the pointers to the FWD and INV transform functions */
Expand Down Expand Up @@ -8476,7 +8488,7 @@ int gmt_grd_project (struct GMT_CTRL *GMT, struct GMT_GRID *I, struct GMT_GRID *

/* On 17-Sep-2007 the slack of GMT_CONV4_LIMIT was added to allow for round-off
errors in the grid limits. */
if (gmt_M_x_is_lon (GMT, GMT_IN) && !gmt_M_is_dnan (x_proj)) {
if (gmt_M_x_is_lon (GMT, GMT_IN) && !gmt_M_is_dnan (x_proj) && !gmt_M_is_dinf(x_proj)) { /* Spilhaus returns Inf instead of NaN */
while (x_proj < I->header->wesn[XLO] - GMT_CONV4_LIMIT) x_proj += 360.0;
while (x_proj > I->header->wesn[XHI] + GMT_CONV4_LIMIT) x_proj -= 360.0;
}
Expand Down
4 changes: 3 additions & 1 deletion src/gmt_plot.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8313,8 +8313,10 @@ char *gmt_importproj4 (struct GMT_CTRL *GMT, char *pStr, int *scale_pos, char *e
//strcat(opt_J, "0/");
}

else /* We don't return yet because we may have a +width/+scale to parse */
else { /* We don't return yet because we may have a +width/+scale to parse */
GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "This projection '%s' is not natively supported in GMT\n", prjcode);
strcat(opt_J, szProj4);
}

#if 0
bool got_a = false, got_b = false;
Expand Down
Loading