Merge branch 'master' of https://github.com/piyushrpt/proj.4

Conflicts:
	nad/tv_out.dist

13 files changed, 517 insertions, 58 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 5d89cc28..83c5f34e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -54,12 +54,12 @@ libproj_la_SOURCES = \
 	PJ_mbt_fps.c PJ_mbtfpp.c PJ_mbtfpq.c PJ_moll.c \
 	PJ_nell.c PJ_nell_h.c PJ_putp2.c PJ_putp3.c \
 	PJ_putp4p.c PJ_putp5.c PJ_putp6.c PJ_qsc.c PJ_robin.c \
-	PJ_sts.c PJ_urm5.c PJ_urmfps.c PJ_wag2.c \
+	PJ_sch.c PJ_sts.c PJ_urm5.c PJ_urmfps.c PJ_wag2.c \
 	PJ_wag3.c PJ_wink1.c PJ_wink2.c pj_latlong.c pj_geocent.c \
 	aasincos.c adjlon.c bch2bps.c bchgen.c \
 	biveval.c dmstor.c mk_cheby.c pj_auth.c \
 	pj_deriv.c pj_ell_set.c pj_ellps.c pj_errno.c \
-	pj_factors.c pj_fwd.c pj_init.c pj_inv.c \
+	pj_factors.c pj_fwd.c pj_init.c pj_inv.c pj_fwd3d.c pj_inv3d.c\
 	pj_list.c pj_malloc.c pj_mlfn.c pj_msfn.c proj_mdist.c \
 	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 \
