1 /******************************************************************************
2 * $Id: gctp_wrap.c 8e5eeb35bf76390e3134a4ea7076dab7d478ea0e 2018-11-14 22:55:13 +0100 Even Rouault $
3 *
4 * Project: Hierarchical Data Format Release 4 (HDF4)
5 * Purpose: This is the wrapper code to use OGR Coordinate Transformation
6 * services instead of GCTP library
7 * Author: Andrey Kiselev, dron@ak4719.spb.edu
8 *
9 ******************************************************************************
10 * Copyright (c) 2004, Andrey Kiselev <dron@ak4719.spb.edu>
11 *
12 * Permission is hereby granted, free of charge, to any person obtaining a
13 * copy of this software and associated documentation files (the "Software"),
14 * to deal in the Software without restriction, including without limitation
15 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons to whom the
17 * Software is furnished to do so, subject to the following conditions:
18 *
19 * The above copyright notice and this permission notice shall be included
20 * in all copies or substantial portions of the Software.
21 *
22 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 * DEALINGS IN THE SOFTWARE.
29 ****************************************************************************/
30
31 #include "ogr_srs_api.h"
32 #include <stdlib.h>
33
34 #include "mfhdf.h"
35
36 #include <math.h>
37
38 #ifndef PI
39 #ifndef M_PI
40 #define PI (3.141592653589793238)
41 #else
42 #define PI (M_PI)
43 #endif
44 #endif
45
46 #define DEG (180.0 / PI)
47 #define RAD (PI / 180.0)
48
49 void for_init(
50 int32 outsys, /* output system code */
51 int32 outzone, /* output zone number */
52 float64 *outparm, /* output array of projection parameters */
53 int32 outdatum, /* output datum */
54 char *fn27, /* NAD 1927 parameter file */
55 char *fn83, /* NAD 1983 parameter file */
56 int32 *iflg, /* status flag */
57 int32 (*for_trans[])(double, double, double *, double *));
58
59 void inv_init(
60 int32 insys, /* input system code */
61 int32 inzone, /* input zone number */
62 float64 *inparm, /* input array of projection parameters */
63 int32 indatum, /* input datum code */
64 char *fn27, /* NAD 1927 parameter file */
65 char *fn83, /* NAD 1983 parameter file */
66 int32 *iflg, /* status flag */
67 int32 (*inv_trans[])(double, double, double*, double*));
68
69 /***** static vars to store the transformers in *****/
70 /***** this is not thread safe *****/
71
72 static OGRCoordinateTransformationH hForCT, hInvCT;
73
74 /******************************************************************************
75 function for forward gctp transformation
76
77 gctp expects Longitude and Latitude values to be in radians
78 ******************************************************************************/
79
osr_for(double lon,double lat,double * x,double * y)80 static int32 osr_for(
81 double lon, /* (I) Longitude */
82 double lat, /* (I) Latitude */
83 double *x, /* (O) X projection coordinate */
84 double *y) /* (O) Y projection coordinate */
85 {
86
87 double dfX, dfY, dfZ = 0.0;
88
89 dfX = lon * DEG;
90 dfY = lat * DEG;
91
92 OCTTransform( hForCT, 1, &dfX, &dfY, &dfZ);
93
94 *x = dfX;
95 *y = dfY;
96
97 return 0;
98 }
99
100 /******************************************************************************
101 function to init a forward gctp transformer
102 ******************************************************************************/
103
for_init(int32 outsys,int32 outzone,float64 * outparm,int32 outdatum,CPL_UNUSED char * fn27,CPL_UNUSED char * fn83,int32 * iflg,int32 (* for_trans[])(double,double,double *,double *))104 void for_init(
105 int32 outsys, /* output system code */
106 int32 outzone, /* output zone number */
107 float64 *outparm, /* output array of projection parameters */
108 int32 outdatum, /* output datum */
109 CPL_UNUSED char *fn27, /* NAD 1927 parameter file */
110 CPL_UNUSED char *fn83, /* NAD 1983 parameter file */
111 int32 *iflg, /* status flag */
112 int32 (*for_trans[])(double, double, double *, double *))
113 /* forward function pointer */
114 {
115 OGRSpatialReferenceH hOutSourceSRS, hLatLong = NULL;
116
117 *iflg = 0;
118
119 hOutSourceSRS = OSRNewSpatialReference( NULL );
120 OSRSetAxisMappingStrategy(hOutSourceSRS, OAMS_TRADITIONAL_GIS_ORDER);
121
122 OSRImportFromUSGS( hOutSourceSRS, outsys, outzone, outparm, outdatum );
123 hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84_LAT_LONG );
124 OSRSetAxisMappingStrategy(hLatLong, OAMS_TRADITIONAL_GIS_ORDER);
125
126 hForCT = OCTNewCoordinateTransformation( hLatLong, hOutSourceSRS );
127
128 OSRDestroySpatialReference( hOutSourceSRS );
129 OSRDestroySpatialReference( hLatLong );
130
131 for_trans[outsys] = osr_for;
132 }
133
134 /******************************************************************************
135 function for inverse gctp transformation
136
137 gctp returns Longitude and Latitude values in radians
138 ******************************************************************************/
139
osr_inv(double x,double y,double * lon,double * lat)140 static int32 osr_inv(
141 double x, /* (I) X projection coordinate */
142 double y, /* (I) Y projection coordinate */
143 double *lon, /* (O) Longitude */
144 double *lat) /* (O) Latitude */
145 {
146
147 double dfX, dfY, dfZ = 0.0;
148
149 dfX = x;
150 dfY = y;
151
152 OCTTransform( hInvCT, 1, &dfX, &dfY, &dfZ );
153
154 *lon = dfX * RAD;
155 *lat = dfY * RAD;
156
157 return 0;
158 }
159
160 /******************************************************************************
161 function to init a inverse gctp transformer
162 ******************************************************************************/
163
inv_init(int32 insys,int32 inzone,float64 * inparm,int32 indatum,CPL_UNUSED char * fn27,CPL_UNUSED char * fn83,int32 * iflg,int32 (* inv_trans[])(double,double,double *,double *))164 void inv_init(
165 int32 insys, /* input system code */
166 int32 inzone, /* input zone number */
167 float64 *inparm, /* input array of projection parameters */
168 int32 indatum, /* input datum code */
169 CPL_UNUSED char *fn27, /* NAD 1927 parameter file */
170 CPL_UNUSED char *fn83, /* NAD 1983 parameter file */
171 int32 *iflg, /* status flag */
172 int32 (*inv_trans[])(double, double, double*, double*))
173 /* inverse function pointer */
174 {
175
176 OGRSpatialReferenceH hInSourceSRS, hLatLong = NULL;
177 *iflg = 0;
178
179 hInSourceSRS = OSRNewSpatialReference( NULL );
180 OSRSetAxisMappingStrategy(hInSourceSRS, OAMS_TRADITIONAL_GIS_ORDER);
181 OSRImportFromUSGS( hInSourceSRS, insys, inzone, inparm, indatum );
182
183 hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84_LAT_LONG );
184 OSRSetAxisMappingStrategy(hLatLong, OAMS_TRADITIONAL_GIS_ORDER);
185
186 hInvCT = OCTNewCoordinateTransformation( hInSourceSRS, hLatLong );
187
188 OSRDestroySpatialReference( hInSourceSRS );
189 OSRDestroySpatialReference( hLatLong );
190
191 inv_trans[insys] = osr_inv;
192 }
193
194 /******************************************************************************
195 function to cleanup the transformers
196
197 note: gctp does not have a function that does this
198 ******************************************************************************/
199 #ifndef GDAL_COMPILATION
gctp_destroy(void)200 void gctp_destroy(void) {
201 OCTDestroyCoordinateTransformation ( hForCT );
202 OCTDestroyCoordinateTransformation ( hInvCT );
203 }
204 #endif
205