aboutsummaryrefslogtreecommitdiff
path: root/src/projections/sch.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/projections/sch.cpp')
-rw-r--r--src/projections/sch.cpp232
1 files changed, 232 insertions, 0 deletions
diff --git a/src/projections/sch.cpp b/src/projections/sch.cpp
new file mode 100644
index 00000000..5a2f944b
--- /dev/null
+++ b/src/projections/sch.cpp
@@ -0,0 +1,232 @@
+/******************************************************************************
+ * 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.
+ ****************************************************************************/
+
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+
+#include "proj.h"
+#include "projects.h"
+#include "geocent.h"
+
+namespace { // anonymous namespace
+struct pj_opaque {
+ 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;
+};
+} // anonymous namespace
+
+PROJ_HEAD(sch, "Spherical Cross-track Height") "\n\tMisc\n\tplat_0= plon_0= phdg_0= [h_0=]";
+
+static LPZ inverse3d(XYZ xyz, PJ *P) {
+ LPZ lpz = {0.0, 0.0, 0.0};
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ double temp[3];
+ double pxyz[3];
+
+ /* Local lat,lon using radius */
+ pxyz[0] = xyz.y * P->a / Q->rcurv;
+ pxyz[1] = xyz.x * P->a / Q->rcurv;
+ pxyz[2] = xyz.z;
+
+ if( pj_Convert_Geodetic_To_Geocentric( &(Q->sph), pxyz[0], pxyz[1], pxyz[2], temp, temp+1, temp+2) != 0) {
+ proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+ return lpz;
+ }
+
+ /* Apply rotation */
+ pxyz[0] = Q->transMat[0] * temp[0] + Q->transMat[1] * temp[1] + Q->transMat[2] * temp[2];
+ pxyz[1] = Q->transMat[3] * temp[0] + Q->transMat[4] * temp[1] + Q->transMat[5] * temp[2];
+ pxyz[2] = Q->transMat[6] * temp[0] + Q->transMat[7] * temp[1] + Q->transMat[8] * temp[2];
+
+ /* Apply offset */
+ pxyz[0] += Q->xyzoff[0];
+ pxyz[1] += Q->xyzoff[1];
+ pxyz[2] += Q->xyzoff[2];
+
+ /* Convert geocentric coordinates to lat lon */
+ pj_Convert_Geocentric_To_Geodetic( &(Q->elp_0), pxyz[0], pxyz[1], pxyz[2],
+ temp, temp+1, temp+2);
+
+
+ lpz.lam = temp[1] ;
+ lpz.phi = temp[0] ;
+ lpz.z = temp[2];
+
+ return lpz;
+}
+
+static XYZ forward3d(LPZ lpz, PJ *P) {
+ XYZ xyz = {0.0, 0.0, 0.0};
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ double temp[3];
+ double pxyz[3];
+
+
+ /* Convert lat lon to geocentric coordinates */
+ if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), lpz.phi, lpz.lam, lpz.z, temp, temp+1, temp+2 ) != 0 ) {
+ proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+ return xyz;
+ }
+
+
+ /* Adjust for offset */
+ temp[0] -= Q->xyzoff[0];
+ temp[1] -= Q->xyzoff[1];
+ temp[2] -= Q->xyzoff[2];
+
+
+ /* Apply rotation */
+ pxyz[0] = Q->transMat[0] * temp[0] + Q->transMat[3] * temp[1] + Q->transMat[6] * temp[2];
+ pxyz[1] = Q->transMat[1] * temp[0] + Q->transMat[4] * temp[1] + Q->transMat[7] * temp[2];
+ pxyz[2] = Q->transMat[2] * temp[0] + Q->transMat[5] * temp[1] + Q->transMat[8] * temp[2];
+
+ /* Convert to local lat,lon */
+ pj_Convert_Geocentric_To_Geodetic( &(Q->sph), pxyz[0], pxyz[1], pxyz[2],
+ temp, temp+1, temp+2);
+
+
+ /* Scale by radius */
+ xyz.x = temp[1] * Q->rcurv / P->a;
+ xyz.y = temp[0] * Q->rcurv / P->a;
+ xyz.z = temp[2];
+
+ return xyz;
+}
+
+
+static PJ *setup(PJ *P) { /* general initialization */
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
+ 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(&(Q->elp_0), P->a, temp) != 0)
+ return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+
+ clt = cos(Q->plat);
+ slt = sin(Q->plat);
+ clo = cos(Q->plon);
+ slo = sin(Q->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(Q->phdg);
+ shdg = sin(Q->phdg);
+
+ Q->rcurv = Q->h0 + (reast*rnorth)/(reast * chdg * chdg + rnorth * shdg * shdg);
+
+ /* Set up local sphere at the given peg point */
+ if ( pj_Set_Geocentric_Parameters(&(Q->sph), Q->rcurv, Q->rcurv) != 0)
+ return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+
+ /* Set up the transformation matrices */
+ Q->transMat[0] = clt * clo;
+ Q->transMat[1] = -shdg*slo - slt*clo * chdg;
+ Q->transMat[2] = slo*chdg - slt*clo*shdg;
+ Q->transMat[3] = clt*slo;
+ Q->transMat[4] = clo*shdg - slt*slo*chdg;
+ Q->transMat[5] = -clo*chdg - slt*slo*shdg;
+ Q->transMat[6] = slt;
+ Q->transMat[7] = clt*chdg;
+ Q->transMat[8] = clt*shdg;
+
+
+ if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), Q->plat, Q->plon, Q->h0,
+ pxyz, pxyz+1, pxyz+2 ) != 0 )
+ return pj_default_destructor(P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
+
+
+ Q->xyzoff[0] = pxyz[0] - (Q->rcurv) * clt * clo;
+ Q->xyzoff[1] = pxyz[1] - (Q->rcurv) * clt * slo;
+ Q->xyzoff[2] = pxyz[2] - (Q->rcurv) * slt;
+
+ P->fwd3d = forward3d;
+ P->inv3d = inverse3d;
+ return P;
+}
+
+
+PJ *PROJECTION(sch) {
+ struct pj_opaque *Q = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof (struct pj_opaque)));
+ if (nullptr==Q)
+ return pj_default_destructor(P, ENOMEM);
+ P->opaque = Q;
+
+ Q->h0 = 0.0;
+
+ /* Check if peg latitude was defined */
+ if (pj_param(P->ctx, P->params, "tplat_0").i)
+ Q->plat = pj_param(P->ctx, P->params, "rplat_0").f;
+ else {
+ return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+ }
+
+ /* Check if peg longitude was defined */
+ if (pj_param(P->ctx, P->params, "tplon_0").i)
+ Q->plon = pj_param(P->ctx, P->params, "rplon_0").f;
+ else {
+ return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+ }
+
+ /* Check if peg latitude is defined */
+ if (pj_param(P->ctx, P->params, "tphdg_0").i)
+ Q->phdg = pj_param(P->ctx, P->params, "rphdg_0").f;
+ else {
+ return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+ }
+
+
+ /* Check if average height was defined - If so read it in */
+ if (pj_param(P->ctx, P->params, "th_0").i)
+ Q->h0 = pj_param(P->ctx, P->params, "dh_0").f;
+
+
+ return setup(P);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 16:26:28 +0000