diff --git a/src/PJ_sch.c b/src/PJ_sch.c
new file mode 100644
index 00000000..5571deb3
--- /dev/null
+++ b/src/PJ_sch.c
@@ -0,0 +1,228 @@
+/******************************************************************************
+ * $Id$
+ *
+ * Project:  SCH Coordinate system
+ * Purpose:  Implementation of SCH Coordinate system
+ * References : 
+ *      1. Hensley. Scott. SCH Coordinates and various transformations. June 15, 2000.  
+ *      2. Buckley, Sean Monroe. Radar interferometry measurement of land subsidence. 2000..
+ *         PhD Thesis. UT Austin. (Appendix)
+ *      3. Hensley, Scott, Elaine Chapin, and T. Michel. "Improved processing of AIRSAR
+ *         data based on the GeoSAR processor." Airsar earth science and applications
+ *         workshop, March. 2002. (http://airsar.jpl.nasa.gov/documents/workshop2002/papers/T3.pdf)
+ *
+ * Author:   Piyush Agram (piyush.agram@jpl.nasa.gov)
+ * Copyright (c) 2015 California Institute of Technology.
+ * Government sponsorship acknowledged.
+ *
+ * NOTE:  The SCH coordinate system is a sensor aligned coordinate system 
+ * developed at JPL for radar mapping missions. Details pertaining to the 
+ * coordinate system have been release in the public domain (see references above).
+ * This code is an independent implementation of the SCH coordinate system
+ * that conforms to the PROJ.4 conventions and uses the details presented in these
+ * publicly released documents. All credit for the development of the coordinate
+ * system and its use should be directed towards the original developers at JPL.
+ ******************************************************************************
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include "geocent.h"
+
+#define PROJ_PARMS__ \
+	double plat; /*Peg Latitude */    \
+        double plon; /*Peg Longitude*/    \
+        double phdg; /*Peg heading  */    \
+        double h0;   /*Average altitude */\
+	double transMat[9]; \
+        double xyzoff[3]; \
+        double rcurv; \
+        GeocentricInfo sph; \
+        GeocentricInfo elp_0; 
+        
+#define PJ_LIB__
+#include	<projects.h>
+
+PROJ_HEAD(sch, "Spherical Cross-track Height") "\n\tMisc\n\tplat_0 = ,plon_0 = , phdg_0 = ,[h_0 = ]";
+
+INVERSE3D(inverse3d);
+        double temp[3];
+        double pxyz[3];
+
+        //Local lat,lon using radius
+        pxyz[0] = xyz.y * P->a / P->rcurv;
+        pxyz[1] = xyz.x * P->a / P->rcurv;
+        pxyz[2] = xyz.z;
+
+
+        if( pj_Convert_Geodetic_To_Geocentric( &(P->sph), pxyz[0], pxyz[1], pxyz[2],
+                    temp, temp+1, temp+2) != 0)
+                I3_ERROR;
+
+        //Apply rotation
+        pxyz[0] = P->transMat[0] * temp[0] + P->transMat[1] * temp[1] + P->transMat[2] * temp[2];
+        pxyz[1] = P->transMat[3] * temp[0] + P->transMat[4] * temp[1] + P->transMat[5] * temp[2];
+        pxyz[2] = P->transMat[6] * temp[0] + P->transMat[7] * temp[1] + P->transMat[8] * temp[2];
+
+        //Apply offset
+        pxyz[0] += P->xyzoff[0];
+        pxyz[1] += P->xyzoff[1];
+        pxyz[2] += P->xyzoff[2];
+
+        //Convert geocentric coordinates to lat lon
+        pj_Convert_Geocentric_To_Geodetic( &(P->elp_0), pxyz[0], pxyz[1], pxyz[2],
+                temp, temp+1, temp+2);
+
+
+        lpz.lam = temp[0] ;
+        lpz.phi = temp[1] ;
+        lpz.z = temp[2];
+
+//        printf("INVERSE: \n");
+//        printf("XYZ: %f %f %f \n", xyz.x, xyz.y, xyz.z);
+//        printf("LPZ: %f %f %f \n", lpz.lam, lpz.phi, lpz.z);
+        return (lpz);
+}
+
+FORWARD3D(forward3d);
+        double temp[3];
+        double pxyz[3];
+
+        //Convert lat lon to geocentric coordinates
+        if( pj_Convert_Geodetic_To_Geocentric( &(P->elp_0), lpz.lam, lpz.phi, lpz.z,
+                    temp, temp+1, temp+2 ) != 0 )
+            F3_ERROR;
+
+        //Adjust for offset
+        temp[0] -= P->xyzoff[0];
+        temp[1] -= P->xyzoff[1];
+        temp[2] -= P->xyzoff[2];
+
+       
+        //Apply rotation 
+        pxyz[0] = P->transMat[0] * temp[0] + P->transMat[3] * temp[1] + P->transMat[6] * temp[2];
+        pxyz[1] = P->transMat[1] * temp[0] + P->transMat[4] * temp[1] + P->transMat[7] * temp[2];
+        pxyz[2] = P->transMat[2] * temp[0] + P->transMat[5] * temp[1] + P->transMat[8] * temp[2];
+
+        //Convert to local lat,lon
+        pj_Convert_Geocentric_To_Geodetic( &(P->sph), pxyz[0], pxyz[1], pxyz[2],
+                temp, temp+1, temp+2);
+
+        //Scale by radius
+        xyz.x = temp[1] * P->rcurv / P->a;
+        xyz.y = temp[0] * P->rcurv / P->a;
+        xyz.z = temp[2];
+
+//        printf("FORWARD: \n");
+//        printf("LPZ: %f %f %f \n", lpz.lam, lpz.phi, lpz.z);
+//        printf("XYZ: %f %f %f \n", xyz.x, xyz.y, xyz.z);
+
+        return (xyz);
+}
+
+FREEUP; if (P) pj_dalloc(P); }
+	static PJ *
+setup(PJ *P) { /* general initialization */
+        double reast, rnorth;
+        double chdg, shdg;
+        double clt, slt;
+        double clo, slo;
+        double temp;
+        double pxyz[3];
+
+        temp = P->a * sqrt(1.0 - P->es);
+
+        //Setup original geocentric system
+        if ( pj_Set_Geocentric_Parameters(&(P->elp_0), P->a, temp) != 0)
+                E_ERROR(-37);
+
+
+        clt = cos(P->plat);
+        slt = sin(P->plat);
+        clo = cos(P->plon);
+        slo = sin(P->plon);
+
+        //Estimate the radius of curvature for given peg
+        temp = sqrt(1.0 - (P->es) * slt * slt);
+        reast = (P->a)/temp;
+        rnorth = (P->a) * (1.0 - (P->es))/pow(temp,3);
+
+        chdg = cos(P->phdg);
+        shdg = sin(P->phdg);
+        
+        P->rcurv = P->h0 + (reast*rnorth)/(reast * chdg * chdg + rnorth * shdg * shdg);
+
+//        printf("North Radius: %f \n", rnorth);
+//        printf("East Radius: %f \n", reast);
+//        printf("Effective Radius: %f \n", P->rcurv);
+
+        //Set up local sphere at the given peg point
+        if ( pj_Set_Geocentric_Parameters(&(P->sph), P->rcurv, P->rcurv) != 0)
+            E_ERROR(-37);
+
+        //Set up the transformation matrices
+        P->transMat[0] = clt * clo; 
+        P->transMat[1] = -shdg*slo - slt*clo * chdg;
+        P->transMat[2] =  slo*chdg - slt*clo*shdg;
+        P->transMat[3] =  clt*slo;
+        P->transMat[4] =  clo*shdg - slt*slo*chdg;
+        P->transMat[5] = -clo*chdg - slt*slo*shdg;
+        P->transMat[6] =  slt;
+        P->transMat[7] =  clt*chdg;
+        P->transMat[8] =  clt*shdg;        
+
+        
+        if( pj_Convert_Geodetic_To_Geocentric( &(P->elp_0), P->plat, P->plon, P->h0, 
+                                               pxyz, pxyz+1, pxyz+2 ) != 0 )
+        {
+            E_ERROR(-14)
+        }
+
+
+        P->xyzoff[0] = pxyz[0] - (P->rcurv) * clt * clo;
+        P->xyzoff[1] = pxyz[1] - (P->rcurv) * clt * slo;
+        P->xyzoff[2] = pxyz[2] - (P->rcurv) * slt;
+
+//        printf("Offset: %f %f %f \n", P->xyzoff[0], P->xyzoff[1], P->xyzoff[2]);
+
+
+        P->fwd3d = forward3d;
+        P->inv3d = inverse3d;
+	return P;
+}
+ENTRY0(sch)
+        P->h0 = 0.0;
+
+        //Check if peg latitude was defined
+        if (pj_param(P->ctx, P->params, "tplat_0").i)
+            P->plat = pj_param(P->ctx, P->params, "rplat_0").f;
+        else
+            E_ERROR(-37);
+
+        //Check if peg longitude was defined
+        if (pj_param(P->ctx, P->params, "tplon_0").i)
+            P->plon = pj_param(P->ctx, P->params, "rplon_0").f;
+        else
+            E_ERROR(-37);
+        
+        //Check if peg latitude is defined
+        if (pj_param(P->ctx, P->params, "tphdg_0").i)
+            P->phdg = pj_param(P->ctx, P->params, "rphdg_0").f;
+        else
+            E_ERROR(-37);
+        
+
+        //Check if average height was defined
+        //If so read it in
+        if (pj_param(P->ctx, P->params, "th_0").i)
+            P->h0 = pj_param(P->ctx, P->params, "dh_0").f;
+
+        //Completed reading in the projection parameters
+//        printf("PSA: Lat = %f Lon = %f Hdg = %f \n", P->plat, P->plon, P->phdg);
+
+ENDENTRY(setup(P))
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 8a6349d8..b7b497cc 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -115,6 +115,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_qsc.c
         PJ_robin.c
         PJ_rpoly.c
