Skip to content

Commit

Permalink
Fix bad effect of roundoff at boarder (#8238)
Browse files Browse the repository at this point in the history
* Fix bad effect of roundoff at boarder

See this foum post for background.  Two issues:

Mollweide should be run on sphere.
Need to allow 10-9 slop at ±180 longitude.

WIth this we get

echo 180 0 | gmt mapproject -Jw0/1:1 -Fe -C --PROJ_ELLIPSOID=Sphere
18019934.021	0
echo 18019934.021 0 | gmt mapproject -Jw0/1:1 -Fe -C --PROJ_ELLIPSOID=Sphere -I
180	0

* Update gmt_proj.c
  • Loading branch information
PaulWessel authored Dec 26, 2023
1 parent d183cd5 commit d8c0d9f
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 2 deletions.
1 change: 1 addition & 0 deletions src/gmt_constants.h
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@

#define GMT_CONV15_LIMIT 1.0e-15 /* Very tight convergence limit or "close to zero" limit */
#define GMT_CONV12_LIMIT 1.0e-12 /* Tight limit for gaps/overlaps in CPT z-values */
#define GMT_CONV9_LIMIT 1.0e-9 /* Fairly tight convergence limit or "close to zero" limit */
#define GMT_CONV8_LIMIT 1.0e-8 /* Fairly tight convergence limit or "close to zero" limit */
#define GMT_CONV6_LIMIT 1.0e-6 /* 1 ppm */
#define GMT_CONV5_LIMIT 1.0e-5 /* 10 ppm */
Expand Down
6 changes: 4 additions & 2 deletions src/gmt_proj.c
Original file line number Diff line number Diff line change
Expand Up @@ -2054,7 +2054,7 @@ GMT_LOCAL void gmtproj_vmollweide (struct GMT_CTRL *GMT, double lon0, double sca

gmtproj_check_R_J (GMT, &lon0);
GMT->current.proj.central_meridian = lon0;
GMT->current.proj.w_x = GMT->current.proj.EQ_RAD * D2R * d_sqrt (8.0) / M_PI;
GMT->current.proj.w_x = GMT->current.proj.EQ_RAD * D2R * 2.0 * M_SQRT2 / M_PI;
GMT->current.proj.w_y = GMT->current.proj.EQ_RAD * M_SQRT2;
GMT->current.proj.w_iy = 1.0 / GMT->current.proj.w_y;
GMT->current.proj.w_r = 0.25 * (scale * GMT->current.proj.M_PR_DEG * 360.0); /* = Half the minor axis */
Expand Down Expand Up @@ -2097,7 +2097,9 @@ GMT_LOCAL void gmtproj_imollweide (struct GMT_CTRL *GMT, double *lon, double *la

phi = asin (y * GMT->current.proj.w_iy);
*lon = x / (GMT->current.proj.w_x * cos(phi));
if (fabs (*lon) > 180.0) { /* Horizon */
if ((fabs (*lon) - 180.0) < GMT_CONV9_LIMIT)
*lon = copysign (180.0, *lon);
if (fabs (*lon) > 180.0) { /* Beyond Horizon */
*lat = *lon = GMT->session.d_NaN;
return;
}
Expand Down

0 comments on commit d8c0d9f

Please sign in to comment.