aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile.am4
-rw-r--r--src/PJ_bertin1953.c94
-rw-r--r--src/PJ_lcc.c33
-rw-r--r--src/PJ_tobmerc.c51
-rw-r--r--src/lib_proj.cmake2
-rw-r--r--src/pj_gridinfo.c4
-rw-r--r--src/pj_list.h2
7 files changed, 168 insertions, 22 deletions
diff --git a/src/Makefile.am b/src/Makefile.am
index ae5e0a2e..359f2f28 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -53,9 +53,9 @@ libproj_la_SOURCES = \
PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c PJ_gall.c PJ_geoc.c \
PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
PJ_mill.c PJ_ocea.c PJ_omerc.c PJ_somerc.c \
- PJ_tcc.c PJ_tcea.c PJ_times.c PJ_tmerc.c \
+ PJ_tcc.c PJ_tcea.c PJ_times.c PJ_tmerc.c PJ_tobmerc.c \
PJ_airy.c PJ_aitoff.c PJ_august.c PJ_bacon.c \
- PJ_chamb.c PJ_hammer.c PJ_lagrng.c PJ_larr.c \
+ PJ_bertin1953.c PJ_chamb.c PJ_hammer.c PJ_lagrng.c PJ_larr.c \
PJ_lask.c PJ_latlong.c PJ_nocol.c PJ_ob_tran.c PJ_oea.c \
PJ_tpeqd.c PJ_vandg.c PJ_vandg2.c PJ_vandg4.c \
PJ_wag7.c PJ_lcca.c PJ_geos.c proj_etmerc.c \
diff --git a/src/PJ_bertin1953.c b/src/PJ_bertin1953.c
new file mode 100644
index 00000000..5a027da2
--- /dev/null
+++ b/src/PJ_bertin1953.c
@@ -0,0 +1,94 @@
+/*
+ Created by Jacques Bertin in 1953, this projection was the go-to choice
+ of the French cartographic school when they wished to represent phenomena
+ on a global scale.
+
+ Formula designed by Philippe Rivière, 2017.
+ https://visionscarto.net/bertin-projection-1953
+
+ Port to PROJ by Philippe Rivière, 21 September 2018
+*/
+
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "projects.h"
+
+PROJ_HEAD(bertin1953, "Bertin 1953")
+ "\n\tMisc Sph no inv.";
+
+struct pj_opaque {
+ double cos_delta_phi, sin_delta_phi, cos_delta_gamma, sin_delta_gamma, deltaLambda;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {
+ XY xy = {0.0,0.0};
+ struct pj_opaque *Q = P->opaque;
+
+ double fu = 1.4, k = 12., w = 1.68, d;
+
+ /* Rotate */
+ double cosphi, x, y, z, z0;
+ lp.lam += PJ_TORAD(-16.5);
+ cosphi = cos(lp.phi);
+ x = cos(lp.lam) * cosphi;
+ y = sin(lp.lam) * cosphi;
+ z = sin(lp.phi);
+ z0 = z * Q->cos_delta_phi + x * Q->sin_delta_phi;
+ lp.lam = atan2(y * Q->cos_delta_gamma - z0 * Q->sin_delta_gamma,
+ x * Q->cos_delta_phi - z * Q->sin_delta_phi);
+ z0 = z0 * Q->cos_delta_gamma + y * Q->sin_delta_gamma;
+ lp.phi = asin(z0);
+
+ lp.lam = adjlon(lp.lam);
+
+ /* Adjust pre-projection */
+ if (lp.lam + lp.phi < -fu) {
+ d = (lp.lam - lp.phi + 1.6) * (lp.lam + lp.phi + fu) / 8.;
+ lp.lam += d;
+ lp.phi -= 0.8 * d * sin(lp.phi + M_PI / 2.);
+ }
+
+ /* Project with Hammer (1.68,2) */
+ cosphi = cos(lp.phi);
+ d = sqrt(2./(1. + cosphi * cos(lp.lam / 2.)));
+ xy.x = w * d * cosphi * sin(lp.lam / 2.);
+ xy.y = d * sin(lp.phi);
+
+ /* Adjust post-projection */
+ d = (1. - cos(lp.lam * lp.phi)) / k;
+ if (xy.y < 0.) {
+ xy.x *= 1. + d;
+ }
+ if (xy.y > 0.) {
+ xy.x *= 1. + d / 1.5 * xy.x * xy.x;
+ }
+
+ return xy;
+}
+
+
+PJ *PROJECTION(bertin1953) {
+ struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ if (0==Q)
+ return pj_default_destructor (P, ENOMEM);
+ P->opaque = Q;
+
+ P->lam0 = 0;
+ P->phi0 = PJ_TORAD(-42.);
+
+ Q->cos_delta_phi = cos(P->phi0);
+ Q->sin_delta_phi = sin(P->phi0);
+ Q->cos_delta_gamma = 1.;
+ Q->sin_delta_gamma = 0.;
+
+ P->es = 0.;
+ P->fwd = s_forward;
+
+ return P;
+}
diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c
index 34ed99cb..96aa3f2a 100644
--- a/src/PJ_lcc.c
+++ b/src/PJ_lcc.c
@@ -5,9 +5,9 @@
#include "proj_math.h"
PROJ_HEAD(lcc, "Lambert Conformal Conic")
- "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0";
+ "\n\tConic, Sph&Ell\n\tlat_1= and lat_2= or lat_0, k_0=";
-# define EPS10 1.e-10
+#define EPS10 1.e-10
struct pj_opaque {
double phi1;
@@ -15,12 +15,11 @@ struct pj_opaque {
double n;
double rho0;
double c;
- int ellips;
};
static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
- XY xy = {0.0,0.0};
+ XY xy = {0., 0.};
struct pj_opaque *Q = P->opaque;
double rho;
@@ -31,18 +30,19 @@ static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
}
rho = 0.;
} else {
- rho = Q->c * (Q->ellips ? pow(pj_tsfn(lp.phi, sin(lp.phi),
- P->e), Q->n) : pow(tan(M_FORTPI + .5 * lp.phi), -Q->n));
+ rho = Q->c * (P->es != 0. ?
+ pow(pj_tsfn(lp.phi, sin(lp.phi), P->e), Q->n) :
+ pow(tan(M_FORTPI + .5 * lp.phi), -Q->n));
}
lp.lam *= Q->n;
- xy.x = P->k0 * (rho * sin( lp.lam) );
- xy.y = P->k0 * (Q->rho0 - rho * cos(lp.lam) );
+ xy.x = P->k0 * (rho * sin(lp.lam));
+ xy.y = P->k0 * (Q->rho0 - rho * cos(lp.lam));
return xy;
}
static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
- LP lp = {0.0,0.0};
+ LP lp = {0., 0.};
struct pj_opaque *Q = P->opaque;
double rho;
@@ -51,13 +51,13 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
xy.y = Q->rho0 - xy.y;
rho = hypot(xy.x, xy.y);
- if (rho != 0.0) {
+ if (rho != 0.) {
if (Q->n < 0.) {
rho = -rho;
xy.x = -xy.x;
xy.y = -xy.y;
}
- if (Q->ellips) {
+ if (P->es != 0.) {
lp.phi = pj_phi2(P->ctx, pow(rho / Q->c, 1./Q->n), P->e);
if (lp.phi == HUGE_VAL) {
proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
@@ -78,13 +78,12 @@ static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
PJ *PROJECTION(lcc) {
double cosphi, sinphi;
int secant;
- struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+ struct pj_opaque *Q = pj_calloc(1, sizeof (struct pj_opaque));
- if (0==Q)
- return pj_default_destructor (P, ENOMEM);
+ if (0 == Q)
+ return pj_default_destructor(P, ENOMEM);
P->opaque = Q;
-
Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
if (pj_param(P->ctx, P->params, "tlat_2").i)
Q->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
@@ -99,10 +98,9 @@ PJ *PROJECTION(lcc) {
Q->n = sinphi = sin(Q->phi1);
cosphi = cos(Q->phi1);
secant = fabs(Q->phi1 - Q->phi2) >= EPS10;
- if( (Q->ellips = (P->es != 0.)) ) {
+ if (P->es != 0.) {
double ml1, m1;
- P->e = sqrt(P->es);
m1 = pj_msfn(sinphi, cosphi, P->es);
ml1 = pj_tsfn(Q->phi1, sinphi, P->e);
if (secant) { /* secant cone */
@@ -128,4 +126,3 @@ PJ *PROJECTION(lcc) {
return P;
}
-
diff --git a/src/PJ_tobmerc.c b/src/PJ_tobmerc.c
new file mode 100644
index 00000000..9c939f0b
--- /dev/null
+++ b/src/PJ_tobmerc.c
@@ -0,0 +1,51 @@
+#define PJ_LIB__
+
+#include <float.h>
+#include <math.h>
+
+#include "proj_internal.h"
+#include "proj.h"
+#include "proj_math.h"
+#include "projects.h"
+
+PROJ_HEAD(tobmerc, "Tobler-Mercator") "\n\tCyl, Sph";
+
+#define EPS10 1.e-10
+static double logtanpfpim1(double x) { /* log(tan(x/2 + M_FORTPI)) */
+ if (fabs(x) <= DBL_EPSILON) {
+ /* tan(M_FORTPI + .5 * x) can be approximated by 1.0 + x */
+ return log1p(x);
+ }
+ return log(tan(M_FORTPI + .5 * x));
+}
+
+static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
+ XY xy = {0.0, 0.0};
+ double cosphi;
+
+ if (fabs(fabs(lp.phi) - M_HALFPI) <= EPS10) {
+ proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+ return xy;
+ }
+
+ cosphi = cos(lp.phi);
+ xy.x = P->k0 * lp.lam * cosphi * cosphi;
+ xy.y = P->k0 * logtanpfpim1(lp.phi);
+ return xy;
+}
+
+static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
+ LP lp = {0.0, 0.0};
+ double cosphi;
+
+ lp.phi = atan(sinh(xy.y / P->k0));
+ cosphi = cos(lp.phi);
+ lp.lam = xy.x / P->k0 / (cosphi * cosphi);
+ return lp;
+}
+
+PJ *PROJECTION(tobmerc) {
+ P->inv = s_inverse;
+ P->fwd = s_forward;
+ return P;
+}
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index a359369e..6a287a43 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -46,6 +46,7 @@ SET(SRC_LIBPROJ_PJ
PJ_august.c
PJ_axisswap.c
PJ_bacon.c
+ PJ_bertin1953.c
PJ_bipc.c
PJ_boggs.c
PJ_bonne.c
@@ -140,6 +141,7 @@ SET(SRC_LIBPROJ_PJ
PJ_tcea.c
PJ_times.c
PJ_tmerc.c
+ PJ_tobmerc.c
PJ_tpeqd.c
PJ_unitconvert.c
PJ_urm5.c
diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
index f201f39e..de0e8d31 100644
--- a/src/pj_gridinfo.c
+++ b/src/pj_gridinfo.c
@@ -660,7 +660,7 @@ static int pj_gridinfo_init_ntv2( projCtx ctx, PAFile fid, PJ_GRIDINFO *gilist )
static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
{
- unsigned char header[176];
+ unsigned char header[192]; /* 12 records of 16 bytes */
struct CTABLE *ct;
LP ur;
@@ -731,7 +731,7 @@ static int pj_gridinfo_init_ntv1( projCtx ctx, PAFile fid, PJ_GRIDINFO *gi )
ct->cvs = NULL;
gi->ct = ct;
- gi->grid_offset = pj_ctx_ftell( ctx, fid );
+ gi->grid_offset = (long) sizeof(header);
gi->format = "ntv1";
return 1;
diff --git a/src/pj_list.h b/src/pj_list.h
index 72a9af7f..8c72dd1f 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -15,6 +15,7 @@ PROJ_HEAD(apian, "Apian Globular I")
PROJ_HEAD(august, "August Epicycloidal")
PROJ_HEAD(axisswap, "Axis ordering")
PROJ_HEAD(bacon, "Bacon Globular")
+PROJ_HEAD(bertin1953, "Bertin 1953")
PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
PROJ_HEAD(boggs, "Boggs Eumorphic")
PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")
@@ -139,6 +140,7 @@ PROJ_HEAD(tcea, "Transverse Cylindrical Equal Area")
PROJ_HEAD(times, "Times Projection")
PROJ_HEAD(tissot, "Tissot Conic")
PROJ_HEAD(tmerc, "Transverse Mercator")
+PROJ_HEAD(tobmerc, "Tobler-Mercator")
PROJ_HEAD(tpeqd, "Two Point Equidistant")
PROJ_HEAD(tpers, "Tilted perspective")
PROJ_HEAD(unitconvert, "Unit conversion")
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-15 02:19:47 +0000