+        PJ_sch.c
         PJ_sconics.c
         PJ_somerc.c
         PJ_sterea.c
@@ -166,6 +167,7 @@ SET(SRC_LIBPROJ_CORE
         pj_errno.c
         pj_factors.c
         pj_fwd.c
+        pj_fwd3d.c
         pj_gauss.c
         pj_gc_reader.c
         pj_geocent.c
@@ -176,6 +178,7 @@ SET(SRC_LIBPROJ_CORE
         pj_init.c
         pj_initcache.c
         pj_inv.c
+        pj_inv3d.c
         pj_latlong.c
         pj_list.c
         pj_list.h
diff --git a/src/makefile.vc b/src/makefile.vc
index 41b7dca8..ed7ab4c4 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -24,9 +24,9 @@ misc = \
 	PJ_airy.obj PJ_aitoff.obj PJ_august.obj PJ_bacon.obj \
 	PJ_chamb.obj PJ_hammer.obj PJ_lagrng.obj PJ_larr.obj \
 	PJ_lask.obj PJ_nocol.obj PJ_ob_tran.obj PJ_oea.obj \
-	PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj PJ_vandg4.obj \
-	PJ_wag7.obj pj_latlong.obj PJ_krovak.obj pj_geocent.obj \
-	PJ_healpix.obj PJ_natearth.obj PJ_qsc.obj
+	PJ_sch.obj PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj \
+        PJ_vandg4.obj PJ_wag7.obj pj_latlong.obj PJ_krovak.obj \
+        pj_geocent.obj PJ_healpix.obj PJ_natearth.obj PJ_qsc.obj
 
 pseudo = \
 	PJ_boggs.obj PJ_collg.obj PJ_crast.obj PJ_denoy.obj \
@@ -45,6 +45,7 @@ support = \
 	biveval.obj dmstor.obj mk_cheby.obj pj_auth.obj \
 	pj_deriv.obj pj_ell_set.obj pj_ellps.obj pj_errno.obj \
 	pj_factors.obj pj_fwd.obj pj_init.obj pj_inv.obj \
+        pj_fwd3d.obj pj_inv3d.obj \
 	pj_list.obj pj_malloc.obj pj_mlfn.obj pj_msfn.obj \
 	pj_open_lib.obj pj_param.obj pj_phi2.obj pj_pr_list.obj \
 	pj_qsfn.obj pj_strerrno.obj pj_tsfn.obj pj_units.obj \
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index b70b4240..1cd002b5 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -24,14 +24,23 @@ pj_fwd(LP lp, PJ *P) {
 		lp.lam -= P->lam0;	/* compute del lp.lam */
 		if (!P->over)
 			lp.lam = adjlon(lp.lam); /* adjust del longitude */
-		xy = (*P->fwd)(lp, P); /* project */
-		if ( P->ctx->last_errno )
+
+                //Check for NULL pointer
+                if (P->fwd != NULL)
+                {
+		    xy = (*P->fwd)(lp, P); /* project */
+		    if ( P->ctx->last_errno )
 			xy.x = xy.y = HUGE_VAL;
-		/* adjust for major axis and easting/northings */
-		else {
+		    /* adjust for major axis and easting/northings */
+		    else {
 			xy.x = P->fr_meter * (P->a * xy.x + P->x0);
 			xy.y = P->fr_meter * (P->a * xy.y + P->y0);
-		}
+		    }
+                }
+                else
+                {
+                    xy.x = xy.y = HUGE_VAL;
+                }
 	}
 	return xy;
 }
diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c
new file mode 100644
index 00000000..834746d4
--- /dev/null
+++ b/src/pj_fwd3d.c
@@ -0,0 +1,47 @@
+/* general forward projection */
+#define PJ_LIB__
+#include <projects.h>
+#include <errno.h>
+# define EPS 1.0e-12
+	XYZ /* forward projection entry */
+pj_fwd3d(LPZ lpz, PJ *P) {
+	XYZ xyz;
+	double t;
+
+	/* check for forward and latitude or longitude overange */
+	if ((t = fabs(lpz.phi)-HALFPI) > EPS || fabs(lpz.lam) > 10.) {
+		xyz.x = xyz.y = xyz.z = HUGE_VAL;
+		pj_ctx_set_errno( P->ctx, -14);
+	} else { /* proceed with projection */
+                P->ctx->last_errno = 0;
+                pj_errno = 0;
+                errno = 0;
+
+		if (fabs(t) <= EPS)
+			lpz.phi = lpz.phi < 0. ? -HALFPI : HALFPI;
+		else if (P->geoc)   //Maybe redundant and never used.
+			lpz.phi = atan(P->rone_es * tan(lpz.phi));
+		lpz.lam -= P->lam0;	/* compute del lp.lam */
+		if (!P->over)
+			lpz.lam = adjlon(lpz.lam); /* adjust del longitude */
+
+                //Check for NULL pointer
+                if (P->fwd3d != NULL)
+                {
+		    xyz = (*P->fwd3d)(lpz, P); /* project */
+		    if ( P->ctx->last_errno )
+			xyz.x = xyz.y = xyz.z = HUGE_VAL;
+		    /* adjust for major axis and easting/northings */
+		    else {
+			xyz.x = P->fr_meter * (P->a * xyz.x + P->x0);
+			xyz.y = P->fr_meter * (P->a * xyz.y + P->y0);
+                        //z is not scaled since this handled by vto_meter outside
+		    }
+                }
+                else
+                {
+                    xyz.x = xyz.y = xyz.z = HUGE_VAL;
+                }
+	}
+	return xyz;
+}
diff --git a/src/pj_inv.c b/src/pj_inv.c
index a418ccd4..d77b4e56 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -19,15 +19,24 @@ pj_inv(XY xy, PJ *P) {
 
 	xy.x = (xy.x * P->to_meter - P->x0) * P->ra; /* descale and de-offset */
 	xy.y = (xy.y * P->to_meter - P->y0) * P->ra;
-	lp = (*P->inv)(xy, P); /* inverse project */
-	if (P->ctx->last_errno )
+
+        //Check for NULL pointer
+        if (P->inv != NULL)
+        {
+	    lp = (*P->inv)(xy, P); /* inverse project */
+	    if (P->ctx->last_errno )
 		lp.lam = lp.phi = HUGE_VAL;
-	else {
+	    else {
 		lp.lam += P->lam0; /* reduce from del lp.lam */
 		if (!P->over)
 			lp.lam = adjlon(lp.lam); /* adjust longitude to CM */
 		if (P->geoc && fabs(fabs(lp.phi)-HALFPI) > EPS)
 			lp.phi = atan(P->one_es * tan(lp.phi));
-	}
+	    }
+        }
+        else
+        {
+           lp.lam = lp.phi = HUGE_VAL;
+        }
 	return lp;
 }
diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c
new file mode 100644
index 00000000..de35e776
--- /dev/null
+++ b/src/pj_inv3d.c
@@ -0,0 +1,45 @@
+/* general inverse projection */
+#define PJ_LIB__
+#include <projects.h>
+#include <errno.h>
+# define EPS 1.0e-12
+	LPZ /* inverse projection entry */
+pj_inv3d(XYZ xyz, PJ *P) {
+	LPZ lpz;
+
+	/* can't do as much preliminary checking as with forward */
+	if (xyz.x == HUGE_VAL || xyz.y == HUGE_VAL || xyz.z == HUGE_VAL ) {
+		lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
+		pj_ctx_set_errno( P->ctx, -15);
+                return lpz;
+	}
+
+	errno = pj_errno = 0;
+        P->ctx->last_errno = 0;
+
+	xyz.x = (xyz.x * P->to_meter - P->x0) * P->ra; /* descale and de-offset */
+	xyz.y = (xyz.y * P->to_meter - P->y0) * P->ra;
+        //z is not scaled since that is handled by vto_meter before we get here
+        
+        //Check for NULL pointer
+        if (P->inv3d != NULL)
+        {
+	    lpz = (*P->inv3d)(xyz, P); /* inverse project */
+	    if (P->ctx->last_errno )
+		lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
+	    else {
+		lpz.lam += P->lam0; /* reduce from del lp.lam */
+		if (!P->over)
+			lpz.lam = adjlon(lpz.lam); /* adjust longitude to CM */
+
+                //This maybe redundant and never user
+		if (P->geoc && fabs(fabs(lpz.phi)-HALFPI) > EPS)
+			lpz.phi = atan(P->one_es * tan(lpz.phi));
+	    }
+        }
+        else
+        {
+            lpz.lam = lpz.phi = lpz.z = HUGE_VAL;
+        }
+	return lpz;
+}
diff --git a/src/pj_list.h b/src/pj_list.h
index 194bcc68..20af2eb2 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -111,6 +111,7 @@ PROJ_HEAD(qsc, "Quadrilateralized Spherical Cube")
 PROJ_HEAD(robin, "Robinson")
 PROJ_HEAD(rouss, "Roussilhe Stereographic")
 PROJ_HEAD(rpoly, "Rectangular Polyconic")
+PROJ_HEAD(sch, "Spherical Cross-track Height")
 PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)")
 PROJ_HEAD(somerc, "Swiss. Obl. Mercator")
 PROJ_HEAD(stere, "Stereographic")
