Skip to content
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.

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

added new projection function for Navionics Mercator #391

Closed
wants to merge 7 commits into from
2 changes: 1 addition & 1 deletion src/Makefile.am
Expand Up @@ -39,7 +39,7 @@ libproj_la_SOURCES = \
PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
PJ_rpoly.c PJ_sconics.c proj_rouss.c \
PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c \
PJ_gall.c PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
PJ_gall.c PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c PJ_nav_merc.c \
PJ_mill.c PJ_ocea.c PJ_omerc.c PJ_somerc.c \
PJ_tcc.c PJ_tcea.c PJ_tmerc.c \
PJ_airy.c PJ_aitoff.c PJ_august.c PJ_bacon.c \
Expand Down
132 changes: 132 additions & 0 deletions src/PJ_nav_merc.c
@@ -0,0 +1,132 @@
/*
***** Navionics Mercator ****************
*
* (co)author pasquini.matteo@gmail.com, mpasquini@navionics.com
*
* datum(barely speaking):
* based on Hayford 1909 ellipsoid but geocentered.
*
* attributes:
* - geocentered
* - a = 6378388
* - b = 6356912
* - (a/b)/a = 1/297.0
* - eccentricity = 0.081991889977
* - KC = 1.00676429271698
*
* orignal functions are in decimal degrees(instead radians)
* this implementation is(should be) the most precise in radians,
* the first outside proprietary Navionics S.p.A. SDK
*
* note the KC is a CONST value that substitute the P->e
* it have some similarity with PJ_merc and PJ_tsfn
* , all maths help are wellcome ...
*
* This few lines represent a "small docs"
* complete documentation is in .docx format, to not pollute this beautiful project
* i prefer that you bother me throught email request.
*
* this is my first piece of C in the real world.
*
*/



#define PJ_LIB__
#include <projects.h>

PROJ_HEAD(nav_merc, "Navionics Mercator") "\n\tEll\n\tlat_ts=";

#define EPS10 1.e-10
#define KC 1.00676429271698


static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
XY xy = {0.0,0.0};
if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10)
F_ERROR;

xy.y = log(tan( .5 * atan (tan(lp.phi)/KC) + M_FORTPI) ) ;
xy.x = P->k0 * lp.lam;
return (xy);
}

static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
LP lp = {0.0,0.0};
if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL)
I_ERROR;

lp.phi = atan2((tan((atan(exp(xy.y)) * 2.) - M_HALFPI) * KC), 1);
lp.lam = xy.x / P->k0;
return lp;
}

static void freeup(PJ *P) { /* Destructor */
pj_dealloc(P);
}


PJ *PROJECTION(nav_merc) {
double phits=0.0;
int is_phits;

if( (is_phits = pj_param(P->ctx, P->params, "tlat_ts").i) ) {
phits = fabs(pj_param(P->ctx, P->params, "rlat_ts").f);
if (phits >= M_HALFPI) E_ERROR(-24);
}

if (P->es) { /* ellipsoid */
if (is_phits)
P->k0 = pj_msfn(sin(phits), cos(phits), P->es);
P->inv = e_inverse;
P->fwd = e_forward;
}

return P;
}


#ifdef PJ_OMIT_SELFTEST
int pj_nav_merc_selftest (void) {return 0;}
#else

int pj_nav_merc_selftest (void) {
double tolerance_lp = 1e-10;
double tolerance_xy = 1e-7;

char e_args[] = {" -w6 +proj=nav_merc +a=6378388.0 +rf=297"};
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

+a=6378388.0 +rf=297 could be replaced with the shorthand +ellps=intl. This is a slight preference.

Here is how I looked that up from the command line.

$ proj -le
[..snip..]
  helmert a=6378200.       rf=298.3         Helmert 1906
    hough a=6378270.0      rf=297.          Hough
     intl a=6378388.0      rf=297.          International 1909 (Hayford)
    krass a=6378245.0      rf=298.3         Krassovsky, 1942
    kaula a=6378163.       rf=298.24        Kaula 1961
[..snip..]


LP fwd_in[] = {
{ 950594.539, 5968306.230 },
{ 1178796.03, 5276463.78 },
{ 18671158.340, -5945578.371 },
{ 20038296.88, -18765191.35 }
};

XY e_fwd_expect[] = {
{ 8.5390000000, 47.3515666667 },
{ 10.5888881683, 42.9574966432 },
{ 167.7192687988, -47.2126350403 },
{ 179.9998888889, -84.0 },
};

XY inv_in[] = {
{ 8.5390000000, 47.3515666667 },
{ 10.5888881683, 42.9574966432 },
{ 167.7192687988, -47.2126350403 },
{ 180.0, -84.0 }
};

LP e_inv_expect[] = {
{ 950594.539, 5968306.230 },
{ 1178796.03, 5276463.78 },
{ 18671158.340, -5945578.371 },
{ 20038296.883, -18765191.347 },
};


return pj_generic_selftest (e_args, 0, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, e_fwd_expect, inv_in, 0, e_inv_expect );
}


#endif
1 change: 1 addition & 0 deletions src/lib_proj.cmake
Expand Up @@ -99,6 +99,7 @@ SET(SRC_LIBPROJ_PJ
PJ_moll.c
PJ_natearth.c
PJ_natearth2.c
PJ_nav_merc.c
PJ_nell.c
PJ_nell_h.c
PJ_nocol.c
Expand Down
1 change: 1 addition & 0 deletions src/pj_list.h
Expand Up @@ -78,6 +78,7 @@ PROJ_HEAD(mbtfpp, "McBride-Thomas Flat-Polar Parabolic")
PROJ_HEAD(mbtfpq, "McBryde-Thomas Flat-Polar Quartic")
PROJ_HEAD(mbtfps, "McBryde-Thomas Flat-Polar Sinusoidal")
PROJ_HEAD(merc, "Mercator")
PROJ_HEAD(nav_merc, "Navionics Mercator")
PROJ_HEAD(mil_os, "Miller Oblated Stereographic")
PROJ_HEAD(mill, "Miller Cylindrical")
PROJ_HEAD(misrsom, "Space oblique for MISR")
Expand Down