From df574ae332d57f556fd56314883b3354cab1d0ff Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 19 Dec 2018 13:00:37 +0100 Subject: cpp conversion: remove useless pj_, PJ_ and proj_ filename prefixes --- src/conversions/PJ_axisswap.cpp | 302 -------------------- src/conversions/PJ_cart.cpp | 219 --------------- src/conversions/PJ_geoc.cpp | 59 ---- src/conversions/PJ_unitconvert.cpp | 552 ------------------------------------- src/conversions/axisswap.cpp | 302 ++++++++++++++++++++ src/conversions/cart.cpp | 219 +++++++++++++++ src/conversions/geoc.cpp | 59 ++++ src/conversions/geocent.cpp | 62 +++++ src/conversions/pj_geocent.cpp | 62 ----- src/conversions/unitconvert.cpp | 552 +++++++++++++++++++++++++++++++++++++ 10 files changed, 1194 insertions(+), 1194 deletions(-) delete mode 100644 src/conversions/PJ_axisswap.cpp delete mode 100644 src/conversions/PJ_cart.cpp delete mode 100644 src/conversions/PJ_geoc.cpp delete mode 100644 src/conversions/PJ_unitconvert.cpp create mode 100644 src/conversions/axisswap.cpp create mode 100644 src/conversions/cart.cpp create mode 100644 src/conversions/geoc.cpp create mode 100644 src/conversions/geocent.cpp delete mode 100644 src/conversions/pj_geocent.cpp create mode 100644 src/conversions/unitconvert.cpp (limited to 'src/conversions') diff --git a/src/conversions/PJ_axisswap.cpp b/src/conversions/PJ_axisswap.cpp deleted file mode 100644 index 8714ec85..00000000 --- a/src/conversions/PJ_axisswap.cpp +++ /dev/null @@ -1,302 +0,0 @@ -/*********************************************************************** - - Axis order operation for use with transformation pipelines. - - Kristian Evers, kreve@sdfe.dk, 2017-10-31 - -************************************************************************ - -Change the order and sign of 2,3 or 4 axes. Each of the possible four -axes are numbered with 1-4, such that the first input axis is 1, the -second is 2 and so on. The output ordering is controlled by a list of the -input axes re-ordered to the new mapping. Examples: - -Reversing the order of the axes: - - +proj=axisswap +order=4,3,2,1 - -Swapping the first two axes (x and y): - - +proj=axisswap +order=2,1,3,4 - -The direction, or sign, of an axis can be changed by adding a minus in -front of the axis-number: - - +proj=axisswap +order=1,-2,3,4 - -It is only necessary to specify the axes that are affected by the swap -operation: - - +proj=axisswap +order=2,1 - -************************************************************************ -* Copyright (c) 2017, Kristian Evers / SDFE -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* to deal in the Software without restriction, including without limitation -* the rights to use, copy, modify, merge, publish, distribute, sublicense, -* and/or sell copies of the Software, and to permit persons to whom the -* Software is furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included -* in all copies or substantial portions of the Software. -* -* 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 -#include -#include - -#include "proj_internal.h" -#include "projects.h" - -PROJ_HEAD(axisswap, "Axis ordering"); - -namespace { // anonymous namespace -struct pj_opaque { - unsigned int axis[4]; - int sign[4]; -}; -} // anonymous namespace - -static int sign(int x) { - return (x > 0) - (x < 0); -} - -static XY forward_2d(LP lp, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - unsigned int i; - PJ_COORD out, in; - - in.lp = lp; - out = proj_coord_error(); - - for (i=0; i<2; i++) - out.v[i] = in.v[Q->axis[i]] * Q->sign[i]; - - return out.xy; -} - - -static LP reverse_2d(XY xy, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - unsigned int i; - PJ_COORD out, in; - - in.xy = xy; - out = proj_coord_error(); - - for (i=0; i<2; i++) - out.v[Q->axis[i]] = in.v[i] * Q->sign[i]; - - return out.lp; -} - - -static XYZ forward_3d(LPZ lpz, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - unsigned int i; - PJ_COORD out, in; - - in.lpz = lpz; - out = proj_coord_error(); - - for (i=0; i<3; i++) - out.v[i] = in.v[Q->axis[i]] * Q->sign[i]; - - return out.xyz; -} - -static LPZ reverse_3d(XYZ xyz, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - unsigned int i; - PJ_COORD in, out; - - out = proj_coord_error(); - in.xyz = xyz; - - for (i=0; i<3; i++) - out.v[Q->axis[i]] = in.v[i] * Q->sign[i]; - - return out.lpz; -} - - -static PJ_COORD forward_4d(PJ_COORD coo, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - unsigned int i; - PJ_COORD out; - - out = proj_coord_error(); - - for (i=0; i<4; i++) - out.v[i] = coo.v[Q->axis[i]] * Q->sign[i]; - - return out; -} - - -static PJ_COORD reverse_4d(PJ_COORD coo, PJ *P) { - struct pj_opaque *Q = (struct pj_opaque *) P->opaque; - unsigned int i; - PJ_COORD out; - - out = proj_coord_error(); - - for (i=0; i<4; i++) - out.v[Q->axis[i]] = coo.v[i] * Q->sign[i]; - - return out; -} - - -/***********************************************************************/ -PJ *CONVERSION(axisswap,0) { -/***********************************************************************/ - struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); - char *s; - unsigned int i, j, n = 0; - - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = (void *) Q; - - - /* +order and +axis are mutually exclusive */ - if ( !pj_param_exists(P->params, "order") == !pj_param_exists(P->params, "axis") ) - return pj_default_destructor(P, PJD_ERR_AXIS); - - /* fill axis list with indices from 4-7 to simplify duplicate search further down */ - for (i=0; i<4; i++) { - Q->axis[i] = i+4; - Q->sign[i] = 1; - } - - /* if the "order" parameter is used */ - if ( pj_param_exists(P->params, "order") ) { - /* read axis order */ - char *order = pj_param(P->ctx, P->params, "sorder").s; - - /* check that all characters are valid */ - for (i=0; iaxis[n] = abs(atoi(s))-1; - if (Q->axis[n] > 3) { - proj_log_error(P, "axisswap: invalid axis '%d'", Q->axis[n]); - return pj_default_destructor(P, PJD_ERR_AXIS); - } - Q->sign[n++] = sign(atoi(s)); - while ( *s != '\0' && *s != ',' ) - s++; - if ( *s == ',' ) - s++; - } - } - - /* if the "axis" parameter is used */ - if ( pj_param_exists(P->params, "axis") ) { - /* parse the classic PROJ.4 enu axis specification */ - for (i=0; i < 3; i++) { - switch(P->axis[i]) { - case 'w': - Q->sign[i] = -1; - Q->axis[i] = 0; - break; - case 'e': - Q->sign[i] = 1; - Q->axis[i] = 0; - break; - case 's': - Q->sign[i] = -1; - Q->axis[i] = 1; - break; - case 'n': - Q->sign[i] = 1; - Q->axis[i] = 1; - break; - case 'd': - Q->sign[i] = -1; - Q->axis[i] = 2; - break; - case 'u': - Q->sign[i] = 1; - Q->axis[i] = 2; - break; - default: - proj_log_error(P, "axisswap: unknown axis '%c'", P->axis[i]); - return pj_default_destructor(P, PJD_ERR_AXIS); - } - } - n = 3; - } - - /* check for duplicate axes */ - for (i=0; i<4; i++) - for (j=0; j<4; j++) { - if (i==j) - continue; - if (Q->axis[i] == Q->axis[j]) { - proj_log_error(P, "swapaxis: duplicate axes specified"); - return pj_default_destructor(P, PJD_ERR_AXIS); - } - } - - - /* only map fwd/inv functions that are possible with the given axis setup */ - if (n == 4) { - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - } - if (n == 3 && Q->axis[0] < 3 && Q->axis[1] < 3 && Q->axis[2] < 3) { - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - } - if (n == 2 && Q->axis[0] < 2 && Q->axis[1] < 2) { - P->fwd = forward_2d; - P->inv = reverse_2d; - } - - - if (P->fwd4d == nullptr && P->fwd3d == nullptr && P->fwd == nullptr) { - proj_log_error(P, "swapaxis: bad axis order"); - return pj_default_destructor(P, PJD_ERR_AXIS); - } - - if (pj_param(P->ctx, P->params, "tangularunits").i) { - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - } else { - P->left = PJ_IO_UNITS_WHATEVER; - P->right = PJ_IO_UNITS_WHATEVER; - } - - - /* Preparation and finalization steps are skipped, since the raison */ - /* d'etre of axisswap is to bring input coordinates in line with the */ - /* the internally expected order (ENU), such that handling of offsets */ - /* etc. can be done correctly in a later step of a pipeline */ - P->skip_fwd_prepare = 1; - P->skip_fwd_finalize = 1; - P->skip_inv_prepare = 1; - P->skip_inv_finalize = 1; - - return P; -} diff --git a/src/conversions/PJ_cart.cpp b/src/conversions/PJ_cart.cpp deleted file mode 100644 index 6fed9985..00000000 --- a/src/conversions/PJ_cart.cpp +++ /dev/null @@ -1,219 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Convert between ellipsoidal, geodetic coordinates and - * cartesian, geocentric coordinates. - * - * Formally, this functionality is also found in the PJ_geocent.c - * code. - * - * Actually, however, the PJ_geocent transformations are carried - * out in concert between 2D stubs in PJ_geocent.c and 3D code - * placed in pj_transform.c. - * - * For pipeline-style datum shifts, we do need direct access - * to the full 3D interface for this functionality. - * - * Hence this code, which may look like "just another PJ_geocent" - * but really is something substantially different. - * - * Author: Thomas Knudsen, thokn@sdfe.dk - * - ****************************************************************************** - * Copyright (c) 2016, Thomas Knudsen / SDFE - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * 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 "proj_internal.h" -#include "projects.h" -#include "proj_math.h" - -PROJ_HEAD(cart, "Geodetic/cartesian conversions"); - - -/************************************************************** - CARTESIAN / GEODETIC CONVERSIONS -*************************************************************** - This material follows: - - Bernhard Hofmann-Wellenhof & Helmut Moritz: - Physical Geodesy, 2nd edition. - Springer, 2005. - - chapter 5.6: Coordinate transformations - (HM, below), - - and - - Wikipedia: Geographic Coordinate Conversion, - https://en.wikipedia.org/wiki/Geographic_coordinate_conversion - - (WP, below). - - The cartesian-to-geodetic conversion is based on Bowring's - celebrated method: - - B. R. Bowring: - Transformation from spatial to geographical coordinates - Survey Review 23(181), pp. 323-327, 1976 - - (BB, below), - - but could probably use some TLC from a newer and faster - algorithm: - - Toshio Fukushima: - Transformation from Cartesian to Geodetic Coordinates - Accelerated by Halley’s Method - Journal of Geodesy, February 2006 - - (TF, below). - - Close to the poles, we avoid singularities by switching to an - approximation requiring knowledge of the geocentric radius - at the given latitude. For this, we use an adaptation of the - formula given in: - - Wikipedia: Earth Radius - https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude - (Derivation and commentary at https://gis.stackexchange.com/q/20200) - - (WP2, below) - - These routines are probably not as robust at those in - geocent.c, at least thay haven't been through as heavy - use as their geocent sisters. Some care has been taken - to avoid singularities, but extreme cases (e.g. setting - es, the squared eccentricity, to 1), will cause havoc. - -**************************************************************/ - - -/*********************************************************************/ -static double normal_radius_of_curvature (double a, double es, double phi) { -/*********************************************************************/ - double s = sin(phi); - if (es==0) - return a; - /* This is from WP. HM formula 2-149 gives an a,b version */ - return a / sqrt (1 - es*s*s); -} - -/*********************************************************************/ -static double geocentric_radius (double a, double b, double phi) { -/********************************************************************* - Return the geocentric radius at latitude phi, of an ellipsoid - with semimajor axis a and semiminor axis b. - - This is from WP2, but uses hypot() for potentially better - numerical robustness -***********************************************************************/ - return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi)); -} - - -/*********************************************************************/ -static XYZ cartesian (LPZ geod, PJ *P) { -/*********************************************************************/ - double N, cosphi = cos(geod.phi); - XYZ xyz; - - N = normal_radius_of_curvature(P->a, P->es, geod.phi); - - /* HM formula 5-27 (z formula follows WP) */ - xyz.x = (N + geod.z) * cosphi * cos(geod.lam); - xyz.y = (N + geod.z) * cosphi * sin(geod.lam); - xyz.z = (N * (1 - P->es) + geod.z) * sin(geod.phi); - - return xyz; -} - - -/*********************************************************************/ -static LPZ geodetic (XYZ cart, PJ *P) { -/*********************************************************************/ - double N, p, theta, c, s; - LPZ lpz; - - /* Perpendicular distance from point to Z-axis (HM eq. 5-28) */ - p = hypot (cart.x, cart.y); - - /* HM eq. (5-37) */ - theta = atan2 (cart.z * P->a, p * P->b); - - /* HM eq. (5-36) (from BB, 1976) */ - c = cos(theta); - s = sin(theta); - lpz.phi = atan2 (cart.z + P->e2s*P->b*s*s*s, p - P->es*P->a*c*c*c); - lpz.lam = atan2 (cart.y, cart.x); - N = normal_radius_of_curvature (P->a, P->es, lpz.phi); - - - c = cos(lpz.phi); - if (fabs(c) < 1e-6) { - /* poleward of 89.99994 deg, we avoid division by zero */ - /* by computing the height as the cartesian z value */ - /* minus the geocentric radius of the Earth at the given */ - /* latitude */ - double r = geocentric_radius (P->a, P->b, lpz.phi); - lpz.z = fabs (cart.z) - r; - } - else - lpz.z = p / c - N; - - return lpz; -} - - - -/* In effect, 2 cartesian coordinates of a point on the ellipsoid. Rather pointless, but... */ -static XY cart_forward (LP lp, PJ *P) { - PJ_COORD point; - point.lp = lp; - point.lpz.z = 0; - - point.xyz = cartesian (point.lpz, P); - return point.xy; -} - -/* And the other way round. Still rather pointless, but... */ -static LP cart_reverse (XY xy, PJ *P) { - PJ_COORD point; - point.xy = xy; - point.xyz.z = 0; - - point.lpz = geodetic (point.xyz, P); - return point.lp; -} - - - -/*********************************************************************/ -PJ *CONVERSION(cart,1) { -/*********************************************************************/ - P->fwd3d = cartesian; - P->inv3d = geodetic; - P->fwd = cart_forward; - P->inv = cart_reverse; - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_CARTESIAN; - return P; -} diff --git a/src/conversions/PJ_geoc.cpp b/src/conversions/PJ_geoc.cpp deleted file mode 100644 index 0455fada..00000000 --- a/src/conversions/PJ_geoc.cpp +++ /dev/null @@ -1,59 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Conversion from geographic to geocentric latitude and back. - * Author: Thomas Knudsen (2017) - * - ****************************************************************************** - * Copyright (c) 2017, SDFE, http://www.sdfe.dk - * Copyright (c) 2017, Thomas Knudsen - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * 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 - -#include "proj.h" -#include "proj_internal.h" -#include "projects.h" - -PROJ_HEAD(geoc, "Geocentric Latitude"); - -/* Geographical to geocentric */ -static PJ_COORD forward(PJ_COORD coo, PJ *P) { - return pj_geocentric_latitude (P, PJ_FWD, coo); -} - -/* Geocentric to geographical */ -static PJ_COORD inverse(PJ_COORD coo, PJ *P) { - return pj_geocentric_latitude (P, PJ_INV, coo); -} - - -static PJ *CONVERSION(geoc, 1) { - P->inv4d = inverse; - P->fwd4d = forward; - - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_ANGULAR; - - P->is_latlong = 1; - return P; -} diff --git a/src/conversions/PJ_unitconvert.cpp b/src/conversions/PJ_unitconvert.cpp deleted file mode 100644 index b25fd5d2..00000000 --- a/src/conversions/PJ_unitconvert.cpp +++ /dev/null @@ -1,552 +0,0 @@ -/*********************************************************************** - - Unit conversion pseudo-projection for use with - transformation pipelines. - - Kristian Evers, 2017-05-16 - -************************************************************************ - -A pseudo-projection that can be used to convert units of input and -output data. Primarily useful in pipelines. - -Unit conversion is performed by means of a pivot unit. The pivot unit -for distance units are the meter and for time we use the modified julian -date. A time unit conversion is performed like - - Unit A -> Modified Julian date -> Unit B - -distance units are converted in the same manner, with meter being the -central unit. - -The modified Julian date is chosen as the pivot unit since it has a -fairly high precision, goes sufficiently long backwards in time, has no -danger of hitting the upper limit in the near future and it is a fairly -common time unit in astronomy and geodesy. Note that we are using the -Julian date and not day. The difference being that the latter is defined -as an integer and is thus limited to days in resolution. This approach -has been extended wherever it makes sense, e.g. the GPS week unit also -has a fractional part that makes it possible to determine the day, hour -and minute of an observation. - -In- and output units are controlled with the parameters - - +xy_in, +xy_out, +z_in, +z_out, +t_in and +t_out - -where xy denotes horizontal units, z vertical units and t time units. - -************************************************************************ - -Kristian Evers, kreve@sdfe.dk, 2017-05-09 -Last update: 2017-05-16 - -************************************************************************ -* Copyright (c) 2017, Kristian Evers / SDFE -* -* Permission is hereby granted, free of charge, to any person obtaining a -* copy of this software and associated documentation files (the "Software"), -* to deal in the Software without restriction, including without limitation -* the rights to use, copy, modify, merge, publish, distribute, sublicense, -* and/or sell copies of the Software, and to permit persons to whom the -* Software is furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included -* in all copies or substantial portions of the Software. -* -* 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 -#include -#include -#include - -#include "proj_internal.h" -#include "proj_math.h" -#include "projects.h" - -PROJ_HEAD(unitconvert, "Unit conversion"); - -typedef double (*tconvert)(double); - -namespace { // anonymous namespace -struct TIME_UNITS { - const char *id; /* units keyword */ - tconvert t_in; /* unit -> mod. julian date function pointer */ - tconvert t_out; /* mod. julian date > unit function pointer */ - const char *name; /* comments */ -}; -} // anonymous namespace - -namespace { // anonymous namespace -struct pj_opaque_unitconvert { - int t_in_id; /* time unit id for the time input unit */ - int t_out_id; /* time unit id for the time output unit */ - double xy_factor; /* unit conversion factor for horizontal components */ - double z_factor; /* unit conversion factor for vertical components */ -}; -} // anonymous namespace - - -/***********************************************************************/ -static int is_leap_year(long year) { -/***********************************************************************/ - return ((year % 4 == 0 && year % 100 != 0) || year % 400 ==0); -} - - -/***********************************************************************/ -static int days_in_year(long year) { -/***********************************************************************/ - return is_leap_year(year) ? 366 : 365; -} - -/***********************************************************************/ -static unsigned int days_in_month(unsigned long year, unsigned long month) { -/***********************************************************************/ - const unsigned int month_table[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}; - unsigned int days; - - if (month > 12) month = 12; - if (month == 0) month = 1; - - days = month_table[month-1]; - if (is_leap_year(year) && month == 2) days++; - - return days; -} - - -/***********************************************************************/ -static int daynumber_in_year(unsigned long year, unsigned long month, unsigned long day) { -/***********************************************************************/ - unsigned int daynumber=0, i; - - if (month > 12) month = 12; - if (month == 0) month = 1; - if (day > days_in_month(year, month)) day = days_in_month(year, month); - - for (i = 1; i < month; i++) - daynumber += days_in_month(year, i); - - daynumber += day; - - return daynumber; - -} - -/***********************************************************************/ -static double mjd_to_mjd(double mjd) { -/*********************************************************************** - Modified julian date no-op function. - - The Julian date is defined as (fractional) days since midnight - on 16th of November in 1858. -************************************************************************/ - return mjd; -} - - -/***********************************************************************/ -static double decimalyear_to_mjd(double decimalyear) { -/*********************************************************************** - Epoch of modified julian date is 1858-11-16 00:00 -************************************************************************/ - long year; - double fractional_year; - double mjd; - - if( decimalyear < -10000 || decimalyear > 10000 ) - return 0; - - year = lround(floor(decimalyear)); - fractional_year = decimalyear - year; - mjd = (year - 1859)*365 + 14 + 31; - mjd += (double)fractional_year*(double)days_in_year(year); - - /* take care of leap days */ - year--; - for (; year > 1858; year--) - if (is_leap_year(year)) - mjd++; - - return mjd; -} - - -/***********************************************************************/ -static double mjd_to_decimalyear(double mjd) { -/*********************************************************************** - Epoch of modified julian date is 1858-11-16 00:00 -************************************************************************/ - double decimalyear = mjd; - double mjd_iter = 14 + 31; - int year = 1859; - - /* a smarter brain than mine could probably to do this more elegantly - - I'll just brute-force my way out of this... */ - for (; mjd >= mjd_iter; year++) { - mjd_iter += days_in_year(year); - } - year--; - mjd_iter -= days_in_year(year); - - decimalyear = year + (mjd-mjd_iter)/days_in_year(year); - return decimalyear; -} - - -/***********************************************************************/ -static double gps_week_to_mjd(double gps_week) { -/*********************************************************************** - GPS weeks are defined as the number of weeks since January the 6th - 1980. - - Epoch of gps weeks is 1980-01-06 00:00, which in modified Julian - date is 44244. -************************************************************************/ - return 44244.0 + gps_week*7.0; -} - - -/***********************************************************************/ -static double mjd_to_gps_week(double mjd) { -/*********************************************************************** - GPS weeks are defined as the number of weeks since January the 6th - 1980. - - Epoch of gps weeks is 1980-01-06 00:00, which in modified Julian - date is 44244. -************************************************************************/ - return (mjd - 44244.0) / 7.0; -} - - -/***********************************************************************/ -static double yyyymmdd_to_mjd(double yyyymmdd) { -/************************************************************************ - Date given in YYYY-MM-DD format. -************************************************************************/ - - long year = lround(floor( yyyymmdd / 10000 )); - long month = lround(floor((yyyymmdd - year*10000) / 100)); - long day = lround(floor( yyyymmdd - year*10000 - month*100)); - double mjd = daynumber_in_year(year, month, day); - - for (year -= 1; year > 1858; year--) - mjd += days_in_year(year); - - return mjd + 13 + 31; -} - - -/***********************************************************************/ -static double mjd_to_yyyymmdd(double mjd) { -/************************************************************************ - Date given in YYYY-MM-DD format. -************************************************************************/ - double mjd_iter = 14 + 31; - int year = 1859, month=0, day=0; - - for (; mjd >= mjd_iter; year++) { - mjd_iter += days_in_year(year); - } - year--; - mjd_iter -= days_in_year(year); - - for (month=1; mjd_iter + days_in_month(year, month) <= mjd; month++) - mjd_iter += days_in_month(year, month); - - day = (int)(mjd - mjd_iter + 1); - - return year*10000.0 + month*100.0 + day; -} - -static const struct TIME_UNITS time_units[] = { - {"mjd", mjd_to_mjd, mjd_to_mjd, "Modified julian date"}, - {"decimalyear", decimalyear_to_mjd, mjd_to_decimalyear, "Decimal year"}, - {"gps_week", gps_week_to_mjd, mjd_to_gps_week, "GPS Week"}, - {"yyyymmdd", yyyymmdd_to_mjd, mjd_to_yyyymmdd, "YYYYMMDD date"}, - {nullptr, nullptr, nullptr, nullptr} -}; - - -/***********************************************************************/ -static XY forward_2d(LP lp, PJ *P) { -/************************************************************************ - Forward unit conversions in the plane -************************************************************************/ - struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - point.lp = lp; - - point.xy.x *= Q->xy_factor; - point.xy.y *= Q->xy_factor; - - return point.xy; -} - - -/***********************************************************************/ -static LP reverse_2d(XY xy, PJ *P) { -/************************************************************************ - Reverse unit conversions in the plane -************************************************************************/ - struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - point.xy = xy; - - point.xy.x /= Q->xy_factor; - point.xy.y /= Q->xy_factor; - - return point.lp; -} - - -/***********************************************************************/ -static XYZ forward_3d(LPZ lpz, PJ *P) { -/************************************************************************ - Forward unit conversions the vertical component -************************************************************************/ - struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - point.lpz = lpz; - - /* take care of the horizontal components in the 2D function */ - point.xy = forward_2d(point.lp, P); - - point.xyz.z *= Q->z_factor; - - return point.xyz; -} - -/***********************************************************************/ -static LPZ reverse_3d(XYZ xyz, PJ *P) { -/************************************************************************ - Reverse unit conversions the vertical component -************************************************************************/ - struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; - PJ_COORD point = {{0,0,0,0}}; - point.xyz = xyz; - - /* take care of the horizontal components in the 2D function */ - point.lp = reverse_2d(point.xy, P); - - point.xyz.z /= Q->z_factor; - - return point.lpz; -} - - -/***********************************************************************/ -static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { -/************************************************************************ - Forward conversion of time units -************************************************************************/ - struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; - PJ_COORD out = obs; - - /* delegate unit conversion of physical dimensions to the 3D function */ - out.xyz = forward_3d(obs.lpz, P); - - if (Q->t_in_id >= 0) - out.xyzt.t = time_units[Q->t_in_id].t_in( obs.xyzt.t ); - if (Q->t_out_id >= 0) - out.xyzt.t = time_units[Q->t_out_id].t_out( out.xyzt.t ); - - return out; -} - - -/***********************************************************************/ -static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { -/************************************************************************ - Reverse conversion of time units -************************************************************************/ - struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; - PJ_COORD out = obs; - - /* delegate unit conversion of physical dimensions to the 3D function */ - out.lpz = reverse_3d(obs.xyz, P); - - if (Q->t_out_id >= 0) - out.xyzt.t = time_units[Q->t_out_id].t_in( obs.xyzt.t ); - if (Q->t_in_id >= 0) - out.xyzt.t = time_units[Q->t_in_id].t_out( out.xyzt.t ); - - return out; -} - -/***********************************************************************/ -static double get_unit_conversion_factor(const char* name, - int* p_is_linear, - const char** p_normalized_name) { -/***********************************************************************/ - int i; - const char* s; - const PJ_UNITS *units; - - units = proj_list_units(); - - /* Try first with linear units */ - for (i = 0; (s = units[i].id) ; ++i) { - if ( strcmp(s, name) == 0 ) { - if( p_normalized_name ) { - *p_normalized_name = units[i].name; - } - if( p_is_linear ) { - *p_is_linear = 1; - } - return units[i].factor; - } - } - - /* And then angular units */ - units = proj_list_angular_units(); - for (i = 0; (s = units[i].id) ; ++i) { - if ( strcmp(s, name) == 0 ) { - if( p_normalized_name ) { - *p_normalized_name = units[i].name; - } - if( p_is_linear ) { - *p_is_linear = 0; - } - return units[i].factor; - } - } - if( p_normalized_name ) { - *p_normalized_name = nullptr; - } - if( p_is_linear ) { - *p_is_linear = -1; - } - return 0.0; -} - -/***********************************************************************/ -PJ *CONVERSION(unitconvert,0) { -/***********************************************************************/ - struct pj_opaque_unitconvert *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_unitconvert))); - const char *s, *name; - int i; - double f; - int xy_in_is_linear = -1; /* unknown */ - int xy_out_is_linear = -1; /* unknown */ - int z_in_is_linear = -1; /* unknown */ - int z_out_is_linear = -1; /* unknown */ - - if (nullptr==Q) - return pj_default_destructor (P, ENOMEM); - P->opaque = (void *) Q; - - P->fwd4d = forward_4d; - P->inv4d = reverse_4d; - P->fwd3d = forward_3d; - P->inv3d = reverse_3d; - P->fwd = forward_2d; - P->inv = reverse_2d; - - P->left = PJ_IO_UNITS_WHATEVER; - P->right = PJ_IO_UNITS_WHATEVER; - - /* if no time input/output unit is specified we can skip them */ - Q->t_in_id = -1; - Q->t_out_id = -1; - - Q->xy_factor = 1.0; - Q->z_factor = 1.0; - - if ((name = pj_param (P->ctx, P->params, "sxy_in").s) != nullptr) { - const char* normalized_name = nullptr; - f = get_unit_conversion_factor(name, &xy_in_is_linear, &normalized_name); - if (f != 0.0) { - proj_log_debug(P, "xy_in unit: %s", normalized_name); - } else { - if ( (f = pj_param (P->ctx, P->params, "dxy_in").f) == 0.0) - return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); - } - if (f != 0.0) - Q->xy_factor *= f; - } - - if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != nullptr) { - const char* normalized_name = nullptr; - f = get_unit_conversion_factor(name, &xy_out_is_linear, &normalized_name); - if (f != 0.0) { - proj_log_debug(P, "xy_out unit: %s", normalized_name); - } else { - if ( (f = pj_param (P->ctx, P->params, "dxy_out").f) == 0.0) - return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); - } - if (f != 0.0) - Q->xy_factor /= f; - } - - if( xy_in_is_linear >= 0 && xy_out_is_linear >= 0 && - xy_in_is_linear != xy_out_is_linear ) { - proj_log_debug(P, "inconsistent unit type between xy_in and xy_out"); - return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT); - } - - if ((name = pj_param (P->ctx, P->params, "sz_in").s) != nullptr) { - const char* normalized_name = nullptr; - f = get_unit_conversion_factor(name, &z_in_is_linear, &normalized_name); - if (f != 0.0) { - proj_log_debug(P, "z_in unit: %s", normalized_name); - } else { - if ( (f = pj_param (P->ctx, P->params, "dz_in").f) == 0.0) - return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); - } - if (f != 0.0) - Q->z_factor *= f; - } - - if ((name = pj_param (P->ctx, P->params, "sz_out").s) != nullptr) { - const char* normalized_name = nullptr; - f = get_unit_conversion_factor(name, &z_out_is_linear, &normalized_name); - if (f != 0.0) { - proj_log_debug(P, "z_out unit: %s", normalized_name); - } else { - if ( (f = pj_param (P->ctx, P->params, "dz_out").f) == 0.0) - return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); - } - if (f != 0.0) - Q->z_factor /= f; - } - - if( z_in_is_linear >= 0 && z_out_is_linear >= 0 && - z_in_is_linear != z_out_is_linear ) { - proj_log_debug(P, "inconsistent unit type between z_in and z_out"); - return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT); - } - - if ((name = pj_param (P->ctx, P->params, "st_in").s) != nullptr) { - for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i); - - if (!s) return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); /* unknown unit conversion id */ - - Q->t_in_id = i; - proj_log_debug(P, "t_in unit: %s", time_units[i].name); - } - - s = nullptr; - if ((name = pj_param (P->ctx, P->params, "st_out").s) != nullptr) { - for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i); - - if (!s) return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); /* unknown unit conversion id */ - - Q->t_out_id = i; - proj_log_debug(P, "t_out unit: %s", time_units[i].name); - } - - return P; -} diff --git a/src/conversions/axisswap.cpp b/src/conversions/axisswap.cpp new file mode 100644 index 00000000..8714ec85 --- /dev/null +++ b/src/conversions/axisswap.cpp @@ -0,0 +1,302 @@ +/*********************************************************************** + + Axis order operation for use with transformation pipelines. + + Kristian Evers, kreve@sdfe.dk, 2017-10-31 + +************************************************************************ + +Change the order and sign of 2,3 or 4 axes. Each of the possible four +axes are numbered with 1-4, such that the first input axis is 1, the +second is 2 and so on. The output ordering is controlled by a list of the +input axes re-ordered to the new mapping. Examples: + +Reversing the order of the axes: + + +proj=axisswap +order=4,3,2,1 + +Swapping the first two axes (x and y): + + +proj=axisswap +order=2,1,3,4 + +The direction, or sign, of an axis can be changed by adding a minus in +front of the axis-number: + + +proj=axisswap +order=1,-2,3,4 + +It is only necessary to specify the axes that are affected by the swap +operation: + + +proj=axisswap +order=2,1 + +************************************************************************ +* Copyright (c) 2017, Kristian Evers / SDFE +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* 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 +#include +#include + +#include "proj_internal.h" +#include "projects.h" + +PROJ_HEAD(axisswap, "Axis ordering"); + +namespace { // anonymous namespace +struct pj_opaque { + unsigned int axis[4]; + int sign[4]; +}; +} // anonymous namespace + +static int sign(int x) { + return (x > 0) - (x < 0); +} + +static XY forward_2d(LP lp, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + unsigned int i; + PJ_COORD out, in; + + in.lp = lp; + out = proj_coord_error(); + + for (i=0; i<2; i++) + out.v[i] = in.v[Q->axis[i]] * Q->sign[i]; + + return out.xy; +} + + +static LP reverse_2d(XY xy, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + unsigned int i; + PJ_COORD out, in; + + in.xy = xy; + out = proj_coord_error(); + + for (i=0; i<2; i++) + out.v[Q->axis[i]] = in.v[i] * Q->sign[i]; + + return out.lp; +} + + +static XYZ forward_3d(LPZ lpz, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + unsigned int i; + PJ_COORD out, in; + + in.lpz = lpz; + out = proj_coord_error(); + + for (i=0; i<3; i++) + out.v[i] = in.v[Q->axis[i]] * Q->sign[i]; + + return out.xyz; +} + +static LPZ reverse_3d(XYZ xyz, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + unsigned int i; + PJ_COORD in, out; + + out = proj_coord_error(); + in.xyz = xyz; + + for (i=0; i<3; i++) + out.v[Q->axis[i]] = in.v[i] * Q->sign[i]; + + return out.lpz; +} + + +static PJ_COORD forward_4d(PJ_COORD coo, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + unsigned int i; + PJ_COORD out; + + out = proj_coord_error(); + + for (i=0; i<4; i++) + out.v[i] = coo.v[Q->axis[i]] * Q->sign[i]; + + return out; +} + + +static PJ_COORD reverse_4d(PJ_COORD coo, PJ *P) { + struct pj_opaque *Q = (struct pj_opaque *) P->opaque; + unsigned int i; + PJ_COORD out; + + out = proj_coord_error(); + + for (i=0; i<4; i++) + out.v[Q->axis[i]] = coo.v[i] * Q->sign[i]; + + return out; +} + + +/***********************************************************************/ +PJ *CONVERSION(axisswap,0) { +/***********************************************************************/ + struct pj_opaque *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque))); + char *s; + unsigned int i, j, n = 0; + + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = (void *) Q; + + + /* +order and +axis are mutually exclusive */ + if ( !pj_param_exists(P->params, "order") == !pj_param_exists(P->params, "axis") ) + return pj_default_destructor(P, PJD_ERR_AXIS); + + /* fill axis list with indices from 4-7 to simplify duplicate search further down */ + for (i=0; i<4; i++) { + Q->axis[i] = i+4; + Q->sign[i] = 1; + } + + /* if the "order" parameter is used */ + if ( pj_param_exists(P->params, "order") ) { + /* read axis order */ + char *order = pj_param(P->ctx, P->params, "sorder").s; + + /* check that all characters are valid */ + for (i=0; iaxis[n] = abs(atoi(s))-1; + if (Q->axis[n] > 3) { + proj_log_error(P, "axisswap: invalid axis '%d'", Q->axis[n]); + return pj_default_destructor(P, PJD_ERR_AXIS); + } + Q->sign[n++] = sign(atoi(s)); + while ( *s != '\0' && *s != ',' ) + s++; + if ( *s == ',' ) + s++; + } + } + + /* if the "axis" parameter is used */ + if ( pj_param_exists(P->params, "axis") ) { + /* parse the classic PROJ.4 enu axis specification */ + for (i=0; i < 3; i++) { + switch(P->axis[i]) { + case 'w': + Q->sign[i] = -1; + Q->axis[i] = 0; + break; + case 'e': + Q->sign[i] = 1; + Q->axis[i] = 0; + break; + case 's': + Q->sign[i] = -1; + Q->axis[i] = 1; + break; + case 'n': + Q->sign[i] = 1; + Q->axis[i] = 1; + break; + case 'd': + Q->sign[i] = -1; + Q->axis[i] = 2; + break; + case 'u': + Q->sign[i] = 1; + Q->axis[i] = 2; + break; + default: + proj_log_error(P, "axisswap: unknown axis '%c'", P->axis[i]); + return pj_default_destructor(P, PJD_ERR_AXIS); + } + } + n = 3; + } + + /* check for duplicate axes */ + for (i=0; i<4; i++) + for (j=0; j<4; j++) { + if (i==j) + continue; + if (Q->axis[i] == Q->axis[j]) { + proj_log_error(P, "swapaxis: duplicate axes specified"); + return pj_default_destructor(P, PJD_ERR_AXIS); + } + } + + + /* only map fwd/inv functions that are possible with the given axis setup */ + if (n == 4) { + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + } + if (n == 3 && Q->axis[0] < 3 && Q->axis[1] < 3 && Q->axis[2] < 3) { + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + } + if (n == 2 && Q->axis[0] < 2 && Q->axis[1] < 2) { + P->fwd = forward_2d; + P->inv = reverse_2d; + } + + + if (P->fwd4d == nullptr && P->fwd3d == nullptr && P->fwd == nullptr) { + proj_log_error(P, "swapaxis: bad axis order"); + return pj_default_destructor(P, PJD_ERR_AXIS); + } + + if (pj_param(P->ctx, P->params, "tangularunits").i) { + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + } else { + P->left = PJ_IO_UNITS_WHATEVER; + P->right = PJ_IO_UNITS_WHATEVER; + } + + + /* Preparation and finalization steps are skipped, since the raison */ + /* d'etre of axisswap is to bring input coordinates in line with the */ + /* the internally expected order (ENU), such that handling of offsets */ + /* etc. can be done correctly in a later step of a pipeline */ + P->skip_fwd_prepare = 1; + P->skip_fwd_finalize = 1; + P->skip_inv_prepare = 1; + P->skip_inv_finalize = 1; + + return P; +} diff --git a/src/conversions/cart.cpp b/src/conversions/cart.cpp new file mode 100644 index 00000000..6fed9985 --- /dev/null +++ b/src/conversions/cart.cpp @@ -0,0 +1,219 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Convert between ellipsoidal, geodetic coordinates and + * cartesian, geocentric coordinates. + * + * Formally, this functionality is also found in the PJ_geocent.c + * code. + * + * Actually, however, the PJ_geocent transformations are carried + * out in concert between 2D stubs in PJ_geocent.c and 3D code + * placed in pj_transform.c. + * + * For pipeline-style datum shifts, we do need direct access + * to the full 3D interface for this functionality. + * + * Hence this code, which may look like "just another PJ_geocent" + * but really is something substantially different. + * + * Author: Thomas Knudsen, thokn@sdfe.dk + * + ****************************************************************************** + * Copyright (c) 2016, Thomas Knudsen / SDFE + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * 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 "proj_internal.h" +#include "projects.h" +#include "proj_math.h" + +PROJ_HEAD(cart, "Geodetic/cartesian conversions"); + + +/************************************************************** + CARTESIAN / GEODETIC CONVERSIONS +*************************************************************** + This material follows: + + Bernhard Hofmann-Wellenhof & Helmut Moritz: + Physical Geodesy, 2nd edition. + Springer, 2005. + + chapter 5.6: Coordinate transformations + (HM, below), + + and + + Wikipedia: Geographic Coordinate Conversion, + https://en.wikipedia.org/wiki/Geographic_coordinate_conversion + + (WP, below). + + The cartesian-to-geodetic conversion is based on Bowring's + celebrated method: + + B. R. Bowring: + Transformation from spatial to geographical coordinates + Survey Review 23(181), pp. 323-327, 1976 + + (BB, below), + + but could probably use some TLC from a newer and faster + algorithm: + + Toshio Fukushima: + Transformation from Cartesian to Geodetic Coordinates + Accelerated by Halley’s Method + Journal of Geodesy, February 2006 + + (TF, below). + + Close to the poles, we avoid singularities by switching to an + approximation requiring knowledge of the geocentric radius + at the given latitude. For this, we use an adaptation of the + formula given in: + + Wikipedia: Earth Radius + https://en.wikipedia.org/wiki/Earth_radius#Radius_at_a_given_geodetic_latitude + (Derivation and commentary at https://gis.stackexchange.com/q/20200) + + (WP2, below) + + These routines are probably not as robust at those in + geocent.c, at least thay haven't been through as heavy + use as their geocent sisters. Some care has been taken + to avoid singularities, but extreme cases (e.g. setting + es, the squared eccentricity, to 1), will cause havoc. + +**************************************************************/ + + +/*********************************************************************/ +static double normal_radius_of_curvature (double a, double es, double phi) { +/*********************************************************************/ + double s = sin(phi); + if (es==0) + return a; + /* This is from WP. HM formula 2-149 gives an a,b version */ + return a / sqrt (1 - es*s*s); +} + +/*********************************************************************/ +static double geocentric_radius (double a, double b, double phi) { +/********************************************************************* + Return the geocentric radius at latitude phi, of an ellipsoid + with semimajor axis a and semiminor axis b. + + This is from WP2, but uses hypot() for potentially better + numerical robustness +***********************************************************************/ + return hypot (a*a*cos (phi), b*b*sin(phi)) / hypot (a*cos(phi), b*sin(phi)); +} + + +/*********************************************************************/ +static XYZ cartesian (LPZ geod, PJ *P) { +/*********************************************************************/ + double N, cosphi = cos(geod.phi); + XYZ xyz; + + N = normal_radius_of_curvature(P->a, P->es, geod.phi); + + /* HM formula 5-27 (z formula follows WP) */ + xyz.x = (N + geod.z) * cosphi * cos(geod.lam); + xyz.y = (N + geod.z) * cosphi * sin(geod.lam); + xyz.z = (N * (1 - P->es) + geod.z) * sin(geod.phi); + + return xyz; +} + + +/*********************************************************************/ +static LPZ geodetic (XYZ cart, PJ *P) { +/*********************************************************************/ + double N, p, theta, c, s; + LPZ lpz; + + /* Perpendicular distance from point to Z-axis (HM eq. 5-28) */ + p = hypot (cart.x, cart.y); + + /* HM eq. (5-37) */ + theta = atan2 (cart.z * P->a, p * P->b); + + /* HM eq. (5-36) (from BB, 1976) */ + c = cos(theta); + s = sin(theta); + lpz.phi = atan2 (cart.z + P->e2s*P->b*s*s*s, p - P->es*P->a*c*c*c); + lpz.lam = atan2 (cart.y, cart.x); + N = normal_radius_of_curvature (P->a, P->es, lpz.phi); + + + c = cos(lpz.phi); + if (fabs(c) < 1e-6) { + /* poleward of 89.99994 deg, we avoid division by zero */ + /* by computing the height as the cartesian z value */ + /* minus the geocentric radius of the Earth at the given */ + /* latitude */ + double r = geocentric_radius (P->a, P->b, lpz.phi); + lpz.z = fabs (cart.z) - r; + } + else + lpz.z = p / c - N; + + return lpz; +} + + + +/* In effect, 2 cartesian coordinates of a point on the ellipsoid. Rather pointless, but... */ +static XY cart_forward (LP lp, PJ *P) { + PJ_COORD point; + point.lp = lp; + point.lpz.z = 0; + + point.xyz = cartesian (point.lpz, P); + return point.xy; +} + +/* And the other way round. Still rather pointless, but... */ +static LP cart_reverse (XY xy, PJ *P) { + PJ_COORD point; + point.xy = xy; + point.xyz.z = 0; + + point.lpz = geodetic (point.xyz, P); + return point.lp; +} + + + +/*********************************************************************/ +PJ *CONVERSION(cart,1) { +/*********************************************************************/ + P->fwd3d = cartesian; + P->inv3d = geodetic; + P->fwd = cart_forward; + P->inv = cart_reverse; + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_CARTESIAN; + return P; +} diff --git a/src/conversions/geoc.cpp b/src/conversions/geoc.cpp new file mode 100644 index 00000000..0455fada --- /dev/null +++ b/src/conversions/geoc.cpp @@ -0,0 +1,59 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Conversion from geographic to geocentric latitude and back. + * Author: Thomas Knudsen (2017) + * + ****************************************************************************** + * Copyright (c) 2017, SDFE, http://www.sdfe.dk + * Copyright (c) 2017, Thomas Knudsen + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * 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 + +#include "proj.h" +#include "proj_internal.h" +#include "projects.h" + +PROJ_HEAD(geoc, "Geocentric Latitude"); + +/* Geographical to geocentric */ +static PJ_COORD forward(PJ_COORD coo, PJ *P) { + return pj_geocentric_latitude (P, PJ_FWD, coo); +} + +/* Geocentric to geographical */ +static PJ_COORD inverse(PJ_COORD coo, PJ *P) { + return pj_geocentric_latitude (P, PJ_INV, coo); +} + + +static PJ *CONVERSION(geoc, 1) { + P->inv4d = inverse; + P->fwd4d = forward; + + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_ANGULAR; + + P->is_latlong = 1; + return P; +} diff --git a/src/conversions/geocent.cpp b/src/conversions/geocent.cpp new file mode 100644 index 00000000..0e9d725e --- /dev/null +++ b/src/conversions/geocent.cpp @@ -0,0 +1,62 @@ +/****************************************************************************** + * Project: PROJ.4 + * Purpose: Stub projection for geocentric. The transformation isn't + * really done here since this code is 2D. The real transformation + * is handled by pj_transform.c. + * Author: Frank Warmerdam, warmerdam@pobox.com + * + ****************************************************************************** + * Copyright (c) 2002, Frank Warmerdam + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * 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 "projects.h" + +PROJ_HEAD(geocent, "Geocentric") "\n\t"; + +static XY forward(LP lp, PJ *P) { + XY xy = {0.0,0.0}; + (void) P; + xy.x = lp.lam; + xy.y = lp.phi; + return xy; +} + +static LP inverse(XY xy, PJ *P) { + LP lp = {0.0,0.0}; + (void) P; + lp.phi = xy.y; + lp.lam = xy.x; + return lp; +} + +PJ *CONVERSION (geocent, 0) { + P->is_geocent = 1; + P->x0 = 0.0; + P->y0 = 0.0; + P->inv = inverse; + P->fwd = forward; + P->left = PJ_IO_UNITS_ANGULAR; + P->right = PJ_IO_UNITS_CARTESIAN; + + return P; +} + diff --git a/src/conversions/pj_geocent.cpp b/src/conversions/pj_geocent.cpp deleted file mode 100644 index 0e9d725e..00000000 --- a/src/conversions/pj_geocent.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Stub projection for geocentric. The transformation isn't - * really done here since this code is 2D. The real transformation - * is handled by pj_transform.c. - * Author: Frank Warmerdam, warmerdam@pobox.com - * - ****************************************************************************** - * Copyright (c) 2002, Frank Warmerdam - * - * Permission is hereby granted, free of charge, to any person obtaining a - * copy of this software and associated documentation files (the "Software"), - * to deal in the Software without restriction, including without limitation - * the rights to use, copy, modify, merge, publish, distribute, sublicense, - * and/or sell copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * 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 "projects.h" - -PROJ_HEAD(geocent, "Geocentric") "\n\t"; - -static XY forward(LP lp, PJ *P) { - XY xy = {0.0,0.0}; - (void) P; - xy.x = lp.lam; - xy.y = lp.phi; - return xy; -} - -static LP inverse(XY xy, PJ *P) { - LP lp = {0.0,0.0}; - (void) P; - lp.phi = xy.y; - lp.lam = xy.x; - return lp; -} - -PJ *CONVERSION (geocent, 0) { - P->is_geocent = 1; - P->x0 = 0.0; - P->y0 = 0.0; - P->inv = inverse; - P->fwd = forward; - P->left = PJ_IO_UNITS_ANGULAR; - P->right = PJ_IO_UNITS_CARTESIAN; - - return P; -} - diff --git a/src/conversions/unitconvert.cpp b/src/conversions/unitconvert.cpp new file mode 100644 index 00000000..b25fd5d2 --- /dev/null +++ b/src/conversions/unitconvert.cpp @@ -0,0 +1,552 @@ +/*********************************************************************** + + Unit conversion pseudo-projection for use with + transformation pipelines. + + Kristian Evers, 2017-05-16 + +************************************************************************ + +A pseudo-projection that can be used to convert units of input and +output data. Primarily useful in pipelines. + +Unit conversion is performed by means of a pivot unit. The pivot unit +for distance units are the meter and for time we use the modified julian +date. A time unit conversion is performed like + + Unit A -> Modified Julian date -> Unit B + +distance units are converted in the same manner, with meter being the +central unit. + +The modified Julian date is chosen as the pivot unit since it has a +fairly high precision, goes sufficiently long backwards in time, has no +danger of hitting the upper limit in the near future and it is a fairly +common time unit in astronomy and geodesy. Note that we are using the +Julian date and not day. The difference being that the latter is defined +as an integer and is thus limited to days in resolution. This approach +has been extended wherever it makes sense, e.g. the GPS week unit also +has a fractional part that makes it possible to determine the day, hour +and minute of an observation. + +In- and output units are controlled with the parameters + + +xy_in, +xy_out, +z_in, +z_out, +t_in and +t_out + +where xy denotes horizontal units, z vertical units and t time units. + +************************************************************************ + +Kristian Evers, kreve@sdfe.dk, 2017-05-09 +Last update: 2017-05-16 + +************************************************************************ +* Copyright (c) 2017, Kristian Evers / SDFE +* +* Permission is hereby granted, free of charge, to any person obtaining a +* copy of this software and associated documentation files (the "Software"), +* to deal in the Software without restriction, including without limitation +* the rights to use, copy, modify, merge, publish, distribute, sublicense, +* and/or sell copies of the Software, and to permit persons to whom the +* Software is furnished to do so, subject to the following conditions: +* +* The above copyright notice and this permission notice shall be included +* in all copies or substantial portions of the Software. +* +* 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 +#include +#include +#include + +#include "proj_internal.h" +#include "proj_math.h" +#include "projects.h" + +PROJ_HEAD(unitconvert, "Unit conversion"); + +typedef double (*tconvert)(double); + +namespace { // anonymous namespace +struct TIME_UNITS { + const char *id; /* units keyword */ + tconvert t_in; /* unit -> mod. julian date function pointer */ + tconvert t_out; /* mod. julian date > unit function pointer */ + const char *name; /* comments */ +}; +} // anonymous namespace + +namespace { // anonymous namespace +struct pj_opaque_unitconvert { + int t_in_id; /* time unit id for the time input unit */ + int t_out_id; /* time unit id for the time output unit */ + double xy_factor; /* unit conversion factor for horizontal components */ + double z_factor; /* unit conversion factor for vertical components */ +}; +} // anonymous namespace + + +/***********************************************************************/ +static int is_leap_year(long year) { +/***********************************************************************/ + return ((year % 4 == 0 && year % 100 != 0) || year % 400 ==0); +} + + +/***********************************************************************/ +static int days_in_year(long year) { +/***********************************************************************/ + return is_leap_year(year) ? 366 : 365; +} + +/***********************************************************************/ +static unsigned int days_in_month(unsigned long year, unsigned long month) { +/***********************************************************************/ + const unsigned int month_table[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}; + unsigned int days; + + if (month > 12) month = 12; + if (month == 0) month = 1; + + days = month_table[month-1]; + if (is_leap_year(year) && month == 2) days++; + + return days; +} + + +/***********************************************************************/ +static int daynumber_in_year(unsigned long year, unsigned long month, unsigned long day) { +/***********************************************************************/ + unsigned int daynumber=0, i; + + if (month > 12) month = 12; + if (month == 0) month = 1; + if (day > days_in_month(year, month)) day = days_in_month(year, month); + + for (i = 1; i < month; i++) + daynumber += days_in_month(year, i); + + daynumber += day; + + return daynumber; + +} + +/***********************************************************************/ +static double mjd_to_mjd(double mjd) { +/*********************************************************************** + Modified julian date no-op function. + + The Julian date is defined as (fractional) days since midnight + on 16th of November in 1858. +************************************************************************/ + return mjd; +} + + +/***********************************************************************/ +static double decimalyear_to_mjd(double decimalyear) { +/*********************************************************************** + Epoch of modified julian date is 1858-11-16 00:00 +************************************************************************/ + long year; + double fractional_year; + double mjd; + + if( decimalyear < -10000 || decimalyear > 10000 ) + return 0; + + year = lround(floor(decimalyear)); + fractional_year = decimalyear - year; + mjd = (year - 1859)*365 + 14 + 31; + mjd += (double)fractional_year*(double)days_in_year(year); + + /* take care of leap days */ + year--; + for (; year > 1858; year--) + if (is_leap_year(year)) + mjd++; + + return mjd; +} + + +/***********************************************************************/ +static double mjd_to_decimalyear(double mjd) { +/*********************************************************************** + Epoch of modified julian date is 1858-11-16 00:00 +************************************************************************/ + double decimalyear = mjd; + double mjd_iter = 14 + 31; + int year = 1859; + + /* a smarter brain than mine could probably to do this more elegantly + - I'll just brute-force my way out of this... */ + for (; mjd >= mjd_iter; year++) { + mjd_iter += days_in_year(year); + } + year--; + mjd_iter -= days_in_year(year); + + decimalyear = year + (mjd-mjd_iter)/days_in_year(year); + return decimalyear; +} + + +/***********************************************************************/ +static double gps_week_to_mjd(double gps_week) { +/*********************************************************************** + GPS weeks are defined as the number of weeks since January the 6th + 1980. + + Epoch of gps weeks is 1980-01-06 00:00, which in modified Julian + date is 44244. +************************************************************************/ + return 44244.0 + gps_week*7.0; +} + + +/***********************************************************************/ +static double mjd_to_gps_week(double mjd) { +/*********************************************************************** + GPS weeks are defined as the number of weeks since January the 6th + 1980. + + Epoch of gps weeks is 1980-01-06 00:00, which in modified Julian + date is 44244. +************************************************************************/ + return (mjd - 44244.0) / 7.0; +} + + +/***********************************************************************/ +static double yyyymmdd_to_mjd(double yyyymmdd) { +/************************************************************************ + Date given in YYYY-MM-DD format. +************************************************************************/ + + long year = lround(floor( yyyymmdd / 10000 )); + long month = lround(floor((yyyymmdd - year*10000) / 100)); + long day = lround(floor( yyyymmdd - year*10000 - month*100)); + double mjd = daynumber_in_year(year, month, day); + + for (year -= 1; year > 1858; year--) + mjd += days_in_year(year); + + return mjd + 13 + 31; +} + + +/***********************************************************************/ +static double mjd_to_yyyymmdd(double mjd) { +/************************************************************************ + Date given in YYYY-MM-DD format. +************************************************************************/ + double mjd_iter = 14 + 31; + int year = 1859, month=0, day=0; + + for (; mjd >= mjd_iter; year++) { + mjd_iter += days_in_year(year); + } + year--; + mjd_iter -= days_in_year(year); + + for (month=1; mjd_iter + days_in_month(year, month) <= mjd; month++) + mjd_iter += days_in_month(year, month); + + day = (int)(mjd - mjd_iter + 1); + + return year*10000.0 + month*100.0 + day; +} + +static const struct TIME_UNITS time_units[] = { + {"mjd", mjd_to_mjd, mjd_to_mjd, "Modified julian date"}, + {"decimalyear", decimalyear_to_mjd, mjd_to_decimalyear, "Decimal year"}, + {"gps_week", gps_week_to_mjd, mjd_to_gps_week, "GPS Week"}, + {"yyyymmdd", yyyymmdd_to_mjd, mjd_to_yyyymmdd, "YYYYMMDD date"}, + {nullptr, nullptr, nullptr, nullptr} +}; + + +/***********************************************************************/ +static XY forward_2d(LP lp, PJ *P) { +/************************************************************************ + Forward unit conversions in the plane +************************************************************************/ + struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + point.lp = lp; + + point.xy.x *= Q->xy_factor; + point.xy.y *= Q->xy_factor; + + return point.xy; +} + + +/***********************************************************************/ +static LP reverse_2d(XY xy, PJ *P) { +/************************************************************************ + Reverse unit conversions in the plane +************************************************************************/ + struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + point.xy = xy; + + point.xy.x /= Q->xy_factor; + point.xy.y /= Q->xy_factor; + + return point.lp; +} + + +/***********************************************************************/ +static XYZ forward_3d(LPZ lpz, PJ *P) { +/************************************************************************ + Forward unit conversions the vertical component +************************************************************************/ + struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + point.lpz = lpz; + + /* take care of the horizontal components in the 2D function */ + point.xy = forward_2d(point.lp, P); + + point.xyz.z *= Q->z_factor; + + return point.xyz; +} + +/***********************************************************************/ +static LPZ reverse_3d(XYZ xyz, PJ *P) { +/************************************************************************ + Reverse unit conversions the vertical component +************************************************************************/ + struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; + PJ_COORD point = {{0,0,0,0}}; + point.xyz = xyz; + + /* take care of the horizontal components in the 2D function */ + point.lp = reverse_2d(point.xy, P); + + point.xyz.z /= Q->z_factor; + + return point.lpz; +} + + +/***********************************************************************/ +static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) { +/************************************************************************ + Forward conversion of time units +************************************************************************/ + struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; + PJ_COORD out = obs; + + /* delegate unit conversion of physical dimensions to the 3D function */ + out.xyz = forward_3d(obs.lpz, P); + + if (Q->t_in_id >= 0) + out.xyzt.t = time_units[Q->t_in_id].t_in( obs.xyzt.t ); + if (Q->t_out_id >= 0) + out.xyzt.t = time_units[Q->t_out_id].t_out( out.xyzt.t ); + + return out; +} + + +/***********************************************************************/ +static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) { +/************************************************************************ + Reverse conversion of time units +************************************************************************/ + struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque; + PJ_COORD out = obs; + + /* delegate unit conversion of physical dimensions to the 3D function */ + out.lpz = reverse_3d(obs.xyz, P); + + if (Q->t_out_id >= 0) + out.xyzt.t = time_units[Q->t_out_id].t_in( obs.xyzt.t ); + if (Q->t_in_id >= 0) + out.xyzt.t = time_units[Q->t_in_id].t_out( out.xyzt.t ); + + return out; +} + +/***********************************************************************/ +static double get_unit_conversion_factor(const char* name, + int* p_is_linear, + const char** p_normalized_name) { +/***********************************************************************/ + int i; + const char* s; + const PJ_UNITS *units; + + units = proj_list_units(); + + /* Try first with linear units */ + for (i = 0; (s = units[i].id) ; ++i) { + if ( strcmp(s, name) == 0 ) { + if( p_normalized_name ) { + *p_normalized_name = units[i].name; + } + if( p_is_linear ) { + *p_is_linear = 1; + } + return units[i].factor; + } + } + + /* And then angular units */ + units = proj_list_angular_units(); + for (i = 0; (s = units[i].id) ; ++i) { + if ( strcmp(s, name) == 0 ) { + if( p_normalized_name ) { + *p_normalized_name = units[i].name; + } + if( p_is_linear ) { + *p_is_linear = 0; + } + return units[i].factor; + } + } + if( p_normalized_name ) { + *p_normalized_name = nullptr; + } + if( p_is_linear ) { + *p_is_linear = -1; + } + return 0.0; +} + +/***********************************************************************/ +PJ *CONVERSION(unitconvert,0) { +/***********************************************************************/ + struct pj_opaque_unitconvert *Q = static_cast(pj_calloc (1, sizeof (struct pj_opaque_unitconvert))); + const char *s, *name; + int i; + double f; + int xy_in_is_linear = -1; /* unknown */ + int xy_out_is_linear = -1; /* unknown */ + int z_in_is_linear = -1; /* unknown */ + int z_out_is_linear = -1; /* unknown */ + + if (nullptr==Q) + return pj_default_destructor (P, ENOMEM); + P->opaque = (void *) Q; + + P->fwd4d = forward_4d; + P->inv4d = reverse_4d; + P->fwd3d = forward_3d; + P->inv3d = reverse_3d; + P->fwd = forward_2d; + P->inv = reverse_2d; + + P->left = PJ_IO_UNITS_WHATEVER; + P->right = PJ_IO_UNITS_WHATEVER; + + /* if no time input/output unit is specified we can skip them */ + Q->t_in_id = -1; + Q->t_out_id = -1; + + Q->xy_factor = 1.0; + Q->z_factor = 1.0; + + if ((name = pj_param (P->ctx, P->params, "sxy_in").s) != nullptr) { + const char* normalized_name = nullptr; + f = get_unit_conversion_factor(name, &xy_in_is_linear, &normalized_name); + if (f != 0.0) { + proj_log_debug(P, "xy_in unit: %s", normalized_name); + } else { + if ( (f = pj_param (P->ctx, P->params, "dxy_in").f) == 0.0) + return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); + } + if (f != 0.0) + Q->xy_factor *= f; + } + + if ((name = pj_param (P->ctx, P->params, "sxy_out").s) != nullptr) { + const char* normalized_name = nullptr; + f = get_unit_conversion_factor(name, &xy_out_is_linear, &normalized_name); + if (f != 0.0) { + proj_log_debug(P, "xy_out unit: %s", normalized_name); + } else { + if ( (f = pj_param (P->ctx, P->params, "dxy_out").f) == 0.0) + return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); + } + if (f != 0.0) + Q->xy_factor /= f; + } + + if( xy_in_is_linear >= 0 && xy_out_is_linear >= 0 && + xy_in_is_linear != xy_out_is_linear ) { + proj_log_debug(P, "inconsistent unit type between xy_in and xy_out"); + return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT); + } + + if ((name = pj_param (P->ctx, P->params, "sz_in").s) != nullptr) { + const char* normalized_name = nullptr; + f = get_unit_conversion_factor(name, &z_in_is_linear, &normalized_name); + if (f != 0.0) { + proj_log_debug(P, "z_in unit: %s", normalized_name); + } else { + if ( (f = pj_param (P->ctx, P->params, "dz_in").f) == 0.0) + return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); + } + if (f != 0.0) + Q->z_factor *= f; + } + + if ((name = pj_param (P->ctx, P->params, "sz_out").s) != nullptr) { + const char* normalized_name = nullptr; + f = get_unit_conversion_factor(name, &z_out_is_linear, &normalized_name); + if (f != 0.0) { + proj_log_debug(P, "z_out unit: %s", normalized_name); + } else { + if ( (f = pj_param (P->ctx, P->params, "dz_out").f) == 0.0) + return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); + } + if (f != 0.0) + Q->z_factor /= f; + } + + if( z_in_is_linear >= 0 && z_out_is_linear >= 0 && + z_in_is_linear != z_out_is_linear ) { + proj_log_debug(P, "inconsistent unit type between z_in and z_out"); + return pj_default_destructor(P, PJD_ERR_INCONSISTENT_UNIT); + } + + if ((name = pj_param (P->ctx, P->params, "st_in").s) != nullptr) { + for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i); + + if (!s) return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); /* unknown unit conversion id */ + + Q->t_in_id = i; + proj_log_debug(P, "t_in unit: %s", time_units[i].name); + } + + s = nullptr; + if ((name = pj_param (P->ctx, P->params, "st_out").s) != nullptr) { + for (i = 0; (s = time_units[i].id) && strcmp(name, s) ; ++i); + + if (!s) return pj_default_destructor(P, PJD_ERR_UNKNOWN_UNIT_ID); /* unknown unit conversion id */ + + Q->t_out_id = i; + proj_log_debug(P, "t_out unit: %s", time_units[i].name); + } + + return P; +} -- cgit v1.2.3