/src/proj/src/projections/sch.cpp

Source (jump to first uncovered line)
/******************************************************************************
 * 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 <errno.h>
#include <math.h>

#include "proj.h"
#include "proj_internal.h"

namespace { // anonymous namespace
struct pj_sch_data {
    double plat; /*Peg Latitude */
    double plon; /*Peg Longitude*/
    double phdg; /*Peg heading  */
    double h0;   /*Average altitude */
    double transMat[9];
    double xyzoff[3];
    double rcurv;
    PJ *cart;
    PJ *cart_sph;
};
} // anonymous namespace

PROJ_HEAD(sch, "Spherical Cross-track Height")
"\n\tMisc\n\tplat_0= plon_0= phdg_0= [h_0=]";

static PJ_LPZ sch_inverse3d(PJ_XYZ xyz, PJ *P) {
    struct pj_sch_data *Q = static_cast<struct pj_sch_data *>(P->opaque);

    PJ_LPZ lpz;
    lpz.lam = xyz.x * (P->a / Q->rcurv);
    lpz.phi = xyz.y * (P->a / Q->rcurv);
    lpz.z = xyz.z;
    xyz = Q->cart_sph->fwd3d(lpz, Q->cart_sph);

    /* Apply rotation */
    xyz = {Q->transMat[0] * xyz.x + Q->transMat[1] * xyz.y +
               Q->transMat[2] * xyz.z,
           Q->transMat[3] * xyz.x + Q->transMat[4] * xyz.y +
               Q->transMat[5] * xyz.z,
           Q->transMat[6] * xyz.x + Q->transMat[7] * xyz.y +
               Q->transMat[8] * xyz.z};

    /* Apply offset */
    xyz.x += Q->xyzoff[0];
    xyz.y += Q->xyzoff[1];
    xyz.z += Q->xyzoff[2];

    /* Convert geocentric coordinates to lat long */
    return Q->cart->inv3d(xyz, Q->cart);
}

static PJ_XYZ sch_forward3d(PJ_LPZ lpz, PJ *P) {
    struct pj_sch_data *Q = static_cast<struct pj_sch_data *>(P->opaque);

    /* Convert lat long to geocentric coordinates */
    PJ_XYZ xyz = Q->cart->fwd3d(lpz, Q->cart);

    /* Adjust for offset */
    xyz.x -= Q->xyzoff[0];
    xyz.y -= Q->xyzoff[1];
    xyz.z -= Q->xyzoff[2];

    /* Apply rotation */
    xyz = {Q->transMat[0] * xyz.x + Q->transMat[3] * xyz.y +
               Q->transMat[6] * xyz.z,
           Q->transMat[1] * xyz.x + Q->transMat[4] * xyz.y +
               Q->transMat[7] * xyz.z,
           Q->transMat[2] * xyz.x + Q->transMat[5] * xyz.y +
               Q->transMat[8] * xyz.z};

    /* Convert to local lat,long */
    lpz = Q->cart_sph->inv3d(xyz, Q->cart_sph);

    /* Scale by radius */
    xyz.x = lpz.lam * (Q->rcurv / P->a);
    xyz.y = lpz.phi * (Q->rcurv / P->a);
    xyz.z = lpz.z;

    return xyz;
}

static PJ *pj_sch_destructor(PJ *P, int errlev) {
    if (nullptr == P)
        return nullptr;

    auto Q = static_cast<struct pj_sch_data *>(P->opaque);
    if (Q) {
        if (Q->cart)
            Q->cart->destructor(Q->cart, errlev);
        if (Q->cart_sph)
            Q->cart_sph->destructor(Q->cart_sph, errlev);
    }

    return pj_default_destructor(P, errlev);
}

