-
Notifications
You must be signed in to change notification settings - Fork 5
/
PJ_calcofi.c
138 lines (135 loc) · 5.1 KB
/
PJ_calcofi.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#define PJ_LIB__
#include <projects.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <proj_api.h>
#include <errno.h>
/* Conversions for the California Cooperative Oceanic Fisheries Investigations
Line/Station coordinate system following the algorithm of:
Eber, L.E., and R.P. Hewitt. 1979. Conversion algorithms for the CALCOFI
station grid. California Cooperative Oceanic Fisheries Investigations Reports
20:135-137. (corrected for typographical errors).
http://www.calcofi.org/publications/calcofireports/v20/Vol_20_Eber___Hewitt.pdf
They assume 1 unit of CalCOFI Line == 1/5 degree in longitude or
meridional units at reference point O, and similarly 1 unit of CalCOFI
Station == 1/15 of a degree at O.
By convention, CalCOFI Line/Station conversions use Clarke 1866 but we use
whatever ellipsoid is provided. */
PROJ_HEAD(calcofi,
"Cal Coop Ocean Fish Invest Lines/Stations") "\n\tCyl, Sph&Ell";
#define EPS10 1.e-10
#define DEG_TO_LINE 5
#define DEG_TO_STATION 15
#define LINE_TO_RAD 0.0034906585039886592
#define STATION_TO_RAD 0.0011635528346628863
#define PT_O_LINE 80 /* reference point O is at line 80, */
#define PT_O_STATION 60 /* station 60, */
#define PT_O_LAMBDA -2.1144663887911301 /* lon -121.15 and */
#define PT_O_PHI 0.59602993955606354 /* lat 34.15 */
#define ROTATION_ANGLE 0.52359877559829882 /*CalCOFI angle of 30 deg in rad */
FORWARD(e_forward); /* ellipsoid */
double oy; /* pt O y value in Mercator */
double l1; /* l1 and l2 are distances calculated using trig that sum
to the east/west distance between point O and point xy */
double l2;
double ry; /* r is the point on the same station as o (60) and the same
line as xy xy, r, o form a right triangle */
/* if the user has specified +lon_0 or +k0 for some reason,
we're going to ignore it so that xy is consistent with point O */
lp.lam = lp.lam + P->lam0;
if (fabs(fabs(lp.phi) - HALFPI) <= EPS10) F_ERROR;
xy.x = lp.lam;
xy.y = -log(pj_tsfn(lp.phi, sin(lp.phi), P->e)); /* Mercator transform xy*/
oy = -log(pj_tsfn(PT_O_PHI, sin(PT_O_PHI), P->e));
l1 = (xy.y - oy) * tan(ROTATION_ANGLE);
l2 = -xy.x - l1 + PT_O_LAMBDA;
ry = l2 * cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE) + xy.y;
ry = pj_phi2(P->ctx, exp(-ry), P->e); /*inverse Mercator*/
xy.x = PT_O_LINE - RAD_TO_DEG *
(ry - PT_O_PHI) * DEG_TO_LINE / cos(ROTATION_ANGLE);
xy.y = PT_O_STATION + RAD_TO_DEG *
(ry - lp.phi) * DEG_TO_STATION / sin(ROTATION_ANGLE);
/* set a = 1, x0 = 0, and y0 = 0 so that no further unit adjustments
are done */
P->a = 1;
P->x0 = 0;
P->y0 = 0;
return (xy);
}
FORWARD(s_forward); /* spheroid */
double oy;
double l1;
double l2;
double ry;
lp.lam = lp.lam + P->lam0;
if (fabs(fabs(lp.phi) - HALFPI) <= EPS10) F_ERROR;
xy.x = lp.lam;
xy.y = log(tan(FORTPI + .5 * lp.phi));
oy = log(tan(FORTPI + .5 * PT_O_PHI));
l1 = (xy.y - oy) * tan(ROTATION_ANGLE);
l2 = -xy.x - l1 + PT_O_LAMBDA;
ry = l2 * cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE) + xy.y;
ry = HALFPI - 2. * atan(exp(-ry));
xy.x = PT_O_LINE - RAD_TO_DEG *
(ry - PT_O_PHI) * DEG_TO_LINE / cos(ROTATION_ANGLE);
xy.y = PT_O_STATION + RAD_TO_DEG *
(ry - lp.phi) * DEG_TO_STATION / sin(ROTATION_ANGLE);
P->a = 1;
P->x0 = 0;
P->y0 = 0;
return (xy);
}
INVERSE(e_inverse); /* ellipsoid */
double ry; /* y value of point r */
double oymctr; /* Mercator-transformed y value of point O */
double rymctr; /* Mercator-transformed ry */
double xymctr; /* Mercator-transformed xy.y */
double l1;
double l2;
/* turn x and y back into Line/Station */
xy.x /= P->ra;
xy.y /= P->ra;
ry = PT_O_PHI - LINE_TO_RAD * (xy.x - PT_O_LINE) *
cos(ROTATION_ANGLE);
lp.phi = ry - STATION_TO_RAD * (xy.y - PT_O_STATION) * sin(ROTATION_ANGLE);
oymctr = -log(pj_tsfn(PT_O_PHI, sin(PT_O_PHI), P->e));
rymctr = -log(pj_tsfn(ry, sin(ry), P->e));
xymctr = -log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
l1 = (xymctr - oymctr) * tan(ROTATION_ANGLE);
l2 = (rymctr - xymctr) / (cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE));
lp.lam = PT_O_LAMBDA - (l1 + l2);
P->over = 1;
return (lp);
}
INVERSE(s_inverse); /* spheroid */
double ry;
double oymctr;
double rymctr;
double xymctr;
double l1;
double l2;
xy.x /= P->ra;
xy.y /= P->ra;
ry = PT_O_PHI - LINE_TO_RAD * (xy.x - PT_O_LINE) *
cos(ROTATION_ANGLE);
lp.phi = ry - STATION_TO_RAD * (xy.y - PT_O_STATION) * sin(ROTATION_ANGLE);
oymctr = log(tan(FORTPI + .5 * PT_O_PHI));
rymctr = log(tan(FORTPI + .5 * ry));
xymctr = log(tan(FORTPI + .5 * lp.phi));
l1 = (xymctr - oymctr) * tan(ROTATION_ANGLE);
l2 = (rymctr - xymctr) / (cos(ROTATION_ANGLE) * sin(ROTATION_ANGLE));
lp.lam = PT_O_LAMBDA - (l1 + l2);
P->over = 1;
return (lp);
}
FREEUP; if (P) pj_dalloc(P); }
ENTRY0(calcofi)
if (P->es) { /* ellipsoid */
P->inv = e_inverse;
P->fwd = e_forward;
} else { /* sphere */
P->inv = s_inverse;
P->fwd = s_forward;
}
ENDENTRY(P)