1 /******************************************************************************
2 * $Id: gctp_wrap.c 22751 2011-07-18 15:25:05Z winkey $
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 /***** static vars to store the transformers in *****/
50 /***** this is not thread safe *****/
51
52 static OGRCoordinateTransformationH hForCT, hInvCT;
53
54 /******************************************************************************
55 function for forward gctp transformation
56
57 gctp expects Longitude and Latitude values to be in radians
58 ******************************************************************************/
59
osr_for(double lon,double lat,double * x,double * y)60 int32 osr_for(
61 double lon, /* (I) Longitude */
62 double lat, /* (I) Latitude */
63 double *x, /* (O) X projection coordinate */
64 double *y) /* (O) Y projection coordinate */
65 {
66
67 double dfX, dfY, dfZ = 0.0;
68
69 dfX = lon * DEG;
70 dfY = lat * DEG;
71
72 OCTTransform( hForCT, 1, &dfX, &dfY, &dfZ);
73
74 *x = dfX;
75 *y = dfY;
76
77 return 0;
78 }
79
80 /******************************************************************************
81 function to init a forward gctp transformer
82 ******************************************************************************/
83
for_init(int32 outsys,int32 outzone,float64 * outparm,int32 outdatum,char * fn27,char * fn83,int32 * iflg,int32 (* for_trans[])(double,double,double *,double *))84 void for_init(
85 int32 outsys, /* output system code */
86 int32 outzone, /* output zone number */
87 float64 *outparm, /* output array of projection parameters */
88 int32 outdatum, /* output datum */
89 char *fn27, /* NAD 1927 parameter file */
90 char *fn83, /* NAD 1983 parameter file */
91 int32 *iflg, /* status flag */
92 int32 (*for_trans[])(double, double, double *, double *))
93 /* forward function pointer */
94 {
95 OGRSpatialReferenceH hOutSourceSRS, hLatLong = NULL;
96
97 *iflg = 0;
98
99 hOutSourceSRS = OSRNewSpatialReference( NULL );
100 OSRImportFromUSGS( hOutSourceSRS, outsys, outzone, outparm, outdatum );
101 hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84 );
102
103 hForCT = OCTNewCoordinateTransformation( hLatLong, hOutSourceSRS );
104
105 OSRDestroySpatialReference( hOutSourceSRS );
106 OSRDestroySpatialReference( hLatLong );
107
108 for_trans[outsys] = osr_for;
109 }
110
111 /******************************************************************************
112 function for inverse gctp transformation
113
114 gctp returns Longitude and Latitude values in radians
115 ******************************************************************************/
116
osr_inv(double x,double y,double * lon,double * lat)117 int32 osr_inv(
118 double x, /* (I) X projection coordinate */
119 double y, /* (I) Y projection coordinate */
120 double *lon, /* (O) Longitude */
121 double *lat) /* (O) Latitude */
122 {
123
124 double dfX, dfY, dfZ = 0.0;
125
126 dfX = x;
127 dfY = y;
128
129 OCTTransform( hInvCT, 1, &dfX, &dfY, &dfZ );
130
131 *lon = dfX * RAD;
132 *lat = dfY * RAD;
133
134 return 0;
135 }
136
137 /******************************************************************************
138 function to init a inverse gctp transformer
139 ******************************************************************************/
140
inv_init(int32 insys,int32 inzone,float64 * inparm,int32 indatum,char * fn27,char * fn83,int32 * iflg,int32 (* inv_trans[])(double,double,double *,double *))141 void inv_init(
142 int32 insys, /* input system code */
143 int32 inzone, /* input zone number */
144 float64 *inparm, /* input array of projection parameters */
145 int32 indatum, /* input datum code */
146 char *fn27, /* NAD 1927 parameter file */
147 char *fn83, /* NAD 1983 parameter file */
148 int32 *iflg, /* status flag */
149 int32 (*inv_trans[])(double, double, double*, double*))
150 /* inverse function pointer */
151 {
152
153 OGRSpatialReferenceH hInSourceSRS, hLatLong = NULL;
154 *iflg = 0;
155
156 hInSourceSRS = OSRNewSpatialReference( NULL );
157 OSRImportFromUSGS( hInSourceSRS, insys, inzone, inparm, indatum );
158
159 hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84 );
160
161 hInvCT = OCTNewCoordinateTransformation( hInSourceSRS, hLatLong );
162
163 OSRDestroySpatialReference( hInSourceSRS );
164 OSRDestroySpatialReference( hLatLong );
165
166 inv_trans[insys] = osr_inv;
167 }
168
169 /******************************************************************************
170 function to cleanup the transformers
171
172 note: gctp does not have a function that does this
173 ******************************************************************************/
174
gctp_destroy(void)175 void gctp_destroy(void) {
176 OCTDestroyCoordinateTransformation ( hForCT );
177 OCTDestroyCoordinateTransformation ( hInvCT );
178 }
179