static PJ *pj_sch_setup(PJ *P) { /* general initialization */
    struct pj_sch_data *Q = static_cast<struct pj_sch_data *>(P->opaque);

    /* Setup original geocentric system */
    // Pass a dummy ellipsoid definition that will be overridden just afterwards
    Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
    if (Q->cart == nullptr)
        return pj_sch_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);

    /* inherit ellipsoid definition from P to Q->cart */
    pj_inherit_ellipsoid_def(P, Q->cart);

    const double clt = cos(Q->plat);
    const double slt = sin(Q->plat);
    const double clo = cos(Q->plon);
    const double slo = sin(Q->plon);

    /* Estimate the radius of curvature for given peg */
    const double temp = sqrt(1.0 - (P->es) * slt * slt);
    const double reast = (P->a) / temp;
    const double rnorth = (P->a) * (1.0 - (P->es)) / pow(temp, 3);

    const double chdg = cos(Q->phdg);
    const double 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 */
    Q->cart_sph = proj_create(P->ctx, "+proj=cart +a=1");
    if (Q->cart_sph == nullptr)
        return pj_sch_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    pj_calc_ellipsoid_params(Q->cart_sph, Q->rcurv, 0);

    /* 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;

    PJ_LPZ lpz;
    lpz.lam = Q->plon;
    lpz.phi = Q->plat;
    lpz.z = Q->h0;
    PJ_XYZ xyz = Q->cart->fwd3d(lpz, Q->cart);
    Q->xyzoff[0] = xyz.x - (Q->rcurv) * clt * clo;
    Q->xyzoff[1] = xyz.y - (Q->rcurv) * clt * slo;
    Q->xyzoff[2] = xyz.z - (Q->rcurv) * slt;

    P->fwd3d = sch_forward3d;
    P->inv3d = sch_inverse3d;
    return P;
}

PJ *PJ_PROJECTION(sch) {
    struct pj_sch_data *Q = static_cast<struct pj_sch_data *>(
        calloc(1, sizeof(struct pj_sch_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_sch_destructor;

    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 {
        proj_log_error(P, _("Missing parameter plat_0."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    /* 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 {
        proj_log_error(P, _("Missing parameter plon_0."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    /* Check if peg heading is defined */
    if (pj_param(P->ctx, P->params, "tphdg_0").i)
        Q->phdg = pj_param(P->ctx, P->params, "rphdg_0").f;
    else {
        proj_log_error(P, _("Missing parameter phdg_0."));
        return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    /* 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 pj_sch_setup(P);
}

Coverage Report

Created: 2025-06-13 06:29

Line	Count	Source (jump to first uncovered line)
1		/******************************************************************************
2		* Project: SCH Coordinate system
3		* Purpose: Implementation of SCH Coordinate system
4		* References :
5		* 1. Hensley. Scott. SCH Coordinates and various transformations. June 15,
6		*2000.
7		* 2. Buckley, Sean Monroe. Radar interferometry measurement of land
8		*subsidence. 2000.. PhD Thesis. UT Austin. (Appendix)
9		* 3. Hensley, Scott, Elaine Chapin, and T. Michel. "Improved processing of
10		*AIRSAR data based on the GeoSAR processor." Airsar earth science and
11		*applications workshop, March. 2002.
12		*(http://airsar.jpl.nasa.gov/documents/workshop2002/papers/T3.pdf)
13		*
14		* Author: Piyush Agram (piyush.agram@jpl.nasa.gov)
15		* Copyright (c) 2015 California Institute of Technology.
16		* Government sponsorship acknowledged.
17		*
18		* NOTE: The SCH coordinate system is a sensor aligned coordinate system
19		* developed at JPL for radar mapping missions. Details pertaining to the
20		* coordinate system have been release in the public domain (see references
21		*above). This code is an independent implementation of the SCH coordinate
22		*system that conforms to the PROJ.4 conventions and uses the details presented
23		*in these publicly released documents. All credit for the development of the
24		*coordinate system and its use should be directed towards the original
25		*developers at JPL.
26		******************************************************************************
27		* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
28		* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29		* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
30		* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31		* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
32		* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
33		* DEALINGS IN THE SOFTWARE.
34		****************************************************************************/
35
36		#include <errno.h>
37		#include <math.h>
38
39		#include "proj.h"
40		#include "proj_internal.h"
41
42		namespace { // anonymous namespace
43		struct pj_sch_data {
44		double plat; /Peg Latitude /
45		double plon; /Peg Longitude/
46		double phdg; /Peg heading /
47		double h0; /Average altitude /
48		double transMat[9];
49		double xyzoff[3];
50		double rcurv;
51		PJ *cart;
52		PJ *cart_sph;
53		};
54		} // anonymous namespace
55
56		PROJ_HEAD(sch, "Spherical Cross-track Height")
57		"\n\tMisc\n\tplat_0= plon_0= phdg_0= [h_0=]";
58
59	0	static PJ_LPZ sch_inverse3d(PJ_XYZ xyz, PJ *P) {
60	0	struct pj_sch_data Q = static_cast<struct pj_sch_data >(P->opaque);
61
62	0	PJ_LPZ lpz;
63	0	lpz.lam = xyz.x * (P->a / Q->rcurv);
64	0	lpz.phi = xyz.y * (P->a / Q->rcurv);
65	0	lpz.z = xyz.z;
66	0	xyz = Q->cart_sph->fwd3d(lpz, Q->cart_sph);
67
68		/* Apply rotation */
69	0	xyz = {Q->transMat[0] * xyz.x + Q->transMat[1] * xyz.y +
70	0	Q->transMat[2] * xyz.z,
71	0	Q->transMat[3] * xyz.x + Q->transMat[4] * xyz.y +
72	0	Q->transMat[5] * xyz.z,
73	0	Q->transMat[6] * xyz.x + Q->transMat[7] * xyz.y +
74	0	Q->transMat[8] * xyz.z};
75
76		/* Apply offset */
77	0	xyz.x += Q->xyzoff[0];
78	0	xyz.y += Q->xyzoff[1];
79	0	xyz.z += Q->xyzoff[2];
80
81		/* Convert geocentric coordinates to lat long */
82	0	return Q->cart->inv3d(xyz, Q->cart);
83	0	}
84
85	0	static PJ_XYZ sch_forward3d(PJ_LPZ lpz, PJ *P) {
86	0	struct pj_sch_data Q = static_cast<struct pj_sch_data >(P->opaque);
87
88		/* Convert lat long to geocentric coordinates */
89	0	PJ_XYZ xyz = Q->cart->fwd3d(lpz, Q->cart);
90
91		/* Adjust for offset */
92	0	xyz.x -= Q->xyzoff[0];
93	0	xyz.y -= Q->xyzoff[1];
94	0	xyz.z -= Q->xyzoff[2];
95
96		/* Apply rotation */
97	0	xyz = {Q->transMat[0] * xyz.x + Q->transMat[3] * xyz.y +
98	0	Q->transMat[6] * xyz.z,
99	0	Q->transMat[1] * xyz.x + Q->transMat[4] * xyz.y +
100	0	Q->transMat[7] * xyz.z,
101	0	Q->transMat[2] * xyz.x + Q->transMat[5] * xyz.y +
102	0	Q->transMat[8] * xyz.z};
103
104		/* Convert to local lat,long */
105	0	lpz = Q->cart_sph->inv3d(xyz, Q->cart_sph);
106
107		/* Scale by radius */
108	0	xyz.x = lpz.lam * (Q->rcurv / P->a);
109	0	xyz.y = lpz.phi * (Q->rcurv / P->a);
110	0	xyz.z = lpz.z;
111
112	0	return xyz;
113	0	}
114
115	0	static PJ pj_sch_destructor(PJ P, int errlev) {
116	0	if (nullptr == P)
117	0	return nullptr;
118
119	0	auto Q = static_cast<struct pj_sch_data *>(P->opaque);
120	0	if (Q) {
121	0	if (Q->cart)
122	0	Q->cart->destructor(Q->cart, errlev);
123	0	if (Q->cart_sph)
124	0	Q->cart_sph->destructor(Q->cart_sph, errlev);
125	0	}
126
127	0	return pj_default_destructor(P, errlev);
128	0	}
129
130	0	static PJ pj_sch_setup(PJ P) { /* general initialization */
131	0	struct pj_sch_data Q = static_cast<struct pj_sch_data >(P->opaque);
132
133		/* Setup original geocentric system */
134		// Pass a dummy ellipsoid definition that will be overridden just afterwards
135	0	Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
136	0	if (Q->cart == nullptr)
137	0	return pj_sch_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
138
139		/* inherit ellipsoid definition from P to Q->cart */
140	0	pj_inherit_ellipsoid_def(P, Q->cart);
141
142	0	const double clt = cos(Q->plat);
143	0	const double slt = sin(Q->plat);
144	0	const double clo = cos(Q->plon);
145	0	const double slo = sin(Q->plon);
146
147		/* Estimate the radius of curvature for given peg */
148	0	const double temp = sqrt(1.0 - (P->es) * slt * slt);
149	0	const double reast = (P->a) / temp;
150	0	const double rnorth = (P->a) * (1.0 - (P->es)) / pow(temp, 3);
151
152	0	const double chdg = cos(Q->phdg);
153	0	const double shdg = sin(Q->phdg);
154
155	0	Q->rcurv =
156	0	Q->h0 + (reast * rnorth) / (reast * chdg * chdg + rnorth * shdg * shdg);
157
158		/* Set up local sphere at the given peg point */
159	0	Q->cart_sph = proj_create(P->ctx, "+proj=cart +a=1");
160	0	if (Q->cart_sph == nullptr)
161	0	return pj_sch_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
162	0	pj_calc_ellipsoid_params(Q->cart_sph, Q->rcurv, 0);
163
164		/* Set up the transformation matrices */
165	0	Q->transMat[0] = clt * clo;
166	0	Q->transMat[1] = -shdg * slo - slt * clo * chdg;
167	0	Q->transMat[2] = slo * chdg - slt * clo * shdg;
168	0	Q->transMat[3] = clt * slo;
169	0	Q->transMat[4] = clo * shdg - slt * slo * chdg;
170	0	Q->transMat[5] = -clo * chdg - slt * slo * shdg;
171	0	Q->transMat[6] = slt;
172	0	Q->transMat[7] = clt * chdg;
173	0	Q->transMat[8] = clt * shdg;
174
175	0	PJ_LPZ lpz;
176	0	lpz.lam = Q->plon;
177	0	lpz.phi = Q->plat;
178	0	lpz.z = Q->h0;
179	0	PJ_XYZ xyz = Q->cart->fwd3d(lpz, Q->cart);
180	0	Q->xyzoff[0] = xyz.x - (Q->rcurv) * clt * clo;
181	0	Q->xyzoff[1] = xyz.y - (Q->rcurv) * clt * slo;
182	0	Q->xyzoff[2] = xyz.z - (Q->rcurv) * slt;
183
184	0	P->fwd3d = sch_forward3d;
185	0	P->inv3d = sch_inverse3d;
186	0	return P;
187	0	}
188
189	4	PJ *PJ_PROJECTION(sch) {
190	4	struct pj_sch_data Q = static_cast<struct pj_sch_data >(
191	4	calloc(1, sizeof(struct pj_sch_data)));
192	4	if (nullptr == Q)
193	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
194	4	P->opaque = Q;
195	4	P->destructor = pj_sch_destructor;
196
197	4	Q->h0 = 0.0;
198
199		/* Check if peg latitude was defined */
200	4	if (pj_param(P->ctx, P->params, "tplat_0").i)
201	0	Q->plat = pj_param(P->ctx, P->params, "rplat_0").f;
202	4	else {
203	4	proj_log_error(P, _("Missing parameter plat_0."));
204	4	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
205	4	}
206
207		/* Check if peg longitude was defined */
208	0	if (pj_param(P->ctx, P->params, "tplon_0").i)
209	0	Q->plon = pj_param(P->ctx, P->params, "rplon_0").f;
210	0	else {
211	0	proj_log_error(P, _("Missing parameter plon_0."));
212	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
213	0	}
214
215		/* Check if peg heading is defined */
216	0	if (pj_param(P->ctx, P->params, "tphdg_0").i)
217	0	Q->phdg = pj_param(P->ctx, P->params, "rphdg_0").f;
218	0	else {
219	0	proj_log_error(P, _("Missing parameter phdg_0."));
220	0	return pj_default_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
221	0	}
222
223		/* Check if average height was defined - If so read it in */
224	0	if (pj_param(P->ctx, P->params, "th_0").i)
225	0	Q->h0 = pj_param(P->ctx, P->params, "dh_0").f;
226
227	0	return pj_sch_setup(P);
228	0	}