@@ -142,3 +143,4 @@ PROJ_HEAD(weren, "Werenskiold I")
 PROJ_HEAD(wink1, "Winkel I")
 PROJ_HEAD(wink2, "Winkel II")
 PROJ_HEAD(wintri, "Winkel Tripel")
+
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 6fe571c6..32f14955 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -152,7 +152,9 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
 /* -------------------------------------------------------------------- */
     else if( !srcdefn->is_latlong )
     {
-        if( srcdefn->inv == NULL )
+
+        //Check first if projection is invertible.
+        if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL))
         {
             pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
             pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR, 
@@ -160,35 +162,85 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
             return -17;
         }
 
-        for( i = 0; i < point_count; i++ )
+        //If invertible - First try inv3d if defined
+        if (srcdefn->inv3d != NULL)
         {
-            XY         projected_loc;
-            LP	       geodetic_loc;
+            //Three dimensions must be defined
+            if ( z == NULL)
+            {
+                pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
+                return PJD_ERR_GEOCENTRIC;
+            }
 
-            projected_loc.u = x[point_offset*i];
-            projected_loc.v = y[point_offset*i];
+            for (i=0; i < point_count; i++)
+            {
+                XYZ projected_loc;
+                XYZ geodetic_loc;
 
-            if( projected_loc.u == HUGE_VAL )
-                continue;
+                projected_loc.u = x[point_offset*i];
+                projected_loc.v = y[point_offset*i];
+                projected_loc.w = z[point_offset*i];
 
-            geodetic_loc = pj_inv( projected_loc, srcdefn );
-            if( srcdefn->ctx->last_errno != 0 )
-            {
-                if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ 
-                     && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
-                    && (srcdefn->ctx->last_errno > 0 
-                        || srcdefn->ctx->last_errno < -44 || point_count == 1
-                        || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
-                    return srcdefn->ctx->last_errno;
-                else
+                if (projected_loc.u == HUGE_VAL)
+                    continue;
+
+                geodetic_loc = pj_inv3d(projected_loc, srcdefn);
+                if( srcdefn->ctx->last_errno != 0 )
                 {
-                    geodetic_loc.u = HUGE_VAL;
-                    geodetic_loc.v = HUGE_VAL;
+                    if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ 
+                         && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
+                        && (srcdefn->ctx->last_errno > 0 
+                            || srcdefn->ctx->last_errno < -44 || point_count == 1
+                            || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
+                        return srcdefn->ctx->last_errno;
+                    else
+                    {
+                        geodetic_loc.u = HUGE_VAL;
+                        geodetic_loc.v = HUGE_VAL;
+                        geodetic_loc.w = HUGE_VAL;
+                    }
                 }
+
+                x[point_offset*i] = geodetic_loc.u;
+                y[point_offset*i] = geodetic_loc.v;
+                z[point_offset*i] = geodetic_loc.w;
+
             }
 
-            x[point_offset*i] = geodetic_loc.u;
-            y[point_offset*i] = geodetic_loc.v;
+        }
+        else
+        {
+            //Fallback to the original PROJ.4 API 2d inversion- inv 
+            for( i = 0; i < point_count; i++ )
+            {
+                XY         projected_loc;
+                LP	       geodetic_loc;
+
+                projected_loc.u = x[point_offset*i];
+                projected_loc.v = y[point_offset*i];
+
+                if( projected_loc.u == HUGE_VAL )
+                    continue;
+
+                geodetic_loc = pj_inv( projected_loc, srcdefn );
+                if( srcdefn->ctx->last_errno != 0 )
+                {
+                    if( (srcdefn->ctx->last_errno != 33 /*EDOM*/ 
+                         && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
+                        && (srcdefn->ctx->last_errno > 0 
+                            || srcdefn->ctx->last_errno < -44 || point_count == 1
+                            || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
+                        return srcdefn->ctx->last_errno;
+                    else
+                    {
+                        geodetic_loc.u = HUGE_VAL;
+                        geodetic_loc.v = HUGE_VAL;
+                    }
+                }
+
+                x[point_offset*i] = geodetic_loc.u;
+                y[point_offset*i] = geodetic_loc.v;
+            }
         }
     }
 /* -------------------------------------------------------------------- */
@@ -289,35 +341,76 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
 /* -------------------------------------------------------------------- */
     else if( !dstdefn->is_latlong )
     {
-        for( i = 0; i < point_count; i++ )
+
+        if( dstdefn->fwd3d != NULL)
         {
-            XY         projected_loc;
-            LP	       geodetic_loc;
+            for( i = 0; i < point_count; i++ )
+            {
+                XYZ projected_loc;
+                LPZ geodetic_loc;
+                
+                geodetic_loc.u = x[point_offset*i];
+                geodetic_loc.v = y[point_offset*i];
+                geodetic_loc.w = z[point_offset*i];
 
-            geodetic_loc.u = x[point_offset*i];
-            geodetic_loc.v = y[point_offset*i];
+                if (geodetic_loc.u == HUGE_VAL)
+                    continue;
 
-            if( geodetic_loc.u == HUGE_VAL )
-                continue;
+                projected_loc = pj_fwd3d( geodetic_loc, dstdefn);
+                if( dstdefn->ctx->last_errno != 0 )
+                {
+                    if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ 
+                         && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
+                        && (dstdefn->ctx->last_errno > 0 
+                            || dstdefn->ctx->last_errno < -44 || point_count == 1
+                            || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
+                        return dstdefn->ctx->last_errno;
+                    else
+                    {
+                        projected_loc.u = HUGE_VAL;
+                        projected_loc.v = HUGE_VAL;
+                        projected_loc.w = HUGE_VAL;
+                    }
+                }   
+
+                x[point_offset*i] = projected_loc.u;
+                y[point_offset*i] = projected_loc.v;
+                z[point_offset*i] = projected_loc.w;
+            }
 
-            projected_loc = pj_fwd( geodetic_loc, dstdefn );
-            if( dstdefn->ctx->last_errno != 0 )
+        }
+        else
+        {
+            for( i = 0; i < point_count; i++ )
             {
-                if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ 
-                     && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
-                    && (dstdefn->ctx->last_errno > 0 
-                        || dstdefn->ctx->last_errno < -44 || point_count == 1
-                        || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
-                    return dstdefn->ctx->last_errno;
-                else
+                XY         projected_loc;
+                LP	       geodetic_loc;
+
+                geodetic_loc.u = x[point_offset*i];
+                geodetic_loc.v = y[point_offset*i];
+
+                if( geodetic_loc.u == HUGE_VAL )
+                    continue;
+
+                projected_loc = pj_fwd( geodetic_loc, dstdefn );
+                if( dstdefn->ctx->last_errno != 0 )
                 {
-                    projected_loc.u = HUGE_VAL;
-                    projected_loc.v = HUGE_VAL;
-                }
+                    if( (dstdefn->ctx->last_errno != 33 /*EDOM*/ 
+                         && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
+                        && (dstdefn->ctx->last_errno > 0 
+                            || dstdefn->ctx->last_errno < -44 || point_count == 1
+                            || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
+                        return dstdefn->ctx->last_errno;
+                    else
+                    {
+                        projected_loc.u = HUGE_VAL;
+                        projected_loc.v = HUGE_VAL;
+                    }
+                }   
+
+                x[point_offset*i] = projected_loc.u;
+                y[point_offset*i] = projected_loc.v;
             }
-
-            x[point_offset*i] = projected_loc.u;
-            y[point_offset*i] = projected_loc.v;
         }
     }
 
diff --git a/src/proj.def b/src/proj.def
index f3146551..c872e6aa 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -74,3 +74,5 @@ EXPORTS
 	pj_open_lib             @72
 	pj_atof                 @73
 	pj_strtod               @74
+        pj_fwd3d                @75
+        pj_inv3d                @76
diff --git a/src/proj_api.h b/src/proj_api.h
index 9e1e49a8..e381815c 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -54,15 +54,20 @@ extern int pj_errno;	/* global error return code */
 
 #if !defined(PROJECTS_H)
     typedef struct { double u, v; } projUV;
+    typedef struct { double u, v, w; } projUVW;
     typedef void *projPJ;
     #define projXY projUV
     #define projLP projUV
+    #define projXYZ projUVW
+    #define projLPZ projUVW
     typedef void *projCtx;
 #else
     typedef PJ *projPJ;
     typedef projCtx_t *projCtx;
 #   define projXY	XY
 #   define projLP       LP
+#   define projXYZ      XYZ
+#   define projLPZ      LPZ
 #endif
 
 /* file reading api, like stdio */
@@ -80,6 +85,9 @@ typedef struct projFileAPI_t {
 projXY pj_fwd(projLP, projPJ);
 projLP pj_inv(projXY, projPJ);
 
+projXYZ pj_fwd3d(projLPZ, projPJ);
+projLPZ pj_inv3d(projXYZ, projPJ);
+
 int pj_transform( projPJ src, projPJ dst, long point_count, int point_offset,
                   double *x, double *y, double *z );
 int pj_datum_transform( projPJ src, projPJ dst, long point_count, int point_offset,
diff --git a/src/projects.h b/src/projects.h
index 43287849..0b1d5b3e 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -147,13 +147,18 @@ typedef struct {
 
 typedef struct { double u, v; } projUV;
 typedef struct { double r, i; }	COMPLEX;
+typedef struct { double u, v, w; } projUVW;
 
 #ifndef PJ_LIB__
 #define XY projUV
 #define LP projUV
+#define XYZ projUVW
+#define LPZ projUVW
 #else
 typedef struct { double x, y; }     XY;
 typedef struct { double lam, phi; } LP;
+typedef struct { double x, y, z; } XYZ;
+typedef struct { double lam, phi, z; } LPZ;
 #endif
 
 typedef union { double  f; int  i; char *s; } PROJVALUE;
@@ -225,6 +230,8 @@ typedef struct PJconsts {
     projCtx_t *ctx;
 	XY  (*fwd)(LP, struct PJconsts *);
 	LP  (*inv)(XY, struct PJconsts *);
+        XYZ (*fwd3d)(LPZ, struct PJconsts *);
+        LPZ (*inv3d)(XYZ, struct PJconsts *); 
 	void (*spc)(LP, struct PJconsts *, struct FACTORS *);
 	void (*pfree)(struct PJconsts *);
 	const char *descr;
@@ -324,6 +331,7 @@ extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[];
 	if( (P = (PJ*) pj_malloc(sizeof(PJ))) != NULL) { \
         memset( P, 0, sizeof(PJ) ); \
 	P->pfree = freeup; P->fwd = 0; P->inv = 0; \
+        P->fwd3d = 0; P->inv3d = 0; \
 	P->spc = 0; P->descr = des_##name;
 #define ENTRYX } return P; } else {
 #define ENTRY0(name) ENTRYA(name) ENTRYX
@@ -333,9 +341,13 @@ extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[];
 #define E_ERROR(err) { pj_ctx_set_errno( P->ctx, err); freeup(P); return(0); }
 #define E_ERROR_0 { freeup(P); return(0); }
 #define F_ERROR { pj_ctx_set_errno( P->ctx, -20); return(xy); }
+#define F3_ERROR { pj_ctx_set_errno( P->ctx, -20); return(xyz); }
 #define I_ERROR { pj_ctx_set_errno( P->ctx, -20); return(lp); }
+#define I3_ERROR { pj_ctx_set_errno( P->ctx, -20); return(lpz); }
 #define FORWARD(name) static XY name(LP lp, PJ *P) { XY xy = {0.0,0.0}
 #define INVERSE(name) static LP name(XY xy, PJ *P) { LP lp = {0.0,0.0}
+#define FORWARD3D(name) static XYZ name(LPZ lpz, PJ *P) {XYZ xyz = {0.0, 0.0, 0.0}
+#define INVERSE3D(name) static LPZ name(XYZ xyz, PJ *P) {LPZ lpz = {0.0, 0.0, 0.0}
 #define FREEUP static void freeup(PJ *P) {
 #define SPECIAL(name) static void name(LP lp, PJ *P, struct FACTORS *fac)
 #endif