Skip to content

Commit

Permalink
Spilhaus: better use lon_0 and lat_0
Browse files Browse the repository at this point in the history
  • Loading branch information
jjimenezshaw committed Feb 18, 2025
1 parent 05f6183 commit 85e7d8e
Show file tree
Hide file tree
Showing 3 changed files with 41 additions and 12 deletions.
4 changes: 2 additions & 2 deletions docs/source/operations/projections/spilhaus.rst
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ Parameters

.. note:: All parameters are optional.

.. option:: +lon_c=<value>
.. option:: +lon_0=<value>

Longitude of projection centre.

Expand All @@ -49,7 +49,7 @@ Parameters

See :ref:`Projection Units <projection_units>` for more information.

.. option:: +lat_c=<value>
.. option:: +lat_0=<value>

Latitude of projection centre.

Expand Down
19 changes: 12 additions & 7 deletions src/projections/spilhaus.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -135,25 +135,30 @@ PJ *PJ_PROJECTION(spilhaus) {
? pj_param(P->ctx, P->params, ("r" + name).c_str()).f
: def * DEG_TO_RAD;
};

// Values from https://github.com/OSGeo/PROJ/issues/1851
const double lat_center = param_rad("lat_c", -49.56371678);
const double lon_center = param_rad("lon_c", 66.94970198);
if (!pj_param(P->ctx, P->params, "tlon_0").i) {
P->lam0 = 66.94970198 * DEG_TO_RAD;
}
if (!pj_param(P->ctx, P->params, "tlat_0").i) {
P->phi0 = -49.56371678 * DEG_TO_RAD;
}
const double azimuth = param_rad("azi", 40.17823482);

const double rotation = param_rad("rot", 45);
Q->cosrot = cos(rotation);
Q->sinrot = sin(rotation);

const double conformal_lat_center = pj_conformal_lat(lat_center, P->e);
const double conformal_lat_center = pj_conformal_lat(P->phi0, P->e);

Q->sinalpha = -cos(conformal_lat_center) * cos(azimuth);
Q->cosalpha = sqrt(1 - Q->sinalpha * Q->sinalpha);
Q->lambda_0 = lon_center + atan2(tan(azimuth), -sin(conformal_lat_center));
Q->lambda_0 = atan2(tan(azimuth), -sin(conformal_lat_center));
Q->beta = M_PI + atan2(-sin(azimuth), -tan(conformal_lat_center));

Q->conformal_distortion =
cos(lat_center) / sqrt(1 - P->es * sin(lat_center) * sin(lat_center)) /
cos(conformal_lat_center);
Q->conformal_distortion = cos(P->phi0) /
sqrt(1 - P->es * sin(P->phi0) * sin(P->phi0)) /
cos(conformal_lat_center);

P->fwd = spilhaus_forward;
P->inv = spilhaus_inverse;
Expand Down
30 changes: 27 additions & 3 deletions test/gie/spilhaus.gie
Original file line number Diff line number Diff line change
Expand Up @@ -367,11 +367,12 @@ expect -15390625 -15546875
------------------------------------------------------------
# Stable for default parameters
------------------------------------------------------------
operation +proj=spilhaus +rot=45 +k_0=1 +lat_c=-49.56371678 +lon_c=66.94970198 +azi=40.17823482
operation +proj=spilhaus +rot=45 +k_0=1 +lat_0=-49.56371678 +lon_0=66.94970198 +azi=40.17823482
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 3733410.0118 -9320.8573
roundtrip 1

------------------------------------------------------------
# Sentitive to input parameters
Expand All @@ -381,36 +382,42 @@ tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 3733410.0118 -9320.8573
roundtrip 1
------------------------------------------------------------
operation +proj=spilhaus +lon_c=10.1
operation +proj=spilhaus +lon_0=10.1
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 4343770.7991 -3701935.6242
roundtrip 1
------------------------------------------------------------
operation +proj=spilhaus +lat_c=30.1
operation +proj=spilhaus +lat_0=30.1
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 3637341.2895 -2571368.8666
roundtrip 1
------------------------------------------------------------
operation +proj=spilhaus +azi=9.1
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 3061806.4542 -1678791.7428
roundtrip 1
------------------------------------------------------------
operation +proj=spilhaus +rot=40.1
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 3720561.6630 309609.603620
roundtrip 1
------------------------------------------------------------
operation +proj=spilhaus +k_0=0.9
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 3360069.0106 -8388.7716
roundtrip 1

------------------------------------------------------------
# Sphere
Expand All @@ -422,4 +429,21 @@ accept 130.4 -16.2
expect 3737644.5177 -7049.7883
roundtrip 1

------------------------------------------------------------
# vs Adams WS2
------------------------------------------------------------
operation +proj=adams_ws2 +R=6378137
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 8199312.0391 -1392652.9172
roundtrip 1
------------------------------------------------------------
operation +proj=spilhaus +R=6378137 +lon_0=0 +lat_0=0 +azi=0 +rot=0
tolerance 1 mm
------------------------------------------------------------
accept 130.4 -16.2
expect 8199312.0391 -1392652.9172
roundtrip 1

</gie-strict>

0 comments on commit 85e7d8e

Please sign in to comment.