From e2dd223c5601b387cd622b88583e857b03e9fade Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sat, 14 Nov 2020 11:28:36 +0100 Subject: Remove src/transform.cpp and related tests --- src/Makefile.am | 2 +- src/lib_proj.cmake | 1 - src/malloc.cpp | 1 - src/proj_internal.h | 12 - src/transform.cpp | 1850 --------------------------------------- test/unit/CMakeLists.txt | 11 - test/unit/Makefile.am | 10 +- test/unit/pj_transform_test.cpp | 737 ---------------- 8 files changed, 3 insertions(+), 2621 deletions(-) delete mode 100644 src/transform.cpp delete mode 100644 test/unit/pj_transform_test.cpp diff --git a/src/Makefile.am b/src/Makefile.am index 62bff403..edfd05ef 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -220,7 +220,7 @@ libproj_la_SOURCES = \ generic_inverse.cpp \ quadtree.hpp \ \ - datums.cpp datum_set.cpp transform.cpp \ + datums.cpp datum_set.cpp \ utils.cpp \ mutex.cpp initcache.cpp geodesic.c \ strtod.cpp \ diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake index 95417004..b0f9a6ba 100644 --- a/src/lib_proj.cmake +++ b/src/lib_proj.cmake @@ -246,7 +246,6 @@ set(SRC_LIBPROJ_CORE rtodms.cpp strerrno.cpp strtod.cpp - transform.cpp tsfn.cpp units.cpp utils.cpp diff --git a/src/malloc.cpp b/src/malloc.cpp index c8de6630..eb64023f 100644 --- a/src/malloc.cpp +++ b/src/malloc.cpp @@ -262,7 +262,6 @@ PJ *pj_default_destructor (PJ *P, int errlev) { /* Destructor */ void proj_cleanup() { /*****************************************************************************/ pj_clear_initcache(); - pj_deallocate_grids(); FileManager::clearMemoryCache(); pj_clear_hgridshift_knowngrids_cache(); pj_clear_vgridshift_knowngrids_cache(); diff --git a/src/proj_internal.h b/src/proj_internal.h index e0cc5170..062f240b 100644 --- a/src/proj_internal.h +++ b/src/proj_internal.h @@ -985,21 +985,9 @@ projXYZ PROJ_DLL pj_fwd3d(projLPZ, projPJ); projLPZ PROJ_DLL pj_inv3d(projXYZ, projPJ); -int PROJ_DLL pj_transform( projPJ src, projPJ dst, long point_count, int point_offset, - double *x, double *y, double *z ); -int PROJ_DLL pj_datum_transform( projPJ src, projPJ dst, long point_count, int point_offset, - double *x, double *y, double *z ); -int PROJ_DLL pj_geocentric_to_geodetic( double a, double es, - long point_count, int point_offset, - double *x, double *y, double *z ); -int PROJ_DLL pj_geodetic_to_geocentric( double a, double es, - long point_count, int point_offset, - double *x, double *y, double *z ); -int PROJ_DLL pj_compare_datums( projPJ srcdefn, projPJ dstdefn ); int PROJ_DLL pj_apply_gridshift( projCtx, const char *, int, long point_count, int point_offset, double *x, double *y, double *z ); -void PROJ_DLL pj_deallocate_grids(void); void PROJ_DLL pj_clear_initcache(void); int PROJ_DLL pj_is_latlong(projPJ); int PROJ_DLL pj_is_geocent(projPJ); diff --git a/src/transform.cpp b/src/transform.cpp deleted file mode 100644 index cff89232..00000000 --- a/src/transform.cpp +++ /dev/null @@ -1,1850 +0,0 @@ -/****************************************************************************** - * Project: PROJ.4 - * Purpose: Perform overall coordinate system to coordinate system - * transformations (pj_transform() function) including reprojection - * and datum shifting. - * Author: Frank Warmerdam, warmerdam@pobox.com - * - ****************************************************************************** - * Copyright (c) 2000, 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. - *****************************************************************************/ - -#include -#include -#include - -#include "proj.h" -#include "proj_internal.h" -#include "grids.hpp" - -using namespace NS_PROJ; - - -/////////////////////////////////////////////////////////////////////////////// -/// From older geocent.h -/////////////////////////////////////////////////////////////////////////////// - - -/***************************************************************************/ -/* RSC IDENTIFIER: GEOCENTRIC - * - * ABSTRACT - * - * This component provides conversions between Geodetic coordinates (latitude, - * longitude in radians and height in meters) and Geocentric coordinates - * (X, Y, Z) in meters. - * - * ERROR HANDLING - * - * This component checks parameters for valid values. If an invalid value - * is found, the error code is combined with the current error code using - * the bitwise or. This combining allows multiple error codes to be - * returned. The possible error codes are: - * - * GEOCENT_NO_ERROR : No errors occurred in function - * GEOCENT_LAT_ERROR : Latitude out of valid range - * (-90 to 90 degrees) - * GEOCENT_LON_ERROR : Longitude out of valid range - * (-180 to 360 degrees) - * GEOCENT_A_ERROR : Semi-major axis less than or equal to zero - * GEOCENT_B_ERROR : Semi-minor axis less than or equal to zero - * GEOCENT_A_LESS_B_ERROR : Semi-major axis less than semi-minor axis - * - * - * REUSE NOTES - * - * GEOCENTRIC is intended for reuse by any application that performs - * coordinate conversions between geodetic coordinates and geocentric - * coordinates. - * - * - * REFERENCES - * - * An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion, - * Ralph Toms, February 1996 UCRL-JC-123138. - * - * Further information on GEOCENTRIC can be found in the Reuse Manual. - * - * GEOCENTRIC originated from : U.S. Army Topographic Engineering Center - * Geospatial Information Division - * 7701 Telegraph Road - * Alexandria, VA 22310-3864 - * - * LICENSES - * - * None apply to this component. - * - * RESTRICTIONS - * - * GEOCENTRIC has no restrictions. - * - * ENVIRONMENT - * - * GEOCENTRIC was tested and certified in the following environments: - * - * 1. Solaris 2.5 with GCC version 2.8.1 - * 2. Windows 95 with MS Visual C++ version 6 - * - * MODIFICATIONS - * - * Date Description - * ---- ----------- - * - * - */ - - -/***************************************************************************/ -/* - * DEFINES - */ -#define GEOCENT_NO_ERROR 0x0000 -#define GEOCENT_LAT_ERROR 0x0001 -#define GEOCENT_LON_ERROR 0x0002 -#define GEOCENT_A_ERROR 0x0004 -#define GEOCENT_B_ERROR 0x0008 -#define GEOCENT_A_LESS_B_ERROR 0x0010 - - -/***************************************************************************/ -/* - * FUNCTION PROTOTYPES - */ - -/* ensure proper linkage to c++ programs */ -#ifdef __cplusplus -extern "C" { -#endif - -typedef struct -{ - double Geocent_a; /* Semi-major axis of ellipsoid in meters */ - double Geocent_b; /* Semi-minor axis of ellipsoid */ - double Geocent_a2; /* Square of semi-major axis */ - double Geocent_b2; /* Square of semi-minor axis */ - double Geocent_e2; /* Eccentricity squared */ - double Geocent_ep2; /* 2nd eccentricity squared */ -} GeocentricInfo; - -void pj_Init_Geocentric( GeocentricInfo *gi ); -long pj_Set_Geocentric_Parameters( GeocentricInfo *gi, - double a, - double b); - -/* - * The function Set_Geocentric_Parameters receives the ellipsoid parameters - * as inputs and sets the corresponding state variables. - * - * a : Semi-major axis, in meters. (input) - * b : Semi-minor axis, in meters. (input) - */ - - -void pj_Get_Geocentric_Parameters ( GeocentricInfo *gi, - double *a, - double *b); - -/* - * The function Get_Geocentric_Parameters returns the ellipsoid parameters - * to be used in geocentric coordinate conversions. - * - * a : Semi-major axis, in meters. (output) - * b : Semi-minor axis, in meters. (output) - */ - - -long pj_Convert_Geodetic_To_Geocentric ( GeocentricInfo *gi, - double Latitude, - double Longitude, - double Height, - double *X, - double *Y, - double *Z); -/* - * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates - * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z), - * according to the current ellipsoid parameters. - * - * Latitude : Geodetic latitude in radians (input) - * Longitude : Geodetic longitude in radians (input) - * Height : Geodetic height, in meters (input) - * X : Calculated Geocentric X coordinate, in meters. (output) - * Y : Calculated Geocentric Y coordinate, in meters. (output) - * Z : Calculated Geocentric Z coordinate, in meters. (output) - * - */ - - -void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi, - double X, - double Y, - double Z, - double *Latitude, - double *Longitude, - double *Height); -/* - * The function Convert_Geocentric_To_Geodetic converts geocentric - * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, - * and height), according to the current ellipsoid parameters. - * - * X : Geocentric X coordinate, in meters. (input) - * Y : Geocentric Y coordinate, in meters. (input) - * Z : Geocentric Z coordinate, in meters. (input) - * Latitude : Calculated latitude value in radians. (output) - * Longitude : Calculated longitude value in radians. (output) - * Height : Calculated height value, in meters. (output) - */ - - -#ifdef __cplusplus -} -#endif - - -/////////////////////////////////////////////////////////////////////////////// -/// From older geocent.cpp -/////////////////////////////////////////////////////////////////////////////// - - -/***************************************************************************/ -/* RSC IDENTIFIER: GEOCENTRIC - * - * ABSTRACT - * - * This component provides conversions between Geodetic coordinates (latitude, - * longitude in radians and height in meters) and Geocentric coordinates - * (X, Y, Z) in meters. - * - * ERROR HANDLING - * - * This component checks parameters for valid values. If an invalid value - * is found, the error code is combined with the current error code using - * the bitwise or. This combining allows multiple error codes to be - * returned. The possible error codes are: - * - * GEOCENT_NO_ERROR : No errors occurred in function - * GEOCENT_LAT_ERROR : Latitude out of valid range - * (-90 to 90 degrees) - * GEOCENT_LON_ERROR : Longitude out of valid range - * (-180 to 360 degrees) - * GEOCENT_A_ERROR : Semi-major axis lessthan or equal to zero - * GEOCENT_B_ERROR : Semi-minor axis lessthan or equal to zero - * GEOCENT_A_LESS_B_ERROR : Semi-major axis less than semi-minor axis - * - * - * REUSE NOTES - * - * GEOCENTRIC is intended for reuse by any application that performs - * coordinate conversions between geodetic coordinates and geocentric - * coordinates. - * - * - * REFERENCES - * - * An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion, - * Ralph Toms, February 1996 UCRL-JC-123138. - * - * Further information on GEOCENTRIC can be found in the Reuse Manual. - * - * GEOCENTRIC originated from : U.S. Army Topographic Engineering Center - * Geospatial Information Division - * 7701 Telegraph Road - * Alexandria, VA 22310-3864 - * - * LICENSES - * - * None apply to this component. - * - * RESTRICTIONS - * - * GEOCENTRIC has no restrictions. - * - * ENVIRONMENT - * - * GEOCENTRIC was tested and certified in the following environments: - * - * 1. Solaris 2.5 with GCC version 2.8.1 - * 2. Windows 95 with MS Visual C++ version 6 - * - * MODIFICATIONS - * - * Date Description - * ---- ----------- - * 25-02-97 Original Code - * - */ - - -/***************************************************************************/ -/* - * INCLUDES - */ -#include -//#include "geocent.h" -/* - * math.h - is needed for calls to sin, cos, tan and sqrt. - * geocent.h - is needed for Error codes and prototype error checking. - */ - - -/***************************************************************************/ -/* - * DEFINES - */ -#define PI 3.14159265358979323e0 -#define PI_OVER_2 (PI / 2.0e0) -#define FALSE 0 -#define TRUE 1 -#define COS_67P5 0.38268343236508977 /* cosine of 67.5 degrees */ -#define AD_C 1.0026000 /* Toms region 1 constant */ - - -/***************************************************************************/ -/* - * FUNCTIONS - */ - - -long pj_Set_Geocentric_Parameters (GeocentricInfo *gi, double a, double b) - -{ /* BEGIN Set_Geocentric_Parameters */ -/* - * The function Set_Geocentric_Parameters receives the ellipsoid parameters - * as inputs and sets the corresponding state variables. - * - * a : Semi-major axis, in meters. (input) - * b : Semi-minor axis, in meters. (input) - */ - long Error_Code = GEOCENT_NO_ERROR; - - if (a <= 0.0) - Error_Code |= GEOCENT_A_ERROR; - if (b <= 0.0) - Error_Code |= GEOCENT_B_ERROR; - if (a < b) - Error_Code |= GEOCENT_A_LESS_B_ERROR; - if (!Error_Code) - { - gi->Geocent_a = a; - gi->Geocent_b = b; - gi->Geocent_a2 = a * a; - gi->Geocent_b2 = b * b; - gi->Geocent_e2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_a2; - gi->Geocent_ep2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_b2; - } - return (Error_Code); -} /* END OF Set_Geocentric_Parameters */ - - -void pj_Get_Geocentric_Parameters (GeocentricInfo *gi, - double *a, - double *b) -{ /* BEGIN Get_Geocentric_Parameters */ -/* - * The function Get_Geocentric_Parameters returns the ellipsoid parameters - * to be used in geocentric coordinate conversions. - * - * a : Semi-major axis, in meters. (output) - * b : Semi-minor axis, in meters. (output) - */ - - *a = gi->Geocent_a; - *b = gi->Geocent_b; -} /* END OF Get_Geocentric_Parameters */ - - -long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo *gi, - double Latitude, - double Longitude, - double Height, - double *X, - double *Y, - double *Z) -{ /* BEGIN Convert_Geodetic_To_Geocentric */ -/* - * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates - * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z), - * according to the current ellipsoid parameters. - * - * Latitude : Geodetic latitude in radians (input) - * Longitude : Geodetic longitude in radians (input) - * Height : Geodetic height, in meters (input) - * X : Calculated Geocentric X coordinate, in meters (output) - * Y : Calculated Geocentric Y coordinate, in meters (output) - * Z : Calculated Geocentric Z coordinate, in meters (output) - * - */ - long Error_Code = GEOCENT_NO_ERROR; - double Rn; /* Earth radius at location */ - double Sin_Lat; /* sin(Latitude) */ - double Sin2_Lat; /* Square of sin(Latitude) */ - double Cos_Lat; /* cos(Latitude) */ - - /* - ** Don't blow up if Latitude is just a little out of the value - ** range as it may just be a rounding issue. Also removed longitude - ** test, it should be wrapped by cos() and sin(). NFW for PROJ.4, Sep/2001. - */ - if( Latitude < -PI_OVER_2 && Latitude > -1.001 * PI_OVER_2 ) - Latitude = -PI_OVER_2; - else if( Latitude > PI_OVER_2 && Latitude < 1.001 * PI_OVER_2 ) - Latitude = PI_OVER_2; - else if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2)) - { /* Latitude out of range */ - Error_Code |= GEOCENT_LAT_ERROR; - } - - if (!Error_Code) - { /* no errors */ - if (Longitude > PI) - Longitude -= (2*PI); - Sin_Lat = sin(Latitude); - Cos_Lat = cos(Latitude); - Sin2_Lat = Sin_Lat * Sin_Lat; - Rn = gi->Geocent_a / (sqrt(1.0e0 - gi->Geocent_e2 * Sin2_Lat)); - *X = (Rn + Height) * Cos_Lat * cos(Longitude); - *Y = (Rn + Height) * Cos_Lat * sin(Longitude); - *Z = ((Rn * (1 - gi->Geocent_e2)) + Height) * Sin_Lat; - - } - return (Error_Code); -} /* END OF Convert_Geodetic_To_Geocentric */ - -/* - * The function Convert_Geocentric_To_Geodetic converts geocentric - * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, - * and height), according to the current ellipsoid parameters. - * - * X : Geocentric X coordinate, in meters. (input) - * Y : Geocentric Y coordinate, in meters. (input) - * Z : Geocentric Z coordinate, in meters. (input) - * Latitude : Calculated latitude value in radians. (output) - * Longitude : Calculated longitude value in radians. (output) - * Height : Calculated height value, in meters. (output) - */ - -#define USE_ITERATIVE_METHOD - -void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi, - double X, - double Y, - double Z, - double *Latitude, - double *Longitude, - double *Height) -{ /* BEGIN Convert_Geocentric_To_Geodetic */ -#if !defined(USE_ITERATIVE_METHOD) -/* - * The method used here is derived from 'An Improved Algorithm for - * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996 - */ - -/* Note: Variable names follow the notation used in Toms, Feb 1996 */ - - double W; /* distance from Z axis */ - double W2; /* square of distance from Z axis */ - double T0; /* initial estimate of vertical component */ - double T1; /* corrected estimate of vertical component */ - double S0; /* initial estimate of horizontal component */ - double S1; /* corrected estimate of horizontal component */ - double Sin_B0; /* sin(B0), B0 is estimate of Bowring aux variable */ - double Sin3_B0; /* cube of sin(B0) */ - double Cos_B0; /* cos(B0) */ - double Sin_p1; /* sin(phi1), phi1 is estimated latitude */ - double Cos_p1; /* cos(phi1) */ - double Rn; /* Earth radius at location */ - double Sum; /* numerator of cos(phi1) */ - int At_Pole; /* indicates location is in polar region */ - - At_Pole = FALSE; - if (X != 0.0) - { - *Longitude = atan2(Y,X); - } - else - { - if (Y > 0) - { - *Longitude = PI_OVER_2; - } - else if (Y < 0) - { - *Longitude = -PI_OVER_2; - } - else - { - At_Pole = TRUE; - *Longitude = 0.0; - if (Z > 0.0) - { /* north pole */ - *Latitude = PI_OVER_2; - } - else if (Z < 0.0) - { /* south pole */ - *Latitude = -PI_OVER_2; - } - else - { /* center of earth */ - *Latitude = PI_OVER_2; - *Height = -Geocent_b; - return; - } - } - } - W2 = X*X + Y*Y; - W = sqrt(W2); - T0 = Z * AD_C; - S0 = sqrt(T0 * T0 + W2); - Sin_B0 = T0 / S0; - Cos_B0 = W / S0; - Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0; - T1 = Z + gi->Geocent_b * gi->Geocent_ep2 * Sin3_B0; - Sum = W - gi->Geocent_a * gi->Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0; - S1 = sqrt(T1*T1 + Sum * Sum); - Sin_p1 = T1 / S1; - Cos_p1 = Sum / S1; - Rn = gi->Geocent_a / sqrt(1.0 - gi->Geocent_e2 * Sin_p1 * Sin_p1); - if (Cos_p1 >= COS_67P5) - { - *Height = W / Cos_p1 - Rn; - } - else if (Cos_p1 <= -COS_67P5) - { - *Height = W / -Cos_p1 - Rn; - } - else - { - *Height = Z / Sin_p1 + Rn * (gi->Geocent_e2 - 1.0); - } - if (At_Pole == FALSE) - { - *Latitude = atan(Sin_p1 / Cos_p1); - } -#else /* defined(USE_ITERATIVE_METHOD) */ -/* -* Reference... -* ============ -* Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für -* das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover -* Nr. 137, p. 130-131. - -* Programmed by GGA- Leibniz-Institute of Applied Geophysics -* Stilleweg 2 -* D-30655 Hannover -* Federal Republic of Germany -* Internet: www.gga-hannover.de -* -* Hannover, March 1999, April 2004. -* see also: comments in statements -* remarks: -* Mathematically exact and because of symmetry of rotation-ellipsoid, -* each point (X,Y,Z) has at least two solutions (Latitude1,Longitude1,Height1) and -* (Latitude2,Longitude2,Height2). Is point=(0.,0.,Z) (P=0.), so you get even -* four solutions, every two symmetrical to the semi-minor axis. -* Here Height1 and Height2 have at least a difference in order of -* radius of curvature (e.g. (0,0,b)=> (90.,0.,0.) or (-90.,0.,-2b); -* (a+100.)*(sqrt(2.)/2.,sqrt(2.)/2.,0.) => (0.,45.,100.) or -* (0.,225.,-(2a+100.))). -* The algorithm always computes (Latitude,Longitude) with smallest |Height|. -* For normal computations, that means |Height|<10000.m, algorithm normally -* converges after to 2-3 steps!!! -* But if |Height| has the amount of length of ellipsoid's axis -* (e.g. -6300000.m), algorithm needs about 15 steps. -*/ - -/* local definitions and variables */ -/* end-criterium of loop, accuracy of sin(Latitude) */ -#define genau 1.E-12 -#define genau2 (genau*genau) -#define maxiter 30 - - double P; /* distance between semi-minor axis and location */ - double RR; /* distance between center and location */ - double CT; /* sin of geocentric latitude */ - double ST; /* cos of geocentric latitude */ - double RX; - double RK; - double RN; /* Earth radius at location */ - double CPHI0; /* cos of start or old geodetic latitude in iterations */ - double SPHI0; /* sin of start or old geodetic latitude in iterations */ - double CPHI; /* cos of searched geodetic latitude */ - double SPHI; /* sin of searched geodetic latitude */ - double SDPHI; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */ - int iter; /* # of continuous iteration, max. 30 is always enough (s.a.) */ - - P = sqrt(X*X+Y*Y); - RR = sqrt(X*X+Y*Y+Z*Z); - -/* special cases for latitude and longitude */ - if (P/gi->Geocent_a < genau) { - -/* special case, if P=0. (X=0., Y=0.) */ - *Longitude = 0.; - -/* if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis - * of ellipsoid (=center of mass), Latitude becomes PI/2 */ - if (RR/gi->Geocent_a < genau) { - *Latitude = PI_OVER_2; - *Height = -gi->Geocent_b; - return ; - - } - } - else { -/* ellipsoidal (geodetic) longitude - * interval: -PI < Longitude <= +PI */ - *Longitude=atan2(Y,X); - } - -/* -------------------------------------------------------------- - * Following iterative algorithm was developed by - * "Institut für Erdmessung", University of Hannover, July 1988. - * Internet: www.ife.uni-hannover.de - * Iterative computation of CPHI,SPHI and Height. - * Iteration of CPHI and SPHI to 10**-12 radian resp. - * 2*10**-7 arcsec. - * -------------------------------------------------------------- - */ - CT = Z/RR; - ST = P/RR; - { - const double denominator = 1.0-gi->Geocent_e2*(2.0-gi->Geocent_e2)*ST*ST; - if( denominator == 0 ) - { - *Latitude = HUGE_VAL; - *Longitude = HUGE_VAL; - *Height = HUGE_VAL; - return; - } - RX = 1.0/sqrt(denominator); - } - CPHI0 = ST*(1.0-gi->Geocent_e2)*RX; - SPHI0 = CT*RX; - iter = 0; - -/* loop to find sin(Latitude) resp. Latitude - * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */ - do - { - iter++; - RN = gi->Geocent_a/sqrt(1.0-gi->Geocent_e2*SPHI0*SPHI0); - -/* ellipsoidal (geodetic) height */ - *Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi->Geocent_e2*SPHI0*SPHI0); - - /* avoid zero division */ - if (RN+*Height==0.0) { - *Latitude = 0.0; - return; - } - RK = gi->Geocent_e2*RN/(RN+*Height); - { - const double denominator = 1.0-RK*(2.0-RK)*ST*ST; - if( denominator == 0 ) - { - *Latitude = HUGE_VAL; - *Longitude = HUGE_VAL; - *Height = HUGE_VAL; - return; - } - RX = 1.0/sqrt(denominator); - } - CPHI = ST*(1.0-RK)*RX; - SPHI = CT*RX; - SDPHI = SPHI*CPHI0-CPHI*SPHI0; - CPHI0 = CPHI; - SPHI0 = SPHI; - } - while (SDPHI*SDPHI > genau2 && iter < maxiter); - -/* ellipsoidal (geodetic) latitude */ - *Latitude=atan2(SPHI, fabs(CPHI)); - -#endif /* defined(USE_ITERATIVE_METHOD) */ -} /* END OF Convert_Geocentric_To_Geodetic */ - - -/////////////////////////////////////////////////////////////////////////////// -/// Main of transform.cpp -/////////////////////////////////////////////////////////////////////////////// - - -static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag, - long point_count, int point_offset, - double *x, double *y, double *z ); - -#ifndef SRS_WGS84_SEMIMAJOR -#define SRS_WGS84_SEMIMAJOR 6378137.0 -#endif - -#ifndef SRS_WGS84_ESQUARED -#define SRS_WGS84_ESQUARED 0.0066943799901413165 -#endif - -#define Dx_BF (defn->datum_params[0]) -#define Dy_BF (defn->datum_params[1]) -#define Dz_BF (defn->datum_params[2]) -#define Rx_BF (defn->datum_params[3]) -#define Ry_BF (defn->datum_params[4]) -#define Rz_BF (defn->datum_params[5]) -#define M_BF (defn->datum_params[6]) - -/* -** This table is intended to indicate for any given error code -** whether that error will occur for all locations (ie. -** it is a problem with the coordinate system as a whole) in which case the -** value would be 0, or if the problem is with the point being transformed -** in which case the value is 1. -** -** At some point we might want to move this array in with the error message -** list or something, but while experimenting with it this should be fine. -** -** -** NOTE (2017-10-01): Non-transient errors really should have resulted in a -** PJ==0 during initialization, and hence should be handled at the level -** before calling pj_transform. The only obvious example of the contrary -** appears to be the PJD_ERR_GRID_AREA case, which may also be taken to -** mean "no grids available" -** -** -*/ - -static const int transient_error[70] = { - /* 0 1 2 3 4 5 6 7 8 9 */ - /* 0 to 9 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, - /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, - /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, - /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0, - /* 60 to 69 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,}; - - -/* -------------------------------------------------------------------- */ -/* Read transient_error[] in a safe way. */ -/* -------------------------------------------------------------------- */ -static int get_transient_error_value(int pos_index) -{ - const int array_size = - (int)(sizeof(transient_error) / sizeof(transient_error[0])); - if( pos_index < 0 || pos_index >= array_size ) { - return 0; - } - return transient_error[pos_index]; -} - - -/* -------------------------------------------------------------------- */ -/* Transform unusual input coordinate axis orientation to */ -/* standard form if needed. */ -/* -------------------------------------------------------------------- */ -static int adjust_axes (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) { - /* Nothing to do? */ - if (0==strcmp(P->axis,"enu")) - return 0; - - return adjust_axis( P->ctx, P->axis, - dir==PJ_FWD ? 1: 0, n, dist, x, y, z ); -} - - - -/* ----------------------------------------------------------------------- */ -/* Transform geographic (lat/long) source coordinates to */ -/* cartesian ("geocentric"), if needed */ -/* ----------------------------------------------------------------------- */ -static int geographic_to_cartesian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) { - int res; - long i; - double fac = P->to_meter; - - /* Nothing to do? */ - if (!P->is_geocent) - return 0; - - if ( z == nullptr ) { - pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC); - return PJD_ERR_GEOCENTRIC; - } - - if (PJ_FWD==dir) { - fac = P->fr_meter; - res = pj_geodetic_to_geocentric( P->a_orig, P->es_orig, n, dist, x, y, z ); - if (res) - return res; - } - - if (fac != 1.0) { - for( i = 0; i < n; i++ ) { - if( x[dist*i] != HUGE_VAL ) { - x[dist*i] *= fac; - y[dist*i] *= fac; - z[dist*i] *= fac; - } - } - } - - if (PJ_FWD==dir) - return 0; - return pj_geocentric_to_geodetic( - P->a_orig, P->es_orig, - n, dist, - x, y, z - ); -} - - - - - - - - - - -/* -------------------------------------------------------------------- */ -/* Transform destination points to projection coordinates, if */ -/* desired. */ -/* */ -/* Ought to fold this into projected_to_geographic */ -/* -------------------------------------------------------------------- */ -static int geographic_to_projected (PJ *P, long n, int dist, double *x, double *y, double *z) { - long i; - - /* Nothing to do? */ - if (P->is_latlong && !P->geoc && P->vto_meter == 1.0) - return 0; - if (P->is_geocent) - return 0; - - if(P->fwd3d != nullptr && !(z == nullptr && P->is_latlong)) - { - /* Three dimensions must be defined */ - if ( z == nullptr) - { - pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC); - return PJD_ERR_GEOCENTRIC; - } - - for( i = 0; i < n; i++ ) - { - PJ_XYZ projected_loc; - PJ_LPZ geodetic_loc; - - geodetic_loc.lam = x[dist*i]; - geodetic_loc.phi = y[dist*i]; - geodetic_loc.z = z[dist*i]; - - if (geodetic_loc.lam == HUGE_VAL) - continue; - - proj_errno_reset( P ); - projected_loc = pj_fwd3d( geodetic_loc, P); - if( P->ctx->last_errno != 0 ) - { - if( (P->ctx->last_errno != EDOM - && P->ctx->last_errno != ERANGE) - && (P->ctx->last_errno > 0 - || P->ctx->last_errno < -44 || n == 1 - || get_transient_error_value(-P->ctx->last_errno) == 0 ) ) - { - return P->ctx->last_errno; - } - else - { - projected_loc.x = HUGE_VAL; - projected_loc.y = HUGE_VAL; - projected_loc.z = HUGE_VAL; - } - } - - x[dist*i] = projected_loc.x; - y[dist*i] = projected_loc.y; - z[dist*i] = projected_loc.z; - } - return 0; - } - - // Ugly hack. See https://github.com/OSGeo/PROJ/issues/1782 - if( P->right == PJ_IO_UNITS_WHATEVER && P->descr && - strncmp(P->descr, "General Oblique Transformation", - strlen("General Oblique Transformation")) == 0 ) - { - P->right = PJ_IO_UNITS_PROJECTED; - } - - for( i = 0; i ctx->last_errno != 0 ) - { - if( (P->ctx->last_errno != EDOM - && P->ctx->last_errno != ERANGE) - && (P->ctx->last_errno > 0 - || P->ctx->last_errno < -44 || n == 1 - || get_transient_error_value(-P->ctx->last_errno) == 0 ) ) - { - return P->ctx->last_errno; - } - else - { - projected_loc.x = HUGE_VAL; - projected_loc.y = HUGE_VAL; - } - } - - x[dist*i] = projected_loc.x; - y[dist*i] = projected_loc.y; - } - return 0; -} - - - - - -/* ----------------------------------------------------------------------- */ -/* Transform projected source coordinates to lat/long, if needed */ -/* ----------------------------------------------------------------------- */ -static int projected_to_geographic (PJ *P, long n, int dist, double *x, double *y, double *z) { - long i; - - /* Nothing to do? */ - if (P->is_latlong && !P->geoc && P->vto_meter == 1.0) - return 0; - if (P->is_geocent) - return 0; - - /* Check first if projection is invertible. */ - if( (P->inv3d == nullptr) && (P->inv == nullptr)) - { - pj_ctx_set_errno(pj_get_ctx(P), PJD_ERR_NON_CONV_INV_MERI_DIST); - pj_log( pj_get_ctx(P), PJ_LOG_ERROR, - "pj_transform(): source projection not invertable" ); - return PJD_ERR_NON_CONV_INV_MERI_DIST; - } - - /* If invertible - First try inv3d if defined */ - if (P->inv3d != nullptr && !(z == nullptr && P->is_latlong)) - { - /* Three dimensions must be defined */ - if ( z == nullptr) - { - pj_ctx_set_errno( pj_get_ctx(P), PJD_ERR_GEOCENTRIC); - return PJD_ERR_GEOCENTRIC; - } - - for (i=0; i < n; i++) - { - PJ_XYZ projected_loc; - PJ_LPZ geodetic_loc; - - projected_loc.x = x[dist*i]; - projected_loc.y = y[dist*i]; - projected_loc.z = z[dist*i]; - - if (projected_loc.x == HUGE_VAL) - continue; - - proj_errno_reset( P ); - geodetic_loc = pj_inv3d(projected_loc, P); - if( P->ctx->last_errno != 0 ) - { - if( (P->ctx->last_errno != EDOM - && P->ctx->last_errno != ERANGE) - && (P->ctx->last_errno > 0 - || P->ctx->last_errno < -44 || n == 1 - || get_transient_error_value(-P->ctx->last_errno) == 0 ) ) - { - return P->ctx->last_errno; - } - else - { - geodetic_loc.lam = HUGE_VAL; - geodetic_loc.phi = HUGE_VAL; - geodetic_loc.z = HUGE_VAL; - } - } - - x[dist*i] = geodetic_loc.lam; - y[dist*i] = geodetic_loc.phi; - z[dist*i] = geodetic_loc.z; - - } - return 0; - } - - // Ugly hack. See https://github.com/OSGeo/PROJ/issues/1782 - if( P->right == PJ_IO_UNITS_WHATEVER && P->descr && - strncmp(P->descr, "General Oblique Transformation", - strlen("General Oblique Transformation")) == 0 ) - { - P->right = PJ_IO_UNITS_PROJECTED; - } - - /* Fallback to the original PROJ.4 API 2d inversion - inv */ - for( i = 0; i < n; i++ ) { - PJ_XY projected_loc; - PJ_LP geodetic_loc; - - projected_loc.x = x[dist*i]; - projected_loc.y = y[dist*i]; - - if( projected_loc.x == HUGE_VAL ) - continue; - - proj_errno_reset( P ); - geodetic_loc = pj_inv( projected_loc, P ); - if( P->ctx->last_errno != 0 ) - { - if( (P->ctx->last_errno != EDOM - && P->ctx->last_errno != ERANGE) - && (P->ctx->last_errno > 0 - || P->ctx->last_errno < -44 || n == 1 - || get_transient_error_value(-P->ctx->last_errno) == 0 ) ) - { - return P->ctx->last_errno; - } - else - { - geodetic_loc.lam = HUGE_VAL; - geodetic_loc.phi = HUGE_VAL; - } - } - - x[dist*i] = geodetic_loc.lam; - y[dist*i] = geodetic_loc.phi; - } - return 0; -} - - - -/* -------------------------------------------------------------------- */ -/* Adjust for the prime meridian if needed. */ -/* -------------------------------------------------------------------- */ -static int prime_meridian (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x) { - int i; - double pm = P->from_greenwich; - - /* Nothing to do? */ - if (pm==0.0) - return 0; - if (!(P->is_geocent || P->is_latlong)) - return 0; - - if (dir==PJ_FWD) - pm = -pm; - - for (i = 0; i < n; i++) - if (x[dist*i] != HUGE_VAL) - x[dist*i] += pm; - - return 0; -} - - - -/* -------------------------------------------------------------------- */ -/* Adjust for vertical scale factor if needed */ -/* -------------------------------------------------------------------- */ -static int height_unit (PJ *P, PJ_DIRECTION dir, long n, int dist, double *z) { - int i; - double fac = P->vto_meter; - - if (PJ_FWD==dir) - fac = P->vfr_meter; - - /* Nothing to do? */ - if (fac==1.0) - return 0; - if (nullptr==z) - return 0; - if (P->is_latlong) - return 0; /* done in pj_inv3d() / pj_fwd3d() */ - - for (i = 0; i < n; i++) - if (z[dist*i] != HUGE_VAL ) - z[dist*i] *= fac; - - return 0; -} - - -/************************************************************************/ -/* pj_apply_vgridshift() */ -/* */ -/* This implementation takes uses the gridlist from a coordinate */ -/* system definition. If the gridlist has not yet been */ -/* populated in the coordinate system definition we set it up */ -/* now. */ -/************************************************************************/ -static int pj_apply_vgridshift( PJ *defn, - int inverse, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - if( defn->vgrids_legacy == nullptr ) - { - defn->vgrids_legacy = new ListOfVGrids; - auto vgrids = pj_vgrid_init(defn, "geoidgrids"); - if( vgrids.empty() ) - return 0; - *static_cast(defn->vgrids_legacy) = std::move(vgrids); - } - if( static_cast(defn->vgrids_legacy)->empty() ) - { - return 0; - } - - for( int i = 0; i < point_count; i++ ) - { - double value; - long io = i * point_offset; - PJ_LP input; - - input.phi = y[io]; - input.lam = x[io]; - - value = pj_vgrid_value(defn, *static_cast(defn->vgrids_legacy), input, 1.0); - - if( inverse ) - z[io] -= value; - else - z[io] += value; - - if( value == HUGE_VAL ) - { - std::string gridlist; - - proj_log_debug(defn, - "pj_apply_vgridshift(): failed to find a grid shift table for\n" - " location (%.7fdW,%.7fdN)", - x[io] * RAD_TO_DEG, - y[io] * RAD_TO_DEG ); - - for( const auto& gridset: *static_cast(defn->vgrids_legacy) ) - { - if( gridlist.empty() ) - gridlist += " tried: "; - else - gridlist += ','; - gridlist += gridset->name(); - } - - proj_log_debug(defn, "%s", gridlist.c_str()); - pj_ctx_set_errno( defn->ctx, PJD_ERR_GRID_AREA ); - - return PJD_ERR_GRID_AREA; - } - } - - return 0; -} - - -/* -------------------------------------------------------------------- */ -/* Transform to ellipsoidal heights if needed */ -/* -------------------------------------------------------------------- */ -static int geometric_to_orthometric (PJ *P, PJ_DIRECTION dir, long n, int dist, double *x, double *y, double *z) { - int err; - if (0==P->has_geoid_vgrids) - return 0; - if (z==nullptr) - return PJD_ERR_GEOCENTRIC; - err = pj_apply_vgridshift (P, dir==PJ_FWD ? 1 : 0, n, dist, x, y, z ); - if (err) - return pj_ctx_get_errno(P->ctx); - return 0; -} - - - -/* -------------------------------------------------------------------- */ -/* Convert datums if needed, and possible. */ -/* -------------------------------------------------------------------- */ -static int datum_transform (PJ *P, PJ *Q, long n, int dist, double *x, double *y, double *z) { - if (0==pj_datum_transform (P, Q, n, dist, x, y, z)) - return 0; - if (P->ctx->last_errno) - return P->ctx->last_errno; - return Q->ctx->last_errno; -} - - - - - -/* -------------------------------------------------------------------- */ -/* If a wrapping center other than 0 is provided, rewrap around */ -/* the suggested center (for latlong coordinate systems only). */ -/* -------------------------------------------------------------------- */ -static int long_wrap (PJ *P, long n, int dist, double *x) { - long i; - - /* Nothing to do? */ - if (P->is_geocent) - return 0; - if (!P->is_long_wrap_set) - return 0; - if (!P->is_latlong) - return 0; - - for (i = 0; i < n; i++ ) { - double val = x[dist*i]; - if (val == HUGE_VAL) - continue; - - /* Get fast in ] -2 PI, 2 PI [ range */ - val = fmod(val, M_TWOPI); - while( val < P->long_wrap_center - M_PI ) - val += M_TWOPI; - while( val > P->long_wrap_center + M_PI ) - val -= M_TWOPI; - x[dist*i] = val; - } - return 0; -} - - - -/************************************************************************/ -/* pj_transform() */ -/* */ -/* Currently this function doesn't recognise if two projections */ -/* are identical (to short circuit reprojection) because it is */ -/* difficult to compare PJ structures (since there are some */ -/* projection specific components). */ -/************************************************************************/ - -int pj_transform( - PJ *src, PJ *dst, - long point_count, int point_offset, - double *x, double *y, double *z -){ - int err; - - src->ctx->last_errno = 0; - dst->ctx->last_errno = 0; - - if( point_offset == 0 ) - point_offset = 1; - - /* Bring input to "normal form": longitude, latitude, ellipsoidal height */ - - err = adjust_axes (src, PJ_INV, point_count, point_offset, x, y, z); - if (err) - return err; - err = geographic_to_cartesian (src, PJ_INV, point_count, point_offset, x, y, z); - if (err) - return err; - err = projected_to_geographic (src, point_count, point_offset, x, y, z); - if (err) - return err; - err = prime_meridian (src, PJ_INV, point_count, point_offset, x); - if (err) - return err; - err = height_unit (src, PJ_INV, point_count, point_offset, z); - if (err) - return err; - err = geometric_to_orthometric (src, PJ_INV, point_count, point_offset, x, y, z); - if (err) - return err; - - /* At the center of the process we do the datum shift (if needed) */ - - err = datum_transform(src, dst, point_count, point_offset, x, y, z ); - if (err) - return err; - - /* Now get out on the other side: Bring "normal form" to output form */ - - err = geometric_to_orthometric (dst, PJ_FWD, point_count, point_offset, x, y, z); - if (err) - return err; - err = height_unit (dst, PJ_FWD, point_count, point_offset, z); - if (err) - return err; - err = prime_meridian (dst, PJ_FWD, point_count, point_offset, x); - if (err) - return err; - err = geographic_to_cartesian (dst, PJ_FWD, point_count, point_offset, x, y, z); - if (err) - return err; - err = geographic_to_projected (dst, point_count, point_offset, x, y, z); - if (err) - return err; - err = long_wrap (dst, point_count, point_offset, x); - if (err) - return err; - err = adjust_axes (dst, PJ_FWD, point_count, point_offset, x, y, z); - if (err) - return err; - - return 0; -} - - - -/************************************************************************/ -/* pj_geodetic_to_geocentric() */ -/************************************************************************/ - -int pj_geodetic_to_geocentric( double a, double es, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - double b; - int i; - GeocentricInfo gi; - int ret_errno = 0; - - if( es == 0.0 ) - b = a; - else - b = a * sqrt(1-es); - - if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 ) - { - return PJD_ERR_GEOCENTRIC; - } - - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - - if( x[io] == HUGE_VAL ) - continue; - - if( pj_Convert_Geodetic_To_Geocentric( &gi, y[io], x[io], z[io], - x+io, y+io, z+io ) != 0 ) - { - ret_errno = PJD_ERR_LAT_OR_LON_EXCEED_LIMIT; - x[io] = y[io] = HUGE_VAL; - /* but keep processing points! */ - } - } - - return ret_errno; -} - -/************************************************************************/ -/* pj_geocentric_to_geodetic() */ -/************************************************************************/ - -int pj_geocentric_to_geodetic( double a, double es, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - double b; - int i; - GeocentricInfo gi; - - if( es == 0.0 ) - b = a; - else - b = a * sqrt(1-es); - - if( pj_Set_Geocentric_Parameters( &gi, a, b ) != 0 ) - { - return PJD_ERR_GEOCENTRIC; - } - - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - - if( x[io] == HUGE_VAL ) - continue; - - pj_Convert_Geocentric_To_Geodetic( &gi, x[io], y[io], z[io], - y+io, x+io, z+io ); - } - - return 0; -} - -/************************************************************************/ -/* pj_compare_datums() */ -/* */ -/* Returns TRUE if the two datums are identical, otherwise */ -/* FALSE. */ -/************************************************************************/ - -int pj_compare_datums( PJ *srcdefn, PJ *dstdefn ) - -{ - if( srcdefn->datum_type != dstdefn->datum_type ) - { - return 0; - } - else if( srcdefn->a_orig != dstdefn->a_orig - || ABS(srcdefn->es_orig - dstdefn->es_orig) > 0.000000000050 ) - { - /* the tolerance for es is to ensure that GRS80 and WGS84 are - considered identical */ - return 0; - } - else if( srcdefn->datum_type == PJD_3PARAM ) - { - return (srcdefn->datum_params[0] == dstdefn->datum_params[0] - && srcdefn->datum_params[1] == dstdefn->datum_params[1] - && srcdefn->datum_params[2] == dstdefn->datum_params[2]); - } - else if( srcdefn->datum_type == PJD_7PARAM ) - { - return (srcdefn->datum_params[0] == dstdefn->datum_params[0] - && srcdefn->datum_params[1] == dstdefn->datum_params[1] - && srcdefn->datum_params[2] == dstdefn->datum_params[2] - && srcdefn->datum_params[3] == dstdefn->datum_params[3] - && srcdefn->datum_params[4] == dstdefn->datum_params[4] - && srcdefn->datum_params[5] == dstdefn->datum_params[5] - && srcdefn->datum_params[6] == dstdefn->datum_params[6]); - } - else if( srcdefn->datum_type == PJD_GRIDSHIFT ) - { - const char* srcnadgrids = - pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s; - const char* dstnadgrids = - pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s; - return srcnadgrids != nullptr && dstnadgrids != nullptr && - strcmp( srcnadgrids, dstnadgrids ) == 0; - } - else - return 1; -} - -/************************************************************************/ -/* pj_geocentic_to_wgs84() */ -/************************************************************************/ - -static -int pj_geocentric_to_wgs84( PJ *defn, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - int i; - - if( defn->datum_type == PJD_3PARAM ) - { - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - - if( x[io] == HUGE_VAL ) - continue; - - x[io] = x[io] + Dx_BF; - y[io] = y[io] + Dy_BF; - z[io] = z[io] + Dz_BF; - } - } - else if( defn->datum_type == PJD_7PARAM ) - { - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - double x_out, y_out, z_out; - - if( x[io] == HUGE_VAL ) - continue; - - x_out = M_BF*( x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF; - y_out = M_BF*( Rz_BF*x[io] + y[io] - Rx_BF*z[io]) + Dy_BF; - z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] + z[io]) + Dz_BF; - - x[io] = x_out; - y[io] = y_out; - z[io] = z_out; - } - } - - return 0; -} - -/************************************************************************/ -/* pj_geocentic_from_wgs84() */ -/************************************************************************/ - -static -int pj_geocentric_from_wgs84( PJ *defn, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - int i; - - if( defn->datum_type == PJD_3PARAM ) - { - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - - if( x[io] == HUGE_VAL ) - continue; - - x[io] = x[io] - Dx_BF; - y[io] = y[io] - Dy_BF; - z[io] = z[io] - Dz_BF; - } - } - else if( defn->datum_type == PJD_7PARAM ) - { - for( i = 0; i < point_count; i++ ) - { - long io = i * point_offset; - double x_tmp, y_tmp, z_tmp; - - if( x[io] == HUGE_VAL ) - continue; - - x_tmp = (x[io] - Dx_BF) / M_BF; - y_tmp = (y[io] - Dy_BF) / M_BF; - z_tmp = (z[io] - Dz_BF) / M_BF; - - x[io] = x_tmp + Rz_BF*y_tmp - Ry_BF*z_tmp; - y[io] = -Rz_BF*x_tmp + y_tmp + Rx_BF*z_tmp; - z[io] = Ry_BF*x_tmp - Rx_BF*y_tmp + z_tmp; - } - } - - return 0; -} - - -/************************************************************************/ -/* pj_apply_gridshift_2() */ -/* */ -/* This implementation uses the gridlist from a coordinate */ -/* system definition. If the gridlist has not yet been */ -/* populated in the coordinate system definition we set it up */ -/* now. */ -/************************************************************************/ -static -int pj_apply_gridshift_2( PJ *defn, int inverse, - long point_count, int point_offset, - double *x, double *y, double * /*z*/ ) - -{ - if( defn->hgrids_legacy == nullptr ) - { - defn->hgrids_legacy = new ListOfHGrids; - auto hgrids = pj_hgrid_init(defn, "nadgrids"); - if( hgrids.empty() ) - return 0; - *static_cast(defn->hgrids_legacy) = std::move(hgrids); - } - if( static_cast(defn->hgrids_legacy)->empty() ) - { - return 0; - } - - for( long i = 0; i < point_count; i++ ) - { - PJ_LP input; - - long io = i * point_offset; - input.phi = y[io]; - input.lam = x[io]; - - auto output = pj_hgrid_apply(defn->ctx, *static_cast(defn->hgrids_legacy), input, inverse ? PJ_INV : PJ_FWD); - - if ( output.lam != HUGE_VAL ) - { - y[io] = output.phi; - x[io] = output.lam; - } - else - { - if( defn->ctx->debug_level >= PJ_LOG_DEBUG_MAJOR ) - { - pj_log( defn->ctx, PJ_LOG_DEBUG_MAJOR, - "pj_apply_gridshift(): failed to find a grid shift table for\n" - " location (%.7fdW,%.7fdN)", - x[io] * RAD_TO_DEG, - y[io] * RAD_TO_DEG ); - } - } - } - - return 0; -} - - -/************************************************************************/ -/* pj_datum_transform() */ -/* */ -/* The input should be long/lat/z coordinates in radians in the */ -/* source datum, and the output should be long/lat/z */ -/* coordinates in radians in the destination datum. */ -/************************************************************************/ - -int pj_datum_transform( PJ *src, PJ *dst, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - double src_a, src_es, dst_a, dst_es; - int z_is_temp = FALSE; - -/* -------------------------------------------------------------------- */ -/* We cannot do any meaningful datum transformation if either */ -/* the source or destination are of an unknown datum type */ -/* (ie. only a +ellps declaration, no +datum). This is new */ -/* behavior for PROJ 4.6.0. */ -/* -------------------------------------------------------------------- */ - if( src->datum_type == PJD_UNKNOWN - || dst->datum_type == PJD_UNKNOWN ) - return 0; - -/* -------------------------------------------------------------------- */ -/* Short cut if the datums are identical. */ -/* -------------------------------------------------------------------- */ - if( pj_compare_datums( src, dst ) ) - return 0; - - src_a = src->a_orig; - src_es = src->es_orig; - - dst_a = dst->a_orig; - dst_es = dst->es_orig; - -/* -------------------------------------------------------------------- */ -/* Create a temporary Z array if one is not provided. */ -/* -------------------------------------------------------------------- */ - if( z == nullptr ) - { - size_t bytes = sizeof(double) * point_count * point_offset; - z = (double *) pj_malloc(bytes); - memset( z, 0, bytes ); - z_is_temp = TRUE; - } - -#define CHECK_RETURN(defn) {if( defn->ctx->last_errno != 0 && (defn->ctx->last_errno > 0 || get_transient_error_value(-defn->ctx->last_errno) == 0) ) { if( z_is_temp ) pj_dalloc(z); return defn->ctx->last_errno; }} - -/* -------------------------------------------------------------------- */ -/* If this datum requires grid shifts, then apply it to geodetic */ -/* coordinates. */ -/* -------------------------------------------------------------------- */ - if( src->datum_type == PJD_GRIDSHIFT ) - { - pj_apply_gridshift_2( src, 0, point_count, point_offset, x, y, z ); - CHECK_RETURN(src); - - src_a = SRS_WGS84_SEMIMAJOR; - src_es = SRS_WGS84_ESQUARED; - } - - if( dst->datum_type == PJD_GRIDSHIFT ) - { - dst_a = SRS_WGS84_SEMIMAJOR; - dst_es = SRS_WGS84_ESQUARED; - } - -/* ==================================================================== */ -/* Do we need to go through geocentric coordinates? */ -/* ==================================================================== */ - if( src_es != dst_es || src_a != dst_a - || src->datum_type == PJD_3PARAM - || src->datum_type == PJD_7PARAM - || dst->datum_type == PJD_3PARAM - || dst->datum_type == PJD_7PARAM) - { -/* -------------------------------------------------------------------- */ -/* Convert to geocentric coordinates. */ -/* -------------------------------------------------------------------- */ - src->ctx->last_errno = - pj_geodetic_to_geocentric( src_a, src_es, - point_count, point_offset, x, y, z ); - CHECK_RETURN(src); - -/* -------------------------------------------------------------------- */ -/* Convert between datums. */ -/* -------------------------------------------------------------------- */ - if( src->datum_type == PJD_3PARAM - || src->datum_type == PJD_7PARAM ) - { - pj_geocentric_to_wgs84( src, point_count, point_offset,x,y,z); - CHECK_RETURN(src); - } - - if( dst->datum_type == PJD_3PARAM - || dst->datum_type == PJD_7PARAM ) - { - pj_geocentric_from_wgs84( dst, point_count,point_offset,x,y,z); - CHECK_RETURN(dst); - } - -/* -------------------------------------------------------------------- */ -/* Convert back to geodetic coordinates. */ -/* -------------------------------------------------------------------- */ - dst->ctx->last_errno = - pj_geocentric_to_geodetic( dst_a, dst_es, - point_count, point_offset, x, y, z ); - CHECK_RETURN(dst); - } - -/* -------------------------------------------------------------------- */ -/* Apply grid shift to destination if required. */ -/* -------------------------------------------------------------------- */ - if( dst->datum_type == PJD_GRIDSHIFT ) - { - pj_apply_gridshift_2( dst, 1, point_count, point_offset, x, y, z ); - CHECK_RETURN(dst); - } - - if( z_is_temp ) - pj_dalloc( z ); - - return 0; -} - -/************************************************************************/ -/* adjust_axis() */ -/* */ -/* Normalize or de-normalized the x/y/z axes. The normal form */ -/* is "enu" (easting, northing, up). */ -/************************************************************************/ -static int adjust_axis( projCtx ctx, - const char *axis, int denormalize_flag, - long point_count, int point_offset, - double *x, double *y, double *z ) - -{ - double x_in, y_in, z_in = 0.0; - int i, i_axis; - - if( !denormalize_flag ) - { - for( i = 0; i < point_count; i++ ) - { - x_in = x[point_offset*i]; - y_in = y[point_offset*i]; - if( z ) - z_in = z[point_offset*i]; - - for( i_axis = 0; i_axis < 3; i_axis++ ) - { - double value; - - if( i_axis == 0 ) - value = x_in; - else if( i_axis == 1 ) - value = y_in; - else - value = z_in; - - switch( axis[i_axis] ) - { - case 'e': - x[point_offset*i] = value; - break; - case 'w': - x[point_offset*i] = -value; - break; - case 'n': - y[point_offset*i] = value; - break; - case 's': - y[point_offset*i] = -value; - break; - case 'u': - if( z ) - z[point_offset*i] = value; - break; - case 'd': - if( z ) - z[point_offset*i] = -value; - break; - default: - pj_ctx_set_errno( ctx, PJD_ERR_AXIS ); - return PJD_ERR_AXIS; - } - } /* i_axis */ - } /* i (point) */ - } - - else /* denormalize */ - { - for( i = 0; i < point_count; i++ ) - { - x_in = x[point_offset*i]; - y_in = y[point_offset*i]; - if( z ) - z_in = z[point_offset*i]; - - for( i_axis = 0; i_axis < 3; i_axis++ ) - { - double *target; - - if( i_axis == 2 && z == nullptr ) - continue; - - if( i_axis == 0 ) - target = x; - else if( i_axis == 1 ) - target = y; - else - target = z; - - switch( axis[i_axis] ) - { - case 'e': - target[point_offset*i] = x_in; break; - case 'w': - target[point_offset*i] = -x_in; break; - case 'n': - target[point_offset*i] = y_in; break; - case 's': - target[point_offset*i] = -y_in; break; - case 'u': - target[point_offset*i] = z_in; break; - case 'd': - target[point_offset*i] = -z_in; break; - default: - pj_ctx_set_errno( ctx, PJD_ERR_AXIS ); - return PJD_ERR_AXIS; - } - } /* i_axis */ - } /* i (point) */ - } - - return 0; -} -// --------------------------------------------------------------------------- - -void pj_deallocate_grids() -{ -} diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt index 2c0c19a9..3924f47d 100644 --- a/test/unit/CMakeLists.txt +++ b/test/unit/CMakeLists.txt @@ -69,17 +69,6 @@ set(PROJ_TEST_ENVIRONMENT "PROJ_SOURCE_DATA=${PROJ_SOURCE_DIR}/data" ) -add_executable(proj_pj_transform_test - main.cpp - pj_transform_test.cpp) -target_link_libraries(proj_pj_transform_test - GTest::gtest - ${PROJ_LIBRARIES}) -add_test(NAME proj_pj_transform_test COMMAND proj_pj_transform_test) -set_property(TEST proj_pj_transform_test - PROPERTY ENVIRONMENT ${PROJ_TEST_ENVIRONMENT}) - - add_executable(proj_errno_string_test main.cpp proj_errno_string_test.cpp) diff --git a/test/unit/Makefile.am b/test/unit/Makefile.am index 483cb0bd..4e931c2b 100644 --- a/test/unit/Makefile.am +++ b/test/unit/Makefile.am @@ -9,8 +9,7 @@ AM_CXXFLAGS = @CXX_WFLAGS@ @NO_ZERO_AS_NULL_POINTER_CONSTANT_FLAG@ PROJ_LIB ?= ../../data/for_tests -noinst_PROGRAMS = pj_transform_test -noinst_PROGRAMS += pj_phi2_test +noinst_PROGRAMS = pj_phi2_test noinst_PROGRAMS += proj_errno_string_test noinst_PROGRAMS += proj_angular_io_test noinst_PROGRAMS += proj_context_test @@ -21,11 +20,6 @@ noinst_PROGRAMS += test_network noinst_PROGRAMS += test_defmodel noinst_PROGRAMS += test_tinshift -pj_transform_test_SOURCES = pj_transform_test.cpp main.cpp -pj_transform_test_LDADD = ../../src/libproj.la @GTEST_LIBS@ - -pj_transform_test-check: pj_transform_test - PROJ_SKIP_READ_USER_WRITABLE_DIRECTORY=YES PROJ_LIB=$(PROJ_LIB) ./pj_transform_test pj_phi2_test_SOURCES = pj_phi2_test.cpp main.cpp pj_phi2_test_LDADD = ../../src/libproj.la @GTEST_LIBS@ @@ -84,6 +78,6 @@ test_tinshift_LDADD = ../../src/libproj.la @GTEST_LIBS@ test_tinshift-check: test_tinshift PROJ_LIB=$(PROJ_LIB) ./test_tinshift -check-local: pj_transform_test-check pj_phi2_test-check proj_errno_string_test-check \ +check-local: pj_phi2_test-check proj_errno_string_test-check \ proj_angular_io_test-check proj_context_test-check test_cpp_api-check \ gie_self_tests-check test_network-check test_defmodel-check test_tinshift-check diff --git a/test/unit/pj_transform_test.cpp b/test/unit/pj_transform_test.cpp deleted file mode 100644 index 6f02ec5a..00000000 --- a/test/unit/pj_transform_test.cpp +++ /dev/null @@ -1,737 +0,0 @@ -/****************************************************************************** - * - * Project: PROJ - * Purpose: Test pj_transform() legacy interface - * Author: Even Rouault - * - ****************************************************************************** - * Copyright (c) 2018, Even Rouault - * - * 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. - ****************************************************************************/ - -#include "gtest_include.h" -#include - -// PROJ include order is sensitive -// clang-format off -#include -#include "proj_internal.h" -// clang-format on - -namespace { - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_to_longlat) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_EQ(x, 2 * DEG_TO_RAD); - EXPECT_EQ(y, 49 * DEG_TO_RAD); - - x = 182 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_EQ(x, 182 * DEG_TO_RAD); - EXPECT_EQ(y, 49 * DEG_TO_RAD); - - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_to_proj) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=utm +zone=31 +datum=WGS84"); - double x = 3 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 500000, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_to_proj_tometer) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=utm +zone=31 +datum=WGS84 +to_meter=1000"); - double x = 3 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 500, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, proj_to_longlat) { - auto src = pj_init_plus("+proj=utm +zone=31 +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 500000; - double y = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 3 * DEG_TO_RAD, 1e-12); - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, proj_to_proj) { - auto src = pj_init_plus("+proj=utm +zone=31 +datum=WGS84"); - auto dst = pj_init_plus("+proj=utm +zone=31 +datum=WGS84"); - double x = 500000; - double y = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 500000, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_to_geocent) { - auto src = pj_init_plus("+proj=longlat +R=2"); - auto dst = pj_init_plus("+proj=geocent +R=2"); - double x = 0; - double y = 0; - double z = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - EXPECT_NEAR(z, 0, 1e-8); - - x = 90 * DEG_TO_RAD; - y = 0; - z = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 0, 1e-8); - EXPECT_NEAR(y, 2, 1e-8); - EXPECT_NEAR(z, 0, 1e-8); - - x = 0; - y = 90 * DEG_TO_RAD; - z = 0.1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 0, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - EXPECT_NEAR(z, 2 + 0.1, 1e-8); - - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_to_geocent_to_meter) { - auto src = pj_init_plus("+proj=longlat +R=2"); - auto dst = pj_init_plus("+proj=geocent +R=2 +to_meter=1000"); - double x = 0; - double y = 0; - double z = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2e-3, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - EXPECT_NEAR(z, 0, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, geocent_to_longlat) { - auto src = pj_init_plus("+proj=geocent +R=2"); - auto dst = pj_init_plus("+proj=longlat +R=2"); - double x = 0; - double y = 2; - double z = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 0, 1e-12); - - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, geocent_to_meter_to_longlat) { - auto src = pj_init_plus("+proj=geocent +to_meter=1000 +R=2"); - auto dst = pj_init_plus("+proj=longlat +R=2"); - double x = 0; - double y = 2e-3; - double z = 0; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 0, 1e-12); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, pm) { - auto src = pj_init_plus("+proj=longlat +pm=3 +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +pm=1 +datum=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, (2 + 3 - 1) * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_EQ(y, 49 * DEG_TO_RAD) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_geoc_to_longlat) { - auto src = pj_init_plus("+proj=longlat +geoc +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 48.809360314691766 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_to_longlat_geoc) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +geoc +datum=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 48.809360314691766 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, ellps_to_ellps_noop) { - auto src = pj_init_plus("+proj=longlat +ellps=clrk66"); - auto dst = pj_init_plus("+proj=longlat +ellps=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, towgs84_3param_noop) { - auto src = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3"); - auto dst = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 10; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, towgs84_7param_noop) { - auto src = - pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7"); - auto dst = - pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=1,2,3,4,5,6,7"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 10; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_towgs84_3param_to_datum) { - auto src = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 90 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - double z = 10; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10 + 1, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_towgs84_3param_to_datum_no_z) { - auto src = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 90 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_towgs84_7param_to_datum) { - auto src = - pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0,0,0,0,0.5"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 90 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - double z = 10; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 14.189073500223458, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, datum_to_longlat_towgs84_3param) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0"); - double x = 90 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - double z = 10 + 1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, datum_to_longlat_towgs84_7param) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = - pj_init_plus("+proj=longlat +ellps=WGS84 +towgs84=0,1,0,0,0,0,0.5"); - double x = 90 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - double z = 14.189073500223458; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 90 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, ellps_grs80_towgs84_to_datum_wgs84) { - auto src = pj_init_plus("+proj=longlat +ellps=GRS80 +towgs84=0,0,0"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 10; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-15) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-15) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, longlat_nadgrids_to_datum) { - auto src = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus"); - auto dst = pj_init_plus("+proj=longlat +datum=NAD83"); - double x = -100 * DEG_TO_RAD; - double y = 40 * DEG_TO_RAD; - double z = 10; - int ret = pj_transform(src, dst, 1, 0, &x, &y, &z); - EXPECT_TRUE(ret == 0 || ret == PJD_ERR_FAILED_TO_LOAD_GRID); - if (ret == 0) { - EXPECT_NEAR(x, -100.00040583667015 * DEG_TO_RAD, 1e-12) - << x / DEG_TO_RAD; - EXPECT_NEAR(y, 40.000005895651363 * DEG_TO_RAD, 1e-12) - << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10.000043224543333, 1e-8); - } - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, nadgrids_noop) { - auto src = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus"); - auto dst = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus"); - double x = -100 * DEG_TO_RAD; - double y = 40 * DEG_TO_RAD; - double z = 10; - int ret = pj_transform(src, dst, 1, 0, &x, &y, &z); - EXPECT_TRUE(ret == 0); - if (ret == 0) { - EXPECT_NEAR(x, -100 * DEG_TO_RAD, 1e-15) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 40 * DEG_TO_RAD, 1e-15) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - } - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, datum_to_longlat_nadgrids) { - auto src = pj_init_plus("+proj=longlat +datum=NAD83"); - auto dst = pj_init_plus("+proj=longlat +ellps=clrk66 +nadgrids=conus"); - double x = -100.00040583667015 * DEG_TO_RAD; - double y = 40.000005895651363 * DEG_TO_RAD; - double z = 10.000043224543333; - int ret = pj_transform(src, dst, 1, 0, &x, &y, &z); - EXPECT_TRUE(ret == 0 || ret == PJD_ERR_FAILED_TO_LOAD_GRID); - if (ret == 0) { - EXPECT_NEAR(x, -100 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 40 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 10, 1e-8); - } - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, long_wrap) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +lon_wrap=180"); - double x = -1 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 359 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 0 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, src_vto_meter) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84 +vto_meter=1000"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 1000, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, dest_vto_meter) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +vto_meter=1000"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 1000; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 1, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, src_axis_neu_with_z) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84 +axis=neu"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 49 * DEG_TO_RAD; - double y = 2 * DEG_TO_RAD; - double z = 1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 1, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, src_axis_neu_without_z) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84 +axis=neu"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 49 * DEG_TO_RAD; - double y = 2 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, src_axis_swd) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84 +axis=swd"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84"); - double x = 49 * DEG_TO_RAD; - double y = 2 * DEG_TO_RAD; - double z = -1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, -2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, -49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 1, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, dst_axis_neu) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +axis=neu"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, 49 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 2 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, 1, 1e-8); - pj_free(src); - pj_free(dst); -} -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, dst_axis_swd) { - auto src = pj_init_plus("+proj=longlat +datum=WGS84"); - auto dst = pj_init_plus("+proj=longlat +datum=WGS84 +axis=swd"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - double z = 1; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, &z), 0); - EXPECT_NEAR(x, -49 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, -2 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - EXPECT_NEAR(z, -1, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, init_epsg) { - auto src = pj_init_plus("+init=epsg:4326"); - ASSERT_TRUE(src != nullptr); - auto dst = pj_init_plus("+init=epsg:32631"); - ASSERT_TRUE(dst != nullptr); - double x = 3 * DEG_TO_RAD; - double y = 0 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 500000, 1e-8); - EXPECT_NEAR(y, 0, 1e-8); - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -TEST(proj_api_h, pj_set_searchpath) { - - const char *path = "/i_do/not/exit"; - pj_set_searchpath(1, &path); - { - auto info = proj_info(); - EXPECT_EQ(info.path_count, 1U); - ASSERT_NE(info.paths, nullptr); - ASSERT_NE(info.paths[0], nullptr); - EXPECT_EQ(std::string(info.paths[0]), path); - } - - pj_set_searchpath(0, nullptr); - { - auto info = proj_info(); - EXPECT_EQ(info.path_count, 0U); - EXPECT_EQ(info.paths, nullptr); - } -} - -// --------------------------------------------------------------------------- - -TEST(proj_api_h, pj_set_finder) { - - const auto myfinder = [](const char *) -> const char * { return nullptr; }; - pj_set_finder(myfinder); - - pj_set_finder(nullptr); -} - -// --------------------------------------------------------------------------- - -TEST(proj_api_h, default_fileapi) { - auto ctx = pj_ctx_alloc(); - auto fid = pj_open_lib(ctx, "proj.db", "rb"); - ASSERT_NE(fid, nullptr); - char header[6]; - ASSERT_EQ(pj_ctx_fread(ctx, header, 1, 6, fid), 6U); - ASSERT_TRUE(memcmp(header, "SQLite", 6) == 0); - ASSERT_EQ(pj_ctx_ftell(ctx, fid), 6); - ASSERT_EQ(pj_ctx_fseek(ctx, fid, 0, SEEK_SET), 0); - ASSERT_EQ(pj_ctx_ftell(ctx, fid), 0); - pj_ctx_fclose(ctx, fid); - pj_ctx_free(ctx); -} - -// --------------------------------------------------------------------------- - -TEST(pj_transform_test, ob_tran_to_meter_as_dest) { - auto src = pj_init_plus( - "+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs"); - auto dst = pj_init_plus("+ellps=WGS84 +proj=ob_tran +o_proj=latlon " - "+o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 " - "+to_meter=0.0174532925199433 +no_defs"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -// --------------------------------------------------------------------------- - -struct Spy { - bool gotInMyFOpen = false; - bool gotInMyFRead = false; - bool gotInMyFSeek = false; - bool gotInMyFTell = false; - bool gotInMyFClose = false; -}; - -struct MyFile { - FILE *fp; - Spy *spy; -}; - -static PAFile myFOpen(projCtx ctx, const char *filename, const char *access) { - FILE *fp = fopen(filename, access); - if (!fp) - return nullptr; - MyFile *myF = new MyFile; - myF->spy = (Spy *)pj_ctx_get_app_data(ctx); - myF->spy->gotInMyFOpen = true; - myF->fp = fp; - return reinterpret_cast(myF); -} - -static size_t myFRead(void *buffer, size_t size, size_t nmemb, PAFile file) { - MyFile *myF = reinterpret_cast(file); - myF->spy->gotInMyFRead = true; - return fread(buffer, size, nmemb, myF->fp); -} - -static int myFSeek(PAFile file, long offset, int whence) { - MyFile *myF = reinterpret_cast(file); - myF->spy->gotInMyFSeek = true; - return fseek(myF->fp, offset, whence); -} - -static long myFTell(PAFile file) { - MyFile *myF = reinterpret_cast(file); - myF->spy->gotInMyFTell = true; - return ftell(myF->fp); -} - -static void myFClose(PAFile file) { - MyFile *myF = reinterpret_cast(file); - myF->spy->gotInMyFClose = true; - fclose(myF->fp); - delete myF; -} - -TEST(proj_api_h, custom_fileapi) { - auto ctx = pj_ctx_alloc(); - Spy spy; - pj_ctx_set_app_data(ctx, &spy); - projFileAPI myAPI = {myFOpen, myFRead, myFSeek, myFTell, myFClose}; - pj_ctx_set_fileapi(ctx, &myAPI); - EXPECT_EQ(pj_ctx_get_fileapi(ctx), &myAPI); - auto fid = pj_open_lib(ctx, "proj.db", "rb"); - ASSERT_NE(fid, nullptr); - char header[6]; - ASSERT_EQ(pj_ctx_fread(ctx, header, 1, 6, fid), 6U); - ASSERT_TRUE(memcmp(header, "SQLite", 6) == 0); - ASSERT_EQ(pj_ctx_ftell(ctx, fid), 6); - ASSERT_EQ(pj_ctx_fseek(ctx, fid, 0, SEEK_SET), 0); - ASSERT_EQ(pj_ctx_ftell(ctx, fid), 0); - pj_ctx_fclose(ctx, fid); - pj_ctx_free(ctx); - EXPECT_TRUE(spy.gotInMyFOpen); - EXPECT_TRUE(spy.gotInMyFRead); - EXPECT_TRUE(spy.gotInMyFSeek); - EXPECT_TRUE(spy.gotInMyFTell); - EXPECT_TRUE(spy.gotInMyFClose); -} - -TEST(pj_transform_test, ob_tran_to_meter_as_srouce) { - auto src = pj_init_plus("+ellps=WGS84 +proj=ob_tran +o_proj=latlon " - "+o_lon_p=0.0 +o_lat_p=90.0 +lon_0=360.0 " - "+to_meter=0.0174532925199433 +no_defs"); - auto dst = pj_init_plus( - "+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs"); - double x = 2 * DEG_TO_RAD; - double y = 49 * DEG_TO_RAD; - EXPECT_EQ(pj_transform(src, dst, 1, 0, &x, &y, nullptr), 0); - EXPECT_NEAR(x, 2 * DEG_TO_RAD, 1e-12) << x / DEG_TO_RAD; - EXPECT_NEAR(y, 49 * DEG_TO_RAD, 1e-12) << y / DEG_TO_RAD; - pj_free(src); - pj_free(dst); -} - -} // namespace -- cgit v1.2.3