add calcofi projection (#135)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@2376 4e78687f-474d-0410-85f9-8d5e500ac6b2

9 files changed, 183 insertions, 9 deletions

diff --git a/ChangeLog b/ChangeLog
index d8d019e0..231490e3 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2013-07-08  Frank Warmerdam  <warmerdam@pobox.com>
+
+	* src/PJ_calcofi.c: Add Cal Coop Ocean Fish Invest Lines/Stations
+	projections (calcofi) (#135)
+
 2013-07-02  Frank Warmerdam  <warmerdam@pobox.com>
 
 	* nad/testvarious, nad/tv_out.dist: add new robinson forward test,
diff --git a/NEWS b/NEWS
index 51242b90..a13f556a 100644
--- a/NEWS
+++ b/NEWS
@@ -12,6 +12,8 @@
  o Removed old (deprecated) Java bindings in favor of the new api introduced
    in 4.8.0.
 
+ o Implement the calcofi (Cal Coop Ocean Fish Invest Lines/Stations) projection
+
  o Various bug fixes and cleanup.
 
 
diff --git a/nad/testvarious b/nad/testvarious
index c1e7f0a2..63ac8423 100755
--- a/nad/testvarious
+++ b/nad/testvarious
@@ -432,6 +432,23 @@ $EXE -f '%.14f' \
  -E >>${OUT} <<EOF
 -6086629.0 4488761.0
 EOF
+echo "##############################################################" >> ${OUT}
+echo "Test forward calcofi projection" >> ${OUT}
+$EXE +proj=latlong +ellps=clrk66 \
+ +to +proj=calcofi +ellps=clrk66 \
+ -E >>${OUT} <<EOF
+120d40'42.273"W	38d56'50.766"N
+121d9'W	34d9'N
+123d59'56.066"W	30d25'4.617"N
+EOF
+echo "Test inverse calcofi projection" >> ${OUT}
+$EXE +proj=calcofi +ellps=clrk66 \
+ +to +proj=longlat +ellps=clrk66 \
+ -E >>${OUT} <<EOF
+60 20
+80 60 
+90 120
+EOF
 ##############################################################################
 # Done! 
 # do 'diff' with distribution results
diff --git a/nad/tv_out.dist b/nad/tv_out.dist
index 8027fca5..7d5a7f28 100644
--- a/nad/tv_out.dist
+++ b/nad/tv_out.dist
@@ -193,3 +193,12 @@ Test pconic (#148)
 ##############################################################
 Test laea
 -6086629.0 4488761.0	156.05863798859943	37.76545829867844 0.00000000000000
+##############################################################
+Test forward calcofi projection
+120d40'42.273"W	38d56'50.766"N	60.00	20.00 0.00
+121d9'W	34d9'N	80.00	60.00 0.00
+123d59'56.066"W	30d25'4.617"N	90.00	120.00 0.00
+Test inverse calcofi projection
+60 20	120d40'42.273"W	38d56'50.766"N 0.000
+80 60	121d9'W	34d9'N 0.000 
+90 120	123d59'56.066"W	30d25'4.617"N 0.000
diff --git a/src/Makefile.am b/src/Makefile.am
index 834bf42a..886075c5 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -57,7 +57,7 @@ libproj_la_SOURCES = \
 	pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \
 	pj_qsfn.c pj_strerrno.c pj_tsfn.c pj_units.c pj_ctx.c pj_log.c \
 	pj_zpoly1.c rtodms.c vector1.c pj_release.c pj_gauss.c \
-	PJ_healpix.c PJ_natearth.c pj_fileapi.c \
+	PJ_healpix.c PJ_natearth.c PJ_calcofi.c pj_fileapi.c \
 	\
 	pj_gc_reader.c pj_gridcatalog.c \
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
diff --git a/src/Makefile.in b/src/Makefile.in
index 019a89d0..1aebcbcf 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -102,12 +102,12 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
 	pj_open_lib.lo pj_param.lo pj_phi2.lo pj_pr_list.lo pj_qsfn.lo \
 	pj_strerrno.lo pj_tsfn.lo pj_units.lo pj_ctx.lo pj_log.lo \
 	pj_zpoly1.lo rtodms.lo vector1.lo pj_release.lo pj_gauss.lo \
-	PJ_healpix.lo PJ_natearth.lo pj_fileapi.lo pj_gc_reader.lo \
-	pj_gridcatalog.lo nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \
-	pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \
-	pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \
-	pj_gridlist.lo jniproj.lo pj_mutex.lo pj_initcache.lo \
-	pj_apply_vgridshift.lo geodesic.lo
+	PJ_healpix.lo PJ_natearth.lo PJ_calcofi.lo pj_fileapi.lo \
+	pj_gc_reader.lo pj_gridcatalog.lo nad_cvt.lo nad_init.lo \
+	nad_intr.lo emess.lo pj_apply_gridshift.lo pj_datums.lo \
+	pj_datum_set.lo pj_transform.lo geocent.lo pj_utils.lo \
+	pj_gridinfo.lo pj_gridlist.lo jniproj.lo pj_mutex.lo \
+	pj_initcache.lo pj_apply_vgridshift.lo geodesic.lo
 libproj_la_OBJECTS = $(am_libproj_la_OBJECTS)
 libproj_la_LINK = $(LIBTOOL) --tag=CC $(AM_LIBTOOLFLAGS) \
 	$(LIBTOOLFLAGS) --mode=link $(CCLD) $(AM_CFLAGS) $(CFLAGS) \
@@ -317,7 +317,7 @@ libproj_la_SOURCES = \
 	pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \
 	pj_qsfn.c pj_strerrno.c pj_tsfn.c pj_units.c pj_ctx.c pj_log.c \
 	pj_zpoly1.c rtodms.c vector1.c pj_release.c pj_gauss.c \
-	PJ_healpix.c PJ_natearth.c pj_fileapi.c \
+	PJ_healpix.c PJ_natearth.c PJ_calcofi.c pj_fileapi.c \
 	\
 	pj_gc_reader.c pj_gridcatalog.c \
 	nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
@@ -484,6 +484,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_bipc.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_boggs.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_bonne.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_calcofi.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_cass.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_cc.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_cea.Plo@am__quote@
diff --git a/src/PJ_calcofi.c b/src/PJ_calcofi.c
new file mode 100644
index 00000000..2d142938
--- /dev/null
+++ b/src/PJ_calcofi.c
@@ -0,0 +1,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)
diff --git a/src/makefile.vc b/src/makefile.vc
index e6d17c57..a420c926 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -38,7 +38,8 @@ pseudo = \
 	PJ_nell.obj PJ_nell_h.obj PJ_putp2.obj PJ_putp3.obj \
 	PJ_putp4p.obj PJ_putp5.obj PJ_putp6.obj PJ_robin.obj \
 	PJ_sts.obj PJ_urm5.obj PJ_urmfps.obj PJ_wag2.obj \
-	PJ_wag3.obj PJ_wink1.obj PJ_wink2.obj PJ_isea.obj 
+	PJ_wag3.obj PJ_wink1.obj PJ_wink2.obj PJ_isea.obj \
+	PJ_calcofi.obj 
 
 support = \
 	aasincos.obj adjlon.obj bch2bps.obj bchgen.obj pj_gauss.obj \
diff --git a/src/pj_list.h b/src/pj_list.h
index 8230f92f..95953758 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -16,6 +16,7 @@ PROJ_HEAD(bacon, "Bacon Globular")
 PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
 PROJ_HEAD(boggs, "Boggs Eumorphic")
 PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")
+PROJ_HEAD(calcofi, "Cal Coop Ocean Fish Invest Lines/Stations")
 PROJ_HEAD(cass, "Cassini")
 PROJ_HEAD(cc, "Central Cylindrical")
 PROJ_HEAD(cea, "Equal Area Cylindrical")