Importance of -R for mapproject

Q&A
1

Hi

(GMT version 6.3.0)

I’m trying to understand the importance of the region specified on the -R option when using mapproject to convert between geographic coordinates and Cartesian. For the same input coordinates and projection it seems that the region affects how the transformation is done, resulting in differences in the result. e.g.

$~>echo "-137000.0 11000.0" | gmt mapproject -Jt-17.927/64.463/1:1 -R-23/-7/60/66 -I -C -F -V
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Using spherical projection with authalic latitudes
mapproject [INFORMATION]: Transform -2.30000000e+01/-7.00000000e+00/6.00000000e+01/6.60000000e+01 <- -2.83125208e+05/6.08385783e+05/-5.14521453e+05/1.98117729e+05 [m]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
-2.07948097e+01	6.46830641e+01
mapproject [INFORMATION]: Projected 1 points
mapproject [INFORMATION]: Input extreme values:  Xmin: -1.37000000e+05 Xmax: -1.37000000e+05 Ymin: 1.10000000e+04 Ymax 1.10000000e+04
mapproject [INFORMATION]: Output extreme values: Xmin: -2.07948097e+01 Xmax: -2.07948097e+01 Ymin: 6.46830641e+01 Ymax 6.46830641e+01
mapproject [INFORMATION]: Mapped 1 x-y pairs [m] to lon-lat pairs

Yields a different results than

$~> echo "-137000.0 11000.0" | gmt mapproject -Jt-17.927/64.463/1:1 -R-23/-8/60/66 -I -C -F -V
mapproject [INFORMATION]: Processing input table data
mapproject [INFORMATION]: Transform -2.30000000e+01/-8.00000000e+00/6.00000000e+01/6.60000000e+01 <- -2.82888271e+05/5.52531638e+05/-4.97396944e+05/2.07039344e+05 [m]
mapproject [INFORMATION]: Reading Data Table from Standard Input stream
mapproject [INFORMATION]: Writing Data Table to Standard Output stream
-2.07821645e+01	6.45340166e+01
mapproject [INFORMATION]: Projected 1 points
mapproject [INFORMATION]: Input extreme values:  Xmin: -1.37000000e+05 Xmax: -1.37000000e+05 Ymin: 1.10000000e+04 Ymax 1.10000000e+04
mapproject [INFORMATION]: Output extreme values: Xmin: -2.07821645e+01 Xmax: -2.07821645e+01 Ymin: 6.45340166e+01 Ymax 6.45340166e+01
mapproject [INFORMATION]: Mapped 1 x-y pairs [m] to lon-lat pairs

The only difference being that the region in the latter call is 1 degree wider in the longitude. I note that the smaller region triggers spherical projection with authalic latitudes so presumably this is the cause of the difference (about 16.64 km) here but I assume that the projection should be unique so asking for some help on how to know which region is appropriate to use.

2

The ellipsoidal series expansion of Snyder breaks down if the longitude range exceeds ~10 degrees. If -R indicates the larger range then we switch to the spherical solution but use authalic latitudes. So
knowing your area range helps set the -R and get the most accurate solution (which is the ellipsoidal expansion as long as you are not too far from the central longitude.

3

Thanks, that makes sense althouhg I’m a bit surprised that the difference is so large (16.6 km difference at ~140 km from the center of projection).

Is this documented in the GMT docs?

4

Not really. 16 km must be a sign of a bug. I will look at your example tomorrow and see what is going on here.

5

Verified and updated master so we are able to match Snyder’s published examples. Authalic was a typo, we do conformal latitudes if spherical TM. Given you are just < 3 degrees from the central meridian you should definitively use the ellipsoidal terms but since you are doing the inverse transformation you need to know which ellipsoid was used and what map scale factor was used (0.9996 ?), then make sure those are set accordingly ( PROJ_ELLIPSOID, PROJ_SCALE_FACTOR)

6

Great, thanks.

I’m basically converting back and forth of the same coordinate set so have been assuming the the default ellipsoid and scale factor have not changed in between, but perhaps safer to set these explicitly just in case.