Line data Source code
1 : /******************************************************************************
2 : * $Id$
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 :
80 0 : 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 0 : double dfX, dfY, dfZ = 0.0;
88 :
89 0 : dfX = lon * DEG;
90 0 : dfY = lat * DEG;
91 :
92 0 : OCTTransform( hForCT, 1, &dfX, &dfY, &dfZ);
93 :
94 0 : *x = dfX;
95 0 : *y = dfY;
96 :
97 0 : return 0;
98 : }
99 :
100 : /******************************************************************************
101 : function to init a forward gctp transformer
102 : ******************************************************************************/
103 :
104 0 : 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 0 : OGRSpatialReferenceH hOutSourceSRS, hLatLong = NULL;
116 :
117 0 : *iflg = 0;
118 :
119 0 : hOutSourceSRS = OSRNewSpatialReference( NULL );
120 0 : OSRSetAxisMappingStrategy(hOutSourceSRS, OAMS_TRADITIONAL_GIS_ORDER);
121 :
122 0 : OSRImportFromUSGS( hOutSourceSRS, outsys, outzone, outparm, outdatum );
123 0 : hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84_LAT_LONG );
124 0 : OSRSetAxisMappingStrategy(hLatLong, OAMS_TRADITIONAL_GIS_ORDER);
125 :
126 0 : hForCT = OCTNewCoordinateTransformation( hLatLong, hOutSourceSRS );
127 :
128 0 : OSRDestroySpatialReference( hOutSourceSRS );
129 0 : OSRDestroySpatialReference( hLatLong );
130 :
131 0 : for_trans[outsys] = osr_for;
132 0 : }
133 :
134 : /******************************************************************************
135 : function for inverse gctp transformation
136 :
137 : gctp returns Longitude and Latitude values in radians
138 : ******************************************************************************/
139 :
140 0 : 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 0 : double dfX, dfY, dfZ = 0.0;
148 :
149 0 : dfX = x;
150 0 : dfY = y;
151 :
152 0 : OCTTransform( hInvCT, 1, &dfX, &dfY, &dfZ );
153 :
154 0 : *lon = dfX * RAD;
155 0 : *lat = dfY * RAD;
156 :
157 0 : return 0;
158 : }
159 :
160 : /******************************************************************************
161 : function to init a inverse gctp transformer
162 : ******************************************************************************/
163 :
164 0 : 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 0 : OGRSpatialReferenceH hInSourceSRS, hLatLong = NULL;
177 0 : *iflg = 0;
178 :
179 0 : hInSourceSRS = OSRNewSpatialReference( NULL );
180 0 : OSRSetAxisMappingStrategy(hInSourceSRS, OAMS_TRADITIONAL_GIS_ORDER);
181 0 : OSRImportFromUSGS( hInSourceSRS, insys, inzone, inparm, indatum );
182 :
183 0 : hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84_LAT_LONG );
184 0 : OSRSetAxisMappingStrategy(hLatLong, OAMS_TRADITIONAL_GIS_ORDER);
185 :
186 0 : hInvCT = OCTNewCoordinateTransformation( hInSourceSRS, hLatLong );
187 :
188 0 : OSRDestroySpatialReference( hInSourceSRS );
189 0 : OSRDestroySpatialReference( hLatLong );
190 :
191 0 : inv_trans[insys] = osr_inv;
192 0 : }
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
200 : void gctp_destroy(void) {
201 : OCTDestroyCoordinateTransformation ( hForCT );
202 : OCTDestroyCoordinateTransformation ( hInvCT );
203 : }
204 : #endif
|
