From 757a2c8f946faccf9d094d76cb79e6ebe0006564 Mon Sep 17 00:00:00 2001 From: Piyush Agram Date: Thu, 10 Sep 2015 21:44:23 -0700 Subject: - API Extension to include fwd3d / inv3d - New SCH coordinate system for radar imaging systems --- nad/testvarious | 26 + nad/tv_out.dist | 14 + src/Makefile.am | 4 +- src/PJ_sch.c | 228 ++++++ src/lib_proj.cmake | 3 + src/makefile.vc | 7 +- src/pj_fwd.c | 19 +- src/pj_fwd3d.c | 47 ++ src/pj_inv.c | 17 +- src/pj_inv3d.c | 45 ++ src/pj_list.h | 2 + src/pj_transform.c | 181 +++-- src/proj.def | 2 + src/proj_api.h | 8 + src/projects.h | 12 + wince/msvc80/projce_dll/projce_dll.vcproj | 1130 +++++++---------------------- wince/msvc80/projce_lib/projce_lib.vcproj | 1098 +++++++--------------------- 17 files changed, 1111 insertions(+), 1732 deletions(-) create mode 100644 src/PJ_sch.c create mode 100644 src/pj_fwd3d.c create mode 100644 src/pj_inv3d.c diff --git a/nad/testvarious b/nad/testvarious index 694d0d19..73a366fd 100755 --- a/nad/testvarious +++ b/nad/testvarious @@ -687,6 +687,32 @@ $EXE -f '%.7f' \ -E >>${OUT} <> ${OUT} +echo "Test SCH forward projection" >> ${OUT} +# +$EXE -f '%.7f' \ + +proj=latlong +datum=WGS84 +to +proj=sch +datum=WGS84 +plat_0=30.0 +plon_0=45.0 \ + +phdg_0=-12.0 +nodefs \ + -E >> ${OUT} <> ${OUT} +echo "Test SCH inverse projection" >> ${OUT} +# +$EXE -f '%.7f' \ + +proj=sch +datum=WGS84 +plat_0=30.0 +plon_0=45.0 +phdg_0=-12.0 +nodefs +to \ + +proj=latlong +datum=WGS84 \ + -E >> ${OUT} < + +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 +#include +# 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 +#include +# 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 99faeafb..fd2d3bdf 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 fe6a44c6..bbc6f2ca 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 diff --git a/wince/msvc80/projce_dll/projce_dll.vcproj b/wince/msvc80/projce_dll/projce_dll.vcproj index 37102afc..fceb5802 100644 --- a/wince/msvc80/projce_dll/projce_dll.vcproj +++ b/wince/msvc80/projce_dll/projce_dll.vcproj @@ -1,987 +1,411 @@ - - - - - - - + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + - + - + - + - + - + - + - + - + - + - + - + - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/wince/msvc80/projce_lib/projce_lib.vcproj b/wince/msvc80/projce_lib/projce_lib.vcproj index 13f00e35..c5c962d7 100644 --- a/wince/msvc80/projce_lib/projce_lib.vcproj +++ b/wince/msvc80/projce_lib/projce_lib.vcproj @@ -1,955 +1,411 @@ - - - - - - - + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + + + + + - + - + - + - + - + - + - + - + - + - + - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + -- cgit v1.2.3