From f39348c2f117ac5be6c0b4e58a8b6b276419c4d1 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sun, 8 Jun 2025 01:32:10 +0100 Subject: [PATCH 1/4] Let us have the Spilhaus projection accessible from pscoast.c For this, it was necessary to add the -g option to pscoast and make a gap detection in gmt_geo_to_xy_line(). We can not yet do the interrupted graticule nor plot a frame. --- src/gmt_init.c | 2 +- src/gmt_map.c | 255 +++++++++++++++++++++++++++------------------- src/gmt_plot.c | 3 + src/gmt_project.h | 5 +- src/grdimage.c | 2 +- src/pscoast.c | 2 +- 6 files changed, 159 insertions(+), 110 deletions(-) diff --git a/src/gmt_init.c b/src/gmt_init.c index c207166726d..3484d3ac1b2 100644 --- a/src/gmt_init.c +++ b/src/gmt_init.c @@ -18546,7 +18546,7 @@ int gmt_parse_common_options (struct GMT_CTRL *GMT, char *list, char option, cha 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.proj.projection = GMT_PROJ4_PROJS; /* This now make it use the proj4 lib */ + 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; diff --git a/src/gmt_map.c b/src/gmt_map.c index 12cff93520e..81ff5e7b71a 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -8026,43 +8026,63 @@ uint64_t gmt_cart_to_xy_line (struct GMT_CTRL *GMT, double *x, double *y, uint64 } /*! . */ -uint64_t gmt_geo_to_xy_line (struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t n) { +uint64_t gmt_geo_to_xy_line(struct GMT_CTRL *GMT, double *lon, double *lat, uint64_t n) { /* Traces the lon/lat array and returns x,y plus appropriate pen moves * Pen moves are caused by breakthroughs of the map boundary or when * a point has lon = NaN or lat = NaN (this means "pick up pen") */ uint64_t j, k, np, n_sections; - bool this_inside, jump, cartesian = gmt_M_is_cartesian (GMT, GMT_IN); + bool this_inside, jump, find_coast_gaps = false, cartesian = gmt_M_is_cartesian(GMT, GMT_IN); /* bool last_inside = false; */ unsigned int sides[4]; unsigned int nx; double xlon[4], xlat[4], xx[4], yy[4]; - double this_x, this_y, last_x, last_y, dummy[4]; + double this_x, this_y, last_x, last_y, dummy[4], dist2_pt; if (n == 0) return 0; /* Absolutely nothing to do */ - while (n > GMT->current.plot.n_alloc) gmt_get_plot_array (GMT); + while (n > GMT->current.plot.n_alloc) gmt_get_plot_array(GMT); /* Here we know n is at least 1 */ np = 0; - gmt_geo_to_xy (GMT, lon[0], lat[0], &last_x, &last_y); - if (!gmt_map_outside (GMT, lon[0], lat[0])) { /* First point is inside the region */ + gmt_geo_to_xy(GMT, lon[0], lat[0], &last_x, &last_y); + if (!gmt_map_outside(GMT, lon[0], lat[0])) { /* First point is inside the region */ GMT->current.plot.x[0] = last_x; GMT->current.plot.y[0] = last_y; GMT->current.plot.pen[np++] = PSL_MOVE; /* last_inside = true; */ } + + if (GMT->common.g.active && strstr(GMT->init.module_name, "coast")) { /* For inter-point gaps detection in pscoast data. */ + dist2_pt = pow(GMT->common.g.gap[0], 2.0); /* Squared gap in plot units */ + find_coast_gaps = true; + } + for (j = 1; j < n; j++) { gmt_geo_to_xy (GMT, lon[j], lat[j], &this_x, &this_y); /* This guy also checks if lon,lat are NaNs */ - this_inside = !gmt_map_outside (GMT, lon[j], lat[j]); /* This also calls gmt_geo_to_xy */ - if (gmt_M_is_dnan (lon[j]) || gmt_M_is_dnan (lat[j])) continue; /* Skip NaN point now. WHY NOT MOVE THIS LINE 2 LINES ABOVE? JL. */ - if (gmt_M_is_dnan (lon[j-1]) || gmt_M_is_dnan (lat[j-1])) { /* Point after NaN needs a move */ + this_inside = !gmt_map_outside(GMT, lon[j], lat[j]); /* This also calls gmt_geo_to_xy */ + if (gmt_M_is_dnan(lon[j]) || gmt_M_is_dnan(lat[j])) continue; /* Skip NaN point now. WHY NOT MOVE THIS LINE 2 LINES ABOVE? JL. */ + if (gmt_M_is_dnan(lon[j-1]) || gmt_M_is_dnan(lat[j-1])) { /* Point after NaN needs a move */ GMT->current.plot.x[np] = this_x; GMT->current.plot.y[np] = this_y; GMT->current.plot.pen[np++] = PSL_MOVE; - if (np == GMT->current.plot.n_alloc) gmt_get_plot_array (GMT); + if (np == GMT->current.plot.n_alloc) gmt_get_plot_array(GMT); last_x = this_x; last_y = this_y; /* last_inside = this_inside; */ continue; } + + if (find_coast_gaps) { /* Pen up if we are in a projected pscoast gap */ + double dx, dy; + dx = this_x - last_x; dy = this_y - last_y; + if ((dx * dx + dy * dy) > dist2_pt) { + GMT->current.plot.x[np] = this_x; GMT->current.plot.y[np] = this_y; + GMT->current.plot.pen[np++] = PSL_MOVE; + if (np == GMT->current.plot.n_alloc) gmt_get_plot_array(GMT); + last_x = this_x; last_y = this_y; + continue; + } + } + if ((nx = gmtmap_crossing (GMT, lon[j-1], lat[j-1], lon[j], lat[j], xlon, xlat, xx, yy, sides))) { /* Do nothing if we get crossings*/ } else if (GMT->current.map.is_world) /* Check global wrapping if 360 range */ nx = (*GMT->current.map.wrap_around_check) (GMT, dummy, last_x, last_y, this_x, this_y, xx, yy, sides); + if (nx == 1) { /* inside-outside or outside-inside; set move&clip vs draw flags */ GMT->current.plot.x[np] = xx[0]; GMT->current.plot.y[np] = yy[0]; /* If next point is inside then we want to move to the crossing, otherwise we want to draw to the crossing */ @@ -8082,6 +8102,7 @@ uint64_t gmt_geo_to_xy_line (struct GMT_CTRL *GMT, double *lon, double *lat, uin if (np == GMT->current.plot.n_alloc) gmt_get_plot_array (GMT); } } + if (this_inside) { if ( np >= GMT->current.plot.n_alloc ) { GMT_Report (GMT->parent, GMT_MSG_ERROR, "bad access: cannot access current.plot.x[%" PRIu64 "], np=%" PRIu64 ", GMT->current.plot.n=%" PRIu64 "\n", np, np, GMT->current.plot.n); @@ -8094,6 +8115,7 @@ uint64_t gmt_geo_to_xy_line (struct GMT_CTRL *GMT, double *lon, double *lat, uin } last_x = this_x; last_y = this_y; /* last_inside = this_inside; */ } + if (np) GMT->current.plot.pen[0] |= PSL_MOVE; /* Sanity override: Gotta start off with new start point */ /* When a line that starts and ends inside the domain exits and reenters, we end up with two pieces. @@ -8104,7 +8126,9 @@ uint64_t gmt_geo_to_xy_line (struct GMT_CTRL *GMT, double *lon, double *lat, uin n_sections++; if (n_sections == 2) k = j; /* Start of 2nd section */ } - if (n_sections == 2 && doubleAlmostEqualZero (GMT->current.plot.x[0], GMT->current.plot.x[np-1]) && doubleAlmostEqualZero (GMT->current.plot.y[0], GMT->current.plot.y[np-1])) { + if (n_sections == 2 && doubleAlmostEqualZero(GMT->current.plot.x[0], GMT->current.plot.x[np-1]) && + doubleAlmostEqualZero(GMT->current.plot.y[0], GMT->current.plot.y[np-1])) { + double *tmp = gmt_M_memory (GMT, NULL, np, double); /* Shuffle x-array safely */ gmt_M_memcpy (tmp, &GMT->current.plot.x[k], np-k, double); @@ -8748,7 +8772,7 @@ double gmt_azim_to_angle (struct GMT_CTRL *GMT, double lon, double lat, double c } /*! . */ -uint64_t gmt_map_clip_path (struct GMT_CTRL *GMT, double **x, double **y, bool *donut) { +uint64_t gmt_map_clip_path(struct GMT_CTRL *GMT, double **x, double **y, bool *donut) { /* This function returns a clip path corresponding to the * extent of the map. */ @@ -8759,6 +8783,9 @@ uint64_t gmt_map_clip_path (struct GMT_CTRL *GMT, double **x, double **y, bool * *donut = false; + if (GMT->current.proj.projection_GMT == -1 && GMT->current.proj.is_proj4) + GMT->current.proj.projection_GMT = GMT->current.proj.projection; + if (GMT->common.R.oblique) /* Rectangular map boundary */ np = 4; else { @@ -8814,6 +8841,9 @@ uint64_t gmt_map_clip_path (struct GMT_CTRL *GMT, double **x, double **y, bool * case GMT_POLYCONIC: np = 2 * (GMT->current.map.n_lon_nodes + GMT->current.map.n_lat_nodes); break; + case GMT_PROJ4_SPILHAUS: + np = 4; + break; default: GMT_Report (GMT->parent, GMT_MSG_ERROR, "Bad case in gmt_map_clip_path (%d)\n", GMT->current.proj.projection_GMT); np = 0; @@ -8974,6 +9004,17 @@ uint64_t gmt_map_clip_path (struct GMT_CTRL *GMT, double **x, double **y, bool * for (i = 1; GMT->common.R.wesn[YLO] != -90.0 && i < GMT->current.map.n_lon_nodes; i++, j++) gmt_geo_to_xy (GMT, GMT->common.R.wesn[XLO] + i * GMT->current.map.dlon, GMT->common.R.wesn[YLO], &work_x[j], &work_y[j]); break; + case GMT_PROJ4_SPILHAUS: + /* The Spilhaus proj is f. The points bellow were obtained by projecting a global grid with gdalwarp and + inverting the UL and LR coordinates. But even those had to be moved a bit to inside (by trial en error). + Couldn't do the same to the other two corners because permanently got Inf's. Anyway, for porposes of + map clipping, this seems good enough. + */ + gmt_geo_to_xy(GMT, -113.0447807067042, 49.56674656682158, &work_x[3], &work_y[3]); /* UL */ + gmt_geo_to_xy(GMT, -113.06804976730443, 49.553963054590234, &work_x[1], &work_y[1]); /* LR */ + work_x[0] = work_x[3]; work_x[2] = work_x[1]; + work_y[0] = work_y[1]; work_y[2] = work_y[3]; + break; } } @@ -9679,131 +9720,133 @@ int gmt_proj_setup (struct GMT_CTRL *GMT, double wesn[]) { gmtmap_reset_oblique_settings (GMT); /* Maybe reset MAP_ANNOT_OBLIQUE defaults depending on situation */ - switch (GMT->current.proj.projection) { - - case GMT_LINEAR: /* Linear transformations */ - error = gmtmap_init_linear (GMT, &search); - break; + if (gmt_M_is_proj4(GMT)) { /* Must process this case first. */ + error = map_init_proj4(GMT, &search); + } + else { + switch (GMT->current.proj.projection) { - case GMT_POLAR: /* Both lon/lat are actually theta, radius */ - error = gmtmap_init_polar (GMT, &search); - break; + case GMT_LINEAR: /* Linear transformations */ + error = gmtmap_init_linear (GMT, &search); + break; - case GMT_MERCATOR: /* Standard Mercator projection */ - error = gmtmap_init_merc (GMT, &search); - break; + case GMT_POLAR: /* Both lon/lat are actually theta, radius */ + error = gmtmap_init_polar (GMT, &search); + break; - case GMT_STEREO: /* Stereographic projection */ - error = gmtmap_init_stereo (GMT, &search); - break; + case GMT_MERCATOR: /* Standard Mercator projection */ + error = gmtmap_init_merc (GMT, &search); + break; - case GMT_LAMBERT: /* Lambert Conformal Conic */ - error = gmtmap_init_lambert (GMT, &search); - break; + case GMT_STEREO: /* Stereographic projection */ + error = gmtmap_init_stereo (GMT, &search); + break; - case GMT_OBLIQUE_MERC: /* Oblique Mercator */ - error = gmtmap_init_oblique (GMT, &search); - break; + case GMT_LAMBERT: /* Lambert Conformal Conic */ + error = gmtmap_init_lambert (GMT, &search); + break; - case GMT_TM: /* Transverse Mercator */ - error = gmtmap_init_tm (GMT, &search); - break; + case GMT_OBLIQUE_MERC: /* Oblique Mercator */ + error = gmtmap_init_oblique (GMT, &search); + break; - case GMT_UTM: /* Universal Transverse Mercator */ - error = gmtmap_init_utm (GMT, &search); - break; + case GMT_TM: /* Transverse Mercator */ + error = gmtmap_init_tm (GMT, &search); + break; - case GMT_LAMB_AZ_EQ: /* Lambert Azimuthal Equal-Area */ - error = gmtmap_init_lambeq (GMT, &search); - break; + case GMT_UTM: /* Universal Transverse Mercator */ + error = gmtmap_init_utm (GMT, &search); + break; - case GMT_ORTHO: /* Orthographic Projection */ - error = gmtmap_init_ortho (GMT, &search); - break; + case GMT_LAMB_AZ_EQ: /* Lambert Azimuthal Equal-Area */ + error = gmtmap_init_lambeq (GMT, &search); + break; - case GMT_GENPER: /* General Perspective Projection */ - error = gmtmap_init_genper (GMT, &search); - break; + case GMT_ORTHO: /* Orthographic Projection */ + error = gmtmap_init_ortho (GMT, &search); + break; - case GMT_AZ_EQDIST: /* Azimuthal Equal-Distance Projection */ - error = gmtmap_init_azeqdist (GMT, &search); - break; + case GMT_GENPER: /* General Perspective Projection */ + error = gmtmap_init_genper (GMT, &search); + break; - case GMT_GNOMONIC: /* Azimuthal Gnomonic Projection */ - error = gmtmap_init_gnomonic (GMT, &search); - break; + case GMT_AZ_EQDIST: /* Azimuthal Equal-Distance Projection */ + error = gmtmap_init_azeqdist (GMT, &search); + break; - case GMT_MOLLWEIDE: /* Mollweide Equal-Area */ - error = gmtmap_init_mollweide (GMT, &search); - break; + case GMT_GNOMONIC: /* Azimuthal Gnomonic Projection */ + error = gmtmap_init_gnomonic (GMT, &search); + break; - case GMT_HAMMER: /* Hammer-Aitoff Equal-Area */ - error = gmtmap_init_hammer (GMT, &search); - break; + case GMT_MOLLWEIDE: /* Mollweide Equal-Area */ + error = gmtmap_init_mollweide (GMT, &search); + break; - case GMT_VANGRINTEN: /* Van der Grinten */ - error = gmtmap_init_grinten (GMT, &search); - break; + case GMT_HAMMER: /* Hammer-Aitoff Equal-Area */ + error = gmtmap_init_hammer (GMT, &search); + break; - case GMT_WINKEL: /* Winkel Tripel */ - error = gmtmap_init_winkel (GMT, &search); - break; + case GMT_VANGRINTEN: /* Van der Grinten */ + error = gmtmap_init_grinten (GMT, &search); + break; - case GMT_ECKERT4: /* Eckert IV */ - error = gmtmap_init_eckert4 (GMT, &search); - break; + case GMT_WINKEL: /* Winkel Tripel */ + error = gmtmap_init_winkel (GMT, &search); + break; - case GMT_ECKERT6: /* Eckert VI */ - error = gmtmap_init_eckert6 (GMT, &search); - break; + case GMT_ECKERT4: /* Eckert IV */ + error = gmtmap_init_eckert4 (GMT, &search); + break; - case GMT_CYL_EQ: /* Cylindrical Equal-Area */ - error = gmtmap_init_cyleq (GMT, &search); - break; + case GMT_ECKERT6: /* Eckert VI */ + error = gmtmap_init_eckert6 (GMT, &search); + break; - case GMT_CYL_STEREO: /* Cylindrical Stereographic */ - error = gmtmap_init_cylstereo (GMT, &search); - break; + case GMT_CYL_EQ: /* Cylindrical Equal-Area */ + error = gmtmap_init_cyleq (GMT, &search); + break; - case GMT_MILLER: /* Miller Cylindrical */ - error = gmtmap_init_miller (GMT, &search); - break; + case GMT_CYL_STEREO: /* Cylindrical Stereographic */ + error = gmtmap_init_cylstereo (GMT, &search); + break; - case GMT_CYL_EQDIST: /* Cylindrical Equidistant */ - error = gmtmap_init_cyleqdist (GMT, &search); - break; + case GMT_MILLER: /* Miller Cylindrical */ + error = gmtmap_init_miller (GMT, &search); + break; - case GMT_ROBINSON: /* Robinson */ - error = gmtmap_init_robinson (GMT, &search); - break; + case GMT_CYL_EQDIST: /* Cylindrical Equidistant */ + error = gmtmap_init_cyleqdist (GMT, &search); + break; - case GMT_SINUSOIDAL: /* Sinusoidal Equal-Area */ - error = gmtmap_init_sinusoidal (GMT, &search); - break; + case GMT_ROBINSON: /* Robinson */ + error = gmtmap_init_robinson (GMT, &search); + break; - case GMT_CASSINI: /* Cassini cylindrical */ - error = gmtmap_init_cassini (GMT, &search); - break; + case GMT_SINUSOIDAL: /* Sinusoidal Equal-Area */ + error = gmtmap_init_sinusoidal (GMT, &search); + break; - case GMT_ALBERS: /* Albers Equal-Area Conic */ - error = gmtmap_init_albers (GMT, &search); - break; + case GMT_CASSINI: /* Cassini cylindrical */ + error = gmtmap_init_cassini (GMT, &search); + break; - case GMT_ECONIC: /* Equidistant Conic */ - error = gmtmap_init_econic (GMT, &search); - break; + case GMT_ALBERS: /* Albers Equal-Area Conic */ + error = gmtmap_init_albers (GMT, &search); + break; - case GMT_POLYCONIC: /* Polyconic */ - error = gmtmap_init_polyconic (GMT, &search); - break; + case GMT_ECONIC: /* Equidistant Conic */ + error = gmtmap_init_econic (GMT, &search); + break; - case GMT_PROJ4_PROJS: /* All proj.4 projections */ - error = map_init_proj4 (GMT, &search); - break; + case GMT_POLYCONIC: /* Polyconic */ + error = gmtmap_init_polyconic (GMT, &search); + break; - default: /* No projection selected, return to a horrible death */ - Error_and_return (GMT_MAP_NO_PROJECTION, GMT_PROJECTION_ERROR); + default: /* No projection selected, return to a horrible death */ + Error_and_return (GMT_MAP_NO_PROJECTION, GMT_PROJECTION_ERROR); + } } + if (error) return (GMT_PROJECTION_ERROR); /* Something went wrong in one of the map_init_* functions */ if (GMT->current.proj.fwd == NULL) /* Some error in projection projection parameters, return to a horrible death */ diff --git a/src/gmt_plot.c b/src/gmt_plot.c index e2e5e4e2463..c08f23a5a3f 100644 --- a/src/gmt_plot.c +++ b/src/gmt_plot.c @@ -3088,6 +3088,9 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) { case GMT_VANGRINTEN: gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n); break; + case GMT_PROJ4_SPILHAUS: + gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n); + break; } PSL_comment (PSL, "End of map frame\n"); } diff --git a/src/gmt_project.h b/src/gmt_project.h index 677757f3edd..f17d03f74fc 100644 --- a/src/gmt_project.h +++ b/src/gmt_project.h @@ -120,7 +120,10 @@ enum gmt_enum_misc {GMT_MOLLWEIDE = 400, GMT_WINKEL}; /* All projections from proj.4 lib */ -enum gmt_enum_allprojs {GMT_PROJ4_PROJS = 500}; +#define gmt_M_is_proj4(C) (C->current.proj.projection / 100 == 5) +enum gmt_enum_allprojs {GMT_PROJ4_PROJS = 500, + GMT_PROJ4_SPILHAUS, +}; /*! The various GMT measurement units */ enum gmt_enum_units {GMT_IS_METER = 0, diff --git a/src/grdimage.c b/src/grdimage.c index b78a723c8c4..70f08629dfe 100644 --- a/src/grdimage.c +++ b/src/grdimage.c @@ -1518,7 +1518,7 @@ EXTERN_MSC int GMT_grdimage(void *V_API, int mode, void *args) { /* Determine if grid/image is to be projected */ need_to_project = (gmt_M_is_nonlinear_graticule (GMT) || Ctrl->E.dpi > 0); - if (need_to_project && GMT->current.proj.projection == GMT_PROJ4_PROJS) { + if (need_to_project && gmt_M_is_proj4(GMT)) { if (GMT->current.proj.is_proj4) if (strstr(GMT->common.J.proj4string, "longlat") || strstr(GMT->common.J.proj4string, "lonlat") || strstr(GMT->common.J.proj4string, "latlon")) need_to_project = false; diff --git a/src/pscoast.c b/src/pscoast.c index d5f6ce0b8aa..768b5dcaa63 100644 --- a/src/pscoast.c +++ b/src/pscoast.c @@ -59,7 +59,7 @@ #define THIS_MODULE_PURPOSE "Plot continents, countries, shorelines, rivers, and borders" #define THIS_MODULE_KEYS ">?},>DE-lL" #define THIS_MODULE_NEEDS "JR" -#define THIS_MODULE_OPTIONS "->BJKOPRUVXYbdptxy" GMT_OPT("Zc") +#define THIS_MODULE_OPTIONS "->BJKOPRUVXYbdgptxy" GMT_OPT("Zc") #define LAKE 0 #define RIVER 1 From 9d2c6532ec7d907404afda19d6b4837e90ae4374 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sun, 8 Jun 2025 12:43:27 +0100 Subject: [PATCH 2/4] Use a default -g in pscoast when proj is Spilhaus --- src/gmt_plot.c | 9 ++++++--- src/pscoast.c | 6 ++++++ 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/gmt_plot.c b/src/gmt_plot.c index c08f23a5a3f..4fe31a9239e 100644 --- a/src/gmt_plot.c +++ b/src/gmt_plot.c @@ -3029,7 +3029,8 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) { if ((way = gmtlib_adjust_we_if_central_lon_set (GMT, &w, &e))) GMT_Report (GMT->parent, GMT_MSG_DEBUG, "W/E boundaries shifted by %d\n", way * 360); - PSL_comment (PSL, "Start of map frame\n"); + if (!PSL->internal.comments) PSL_command(PSL, "\n%%Start of map frame\n"); + PSL_comment(PSL, "Start of map frame\n"); gmt_setpen (GMT, &GMT->current.setting.map_frame_pen); PSL_setcolor (PSL, GMT->current.setting.map_frame_pen.rgb, PSL_IS_STROKE); @@ -3089,10 +3090,12 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) { gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n); break; case GMT_PROJ4_SPILHAUS: - gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n); + //gmtplot_linear_map_boundary(GMT, PSL, w, e, s, n); break; } - PSL_comment (PSL, "End of map frame\n"); + + PSL_comment(PSL, "End of map frame\n"); + if (!PSL->internal.comments) PSL_command(PSL, "\n%%End of map fram\n"); } /* gmt_map_basemap will create a basemap for the given area. diff --git a/src/pscoast.c b/src/pscoast.c index 768b5dcaa63..fb0dff977ff 100644 --- a/src/pscoast.c +++ b/src/pscoast.c @@ -533,6 +533,12 @@ static int parse (struct GMT_CTRL *GMT, struct PSCOAST_CTRL *Ctrl, struct GMT_OP } } + /* If Spilhaus projection and no -g used, set a default value. */ + if (!GMT->common.g.active && GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) { + gmt_parse_g_option(GMT, "D1"); + GMT->common.g.active = true; + } + if ((error = gmt_DCW_list (GMT, &(Ctrl->E.info)))) { /* This is either success or failure... */ if (error != GMT_DCW_LIST) return (1); /* Not good */ From 931f8017957753f269919f29be0cf0e83f10225d Mon Sep 17 00:00:00 2001 From: Joaquim Date: Sun, 8 Jun 2025 23:53:03 +0100 Subject: [PATCH 3/4] Fix issues with scale and frame drawing. Now we can use -B --- src/gmt_map.c | 7 +++++++ src/gmt_plot.c | 2 +- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/src/gmt_map.c b/src/gmt_map.c index 81ff5e7b71a..960ea484603 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -2187,6 +2187,13 @@ GMT_LOCAL void gmtmap_xy_search (struct GMT_CTRL *GMT, double *x0, double *x1, d /* Find min/max forward values */ + if (GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) { /* The Spilhaus proj is very tricky. */ + (*GMT->current.proj.fwd) (GMT, -113.0447807067042, 49.56674656682158, &xmin, &ymax); /* UL */ + (*GMT->current.proj.fwd) (GMT, -113.06804976730443, 49.553963054590234, &xmax, &ymin); /* LR */ + *x0 = xmin; *x1 = xmax; *y0 = ymin; *y1 = ymax; + return; + } + xmax = ymax = -DBL_MAX; xmin = ymin = DBL_MAX; dlon = fabs (e0 - w0) / GMTMAP_N_STEPS; diff --git a/src/gmt_plot.c b/src/gmt_plot.c index 4fe31a9239e..1b0309d546e 100644 --- a/src/gmt_plot.c +++ b/src/gmt_plot.c @@ -3090,7 +3090,7 @@ GMT_LOCAL void gmtplot_map_boundary (struct GMT_CTRL *GMT) { gmtplot_basic_map_boundary (GMT, PSL, w, e, s, n); break; case GMT_PROJ4_SPILHAUS: - //gmtplot_linear_map_boundary(GMT, PSL, w, e, s, n); + gmtplot_linear_map_boundary(GMT, PSL, w, e, s, n); break; } From 205a13e145fcb8929beda298b6e874cb1a7540f6 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Mon, 9 Jun 2025 15:51:09 +0100 Subject: [PATCH 4/4] Document -g in pscoast. Make Spilhaus parameters as #defines. Shut up (avoid them) PROJ warnings. --- doc/rst/source/coast.rst | 10 +++++++++ doc/rst/source/pscoast.rst | 1 + src/gmt_map.c | 43 ++++++++++++++++++++++++-------------- src/gmt_proj.c | 14 +++++++++++++ src/pscoast.c | 9 ++++++-- 5 files changed, 59 insertions(+), 18 deletions(-) diff --git a/doc/rst/source/coast.rst b/doc/rst/source/coast.rst index 32568b22fff..133e0e9af48 100644 --- a/doc/rst/source/coast.rst +++ b/doc/rst/source/coast.rst @@ -36,6 +36,7 @@ Synopsis [ |SYN_OPT-Y| ] [ |SYN_OPT-bo| ] [ |SYN_OPT-d| ] +[ |SYN_OPT-g| ] [ |SYN_OPT-p| ] [ |SYN_OPT-t| ] [ |SYN_OPT--| ] @@ -290,6 +291,15 @@ Optional Arguments .. |Add_-d| unicode:: 0x20 .. just an invisible code .. include:: explain_-d.rst_ +.. _-g: + +**-gD**\ *dist* + Short version of the global **-g** option. Here it is used to set the _dist_ distance, in map units, + higher then which we consider to have a gap. Useful for the Spilhaus projection (though we automatically set it) + and when line wrapping on dateline was not correctly detected. As an example of this, see for instance + the old issue `#3667 `_ that can be fixed + by adding **-gD**\ 1 to the command line. + .. |Add_perspective| unicode:: 0x20 .. just an invisible code .. include:: explain_perspective.rst_ diff --git a/doc/rst/source/pscoast.rst b/doc/rst/source/pscoast.rst index 654a74e5988..854d36dd20d 100644 --- a/doc/rst/source/pscoast.rst +++ b/doc/rst/source/pscoast.rst @@ -37,6 +37,7 @@ Synopsis [ |SYN_OPT-X| ] [ |SYN_OPT-Y| ] [ |SYN_OPT-bo| ] +[ |SYN_OPT-g| ] [ |SYN_OPT-p| ] [ |SYN_OPT-t| ] [ |SYN_OPT--| ] diff --git a/src/gmt_map.c b/src/gmt_map.c index 960ea484603..a88e358708f 100644 --- a/src/gmt_map.c +++ b/src/gmt_map.c @@ -2188,8 +2188,8 @@ GMT_LOCAL void gmtmap_xy_search (struct GMT_CTRL *GMT, double *x0, double *x1, d /* Find min/max forward values */ if (GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) { /* The Spilhaus proj is very tricky. */ - (*GMT->current.proj.fwd) (GMT, -113.0447807067042, 49.56674656682158, &xmin, &ymax); /* UL */ - (*GMT->current.proj.fwd) (GMT, -113.06804976730443, 49.553963054590234, &xmax, &ymin); /* LR */ + (*GMT->current.proj.fwd) (GMT, GMT_PROJ_SPILH_LON_UL, GMT_PROJ_SPILH_LAT_UL, &xmin, &ymax); /* UL */ + (*GMT->current.proj.fwd) (GMT, GMT_PROJ_SPILH_LON_LR, GMT_PROJ_SPILH_LAT_LR, &xmax, &ymin); /* LR */ *x0 = xmin; *x1 = xmax; *y0 = ymin; *y1 = ymax; return; } @@ -6387,18 +6387,29 @@ GMT_LOCAL int gmtmap_init_three_D (struct GMT_CTRL *GMT) { xx[0] = xx[3] = GMT->current.proj.rect[XLO]; xx[1] = xx[2] = GMT->current.proj.rect[XHI]; yy[0] = yy[1] = GMT->current.proj.rect[YLO]; yy[2] = yy[3] = GMT->current.proj.rect[YHI]; - for (i = 0; i < 4; i++) { - gmt_xy_to_geo (GMT, &GMT->current.proj.z_project.corner_x[i], &GMT->current.proj.z_project.corner_y[i], xx[i], yy[i]); - gmt_xyz_to_xy (GMT, xx[i], yy[i], gmt_z_to_zz(GMT, GMT->common.R.wesn[ZLO]), &x, &y); - GMT->current.proj.z_project.xmin = MIN (GMT->current.proj.z_project.xmin, x); - GMT->current.proj.z_project.xmax = MAX (GMT->current.proj.z_project.xmax, x); - GMT->current.proj.z_project.ymin = MIN (GMT->current.proj.z_project.ymin, y); - GMT->current.proj.z_project.ymax = MAX (GMT->current.proj.z_project.ymax, y); - gmt_xyz_to_xy (GMT, xx[i], yy[i], gmt_z_to_zz(GMT, GMT->common.R.wesn[ZHI]), &x, &y); - GMT->current.proj.z_project.xmin = MIN (GMT->current.proj.z_project.xmin, x); - GMT->current.proj.z_project.xmax = MAX (GMT->current.proj.z_project.xmax, x); - GMT->current.proj.z_project.ymin = MIN (GMT->current.proj.z_project.ymin, y); - GMT->current.proj.z_project.ymax = MAX (GMT->current.proj.z_project.ymax, y); + if (GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) { + /* Again, Spilhaus is not invertible at all 4 corners, so we set the -Rd global values to avoid annoying warnings */ + GMT->current.proj.z_project.corner_x[0] = GMT->current.proj.z_project.corner_x[3] = GMT_PROJ_SPILH_LON_UL; + GMT->current.proj.z_project.corner_x[1] = GMT->current.proj.z_project.corner_x[2] = GMT_PROJ_SPILH_LON_LR; + GMT->current.proj.z_project.corner_y[0] = GMT->current.proj.z_project.corner_y[1] = GMT_PROJ_SPILH_LAT_LR; + GMT->current.proj.z_project.corner_y[2] = GMT->current.proj.z_project.corner_y[3] = GMT_PROJ_SPILH_LAT_UL; + GMT->current.proj.z_project.xmin = GMT->current.proj.z_project.ymin = 0.0; + GMT->current.proj.z_project.xmax = GMT->current.proj.z_project.ymax = GMT->current.proj.rect[XHI]; + } + else { + for (i = 0; i < 4; i++) { + gmt_xy_to_geo (GMT, &GMT->current.proj.z_project.corner_x[i], &GMT->current.proj.z_project.corner_y[i], xx[i], yy[i]); + gmt_xyz_to_xy (GMT, xx[i], yy[i], gmt_z_to_zz(GMT, GMT->common.R.wesn[ZLO]), &x, &y); + GMT->current.proj.z_project.xmin = MIN (GMT->current.proj.z_project.xmin, x); + GMT->current.proj.z_project.xmax = MAX (GMT->current.proj.z_project.xmax, x); + GMT->current.proj.z_project.ymin = MIN (GMT->current.proj.z_project.ymin, y); + GMT->current.proj.z_project.ymax = MAX (GMT->current.proj.z_project.ymax, y); + gmt_xyz_to_xy (GMT, xx[i], yy[i], gmt_z_to_zz(GMT, GMT->common.R.wesn[ZHI]), &x, &y); + GMT->current.proj.z_project.xmin = MIN (GMT->current.proj.z_project.xmin, x); + GMT->current.proj.z_project.xmax = MAX (GMT->current.proj.z_project.xmax, x); + GMT->current.proj.z_project.ymin = MIN (GMT->current.proj.z_project.ymin, y); + GMT->current.proj.z_project.ymax = MAX (GMT->current.proj.z_project.ymax, y); + } } } else if (GMT->current.proj.r > 0.0) { /* Do not think the next four lines mean anything in this case, just copied from the general case */ @@ -9017,8 +9028,8 @@ uint64_t gmt_map_clip_path(struct GMT_CTRL *GMT, double **x, double **y, bool *d Couldn't do the same to the other two corners because permanently got Inf's. Anyway, for porposes of map clipping, this seems good enough. */ - gmt_geo_to_xy(GMT, -113.0447807067042, 49.56674656682158, &work_x[3], &work_y[3]); /* UL */ - gmt_geo_to_xy(GMT, -113.06804976730443, 49.553963054590234, &work_x[1], &work_y[1]); /* LR */ + gmt_geo_to_xy(GMT, GMT_PROJ_SPILH_LON_UL, GMT_PROJ_SPILH_LAT_UL, &work_x[3], &work_y[3]); /* UL */ + gmt_geo_to_xy(GMT, GMT_PROJ_SPILH_LON_LR, GMT_PROJ_SPILH_LAT_LR, &work_x[1], &work_y[1]); /* LR */ work_x[0] = work_x[3]; work_x[2] = work_x[1]; work_y[0] = work_y[1]; work_y[2] = work_y[3]; break; diff --git a/src/gmt_proj.c b/src/gmt_proj.c index 56bbd0182c9..2947f36e22b 100644 --- a/src/gmt_proj.c +++ b/src/gmt_proj.c @@ -90,6 +90,20 @@ #define GMT_PROJ_CONV_LIMIT 1e-9 #define gmt_M_proj_is_zero(x) (fabs (x) < GMT_PROJ_CONV_LIMIT) +/* The Spilhaus proj is f. The points bellow were obtained by projecting a global grid at 1 + deg resolution with gdalwarp and inverting the UL and LR coordinates. But even those + had to be moved a bit to inside (by trial en error). Couldn't do the same to the other + two corners because permanently got Inf's. This could probably still be refined a tinny bit. + They are good for global maps, but will be in excess for some local -R. However, given + the complexity of the projection, it is not worth the effort to determine them for all + cases. Anyway, the beauty of the Spilhaus projection is when applyied to the whole world, + so it is not a big deal. +*/ +#define GMT_PROJ_SPILH_LON_UL -113.0447807067042 +#define GMT_PROJ_SPILH_LAT_UL 49.56674656682158 +#define GMT_PROJ_SPILH_LON_LR -113.06804976730443 +#define GMT_PROJ_SPILH_LAT_LR 49.553963054590234 + GMT_LOCAL double gmtproj_robinson_spline (struct GMT_CTRL *GMT, double xp, double *x, double *y, double *c) { /* Returns the interpolated value y(xp) from the Robinson coefficients */ diff --git a/src/pscoast.c b/src/pscoast.c index fb0dff977ff..25456ca2b95 100644 --- a/src/pscoast.c +++ b/src/pscoast.c @@ -199,10 +199,10 @@ static int usage (struct GMTAPI_CTRL *API, int level) { if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Usage (API, 0, "usage: %s %s %s [%s] [%s] [-C[+l|r]] [-D[+f]] [-E%s] " "[-G[]] [-F%s] [-I[/]] %s [-L%s] [-M] [-N[/]] %s%s[-Q] [-S[]] " - "[-Td%s] [-Tm%s] [%s] [%s] [-W[/][]] [%s] [%s] [%s] %s[%s] [%s] [%s] [%s]%s [%s]\n", + "[-Td%s] [-Tm%s] [%s] [%s] [-W[/][]] [%s] [%s] [%s] [%s] %s[%s] [%s] [%s] [%s]%s [%s]\n", name, GMT_J_OPT, GMT_Rgeoz_OPT, GMT_A_OPT, GMT_B_OPT, DCW_OPT, GMT_PANEL, API->K_OPT, GMT_SCALE, API->O_OPT, API->P_OPT, GMT_TROSE_DIR, GMT_TROSE_MAG, GMT_U_OPT, GMT_V_OPT, GMT_X_OPT, GMT_Y_OPT, GMT_bo_OPT, API->c_OPT, - GMT_do_OPT, GMT_p_OPT, GMT_t_OPT, GMT_colon_OPT, dbg, GMT_PAR_OPT); + GMT_do_OPT, GMT_g_OPT, GMT_p_OPT, GMT_t_OPT, GMT_colon_OPT, dbg, GMT_PAR_OPT); if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS); @@ -268,6 +268,9 @@ static int usage (struct GMTAPI_CTRL *API, int level) { GMT_Usage (API, 3, "4: Lake in island in lake shores."); GMT_Usage (API, -2, "Note: When feature-specific pens are used, those not set are deactivated."); GMT_Option (API, "X,bo,c,do,p,t"); + GMT_Usage (API, 1, "\n-gD Short version of the global -g option. Here it is used to set the distance, in map " + "units, above which we consider to have a gap. Useful for the Spilhaus projection (though we automatically set it) " + "and when line wrapping on dateline was not correctly detected."); #ifdef DEBUG GMT_Usage (API, 1, "\n-+ (repeatable up to 16 times)"); GMT_Usage (API, -2, "Plot only the specified bins (debug option)."); @@ -537,6 +540,8 @@ static int parse (struct GMT_CTRL *GMT, struct PSCOAST_CTRL *Ctrl, struct GMT_OP if (!GMT->common.g.active && GMT->current.proj.projection == GMT_PROJ4_SPILHAUS) { gmt_parse_g_option(GMT, "D1"); GMT->common.g.active = true; + if (!((GMT->common.R.wesn[XHI] - GMT->common.R.wesn[XLO]) == 360 && GMT->common.R.wesn[YLO] == -90 && GMT->common.R.wesn[YHI] == 90)) + GMT_Report(API, GMT_MSG_WARNING, "Using a non-global region with Spilhaus projection has unknown effects.\n"); } if ((error = gmt_DCW_list (GMT, &(Ctrl->E.info)))) { /* This is either success or failure... */