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 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "ogr_srs_api.h"
16 : #include <stdlib.h>
17 :
18 : #include "mfhdf.h"
19 :
20 : #include <math.h>
21 :
22 : #ifndef PI
23 : #ifndef M_PI
24 : #define PI (3.141592653589793238)
25 : #else
26 : #define PI (M_PI)
27 : #endif
28 : #endif
29 :
30 : #define DEG (180.0 / PI)
31 : #define RAD (PI / 180.0)
32 :
33 : void for_init(
34 : int32 outsys, /* output system code */
35 : int32 outzone, /* output zone number */
36 : float64 *outparm, /* output array of projection parameters */
37 : int32 outdatum, /* output datum */
38 : char *fn27, /* NAD 1927 parameter file */
39 : char *fn83, /* NAD 1983 parameter file */
40 : int32 *iflg, /* status flag */
41 : int32 (*for_trans[])(double, double, double *, double *));
42 :
43 : void inv_init(
44 : int32 insys, /* input system code */
45 : int32 inzone, /* input zone number */
46 : float64 *inparm, /* input array of projection parameters */
47 : int32 indatum, /* input datum code */
48 : char *fn27, /* NAD 1927 parameter file */
49 : char *fn83, /* NAD 1983 parameter file */
50 : int32 *iflg, /* status flag */
51 : int32 (*inv_trans[])(double, double, double*, double*));
52 :
53 : /***** static vars to store the transformers in *****/
54 : /***** this is not thread safe *****/
55 :
56 : static OGRCoordinateTransformationH hForCT, hInvCT;
57 :
58 : /******************************************************************************
59 : function for forward gctp transformation
60 :
61 : gctp expects Longitude and Latitude values to be in radians
62 : ******************************************************************************/
63 :
64 0 : static int32 osr_for(
65 : double lon, /* (I) Longitude */
66 : double lat, /* (I) Latitude */
67 : double *x, /* (O) X projection coordinate */
68 : double *y) /* (O) Y projection coordinate */
69 : {
70 :
71 0 : double dfX, dfY, dfZ = 0.0;
72 :
73 0 : dfX = lon * DEG;
74 0 : dfY = lat * DEG;
75 :
76 0 : OCTTransform( hForCT, 1, &dfX, &dfY, &dfZ);
77 :
78 0 : *x = dfX;
79 0 : *y = dfY;
80 :
81 0 : return 0;
82 : }
83 :
84 : /******************************************************************************
85 : function to init a forward gctp transformer
86 : ******************************************************************************/
87 :
88 0 : void for_init(
89 : int32 outsys, /* output system code */
90 : int32 outzone, /* output zone number */
91 : float64 *outparm, /* output array of projection parameters */
92 : int32 outdatum, /* output datum */
93 : CPL_UNUSED char *fn27, /* NAD 1927 parameter file */
94 : CPL_UNUSED char *fn83, /* NAD 1983 parameter file */
95 : int32 *iflg, /* status flag */
96 : int32 (*for_trans[])(double, double, double *, double *))
97 : /* forward function pointer */
98 : {
99 0 : OGRSpatialReferenceH hOutSourceSRS, hLatLong = NULL;
100 :
101 0 : *iflg = 0;
102 :
103 0 : hOutSourceSRS = OSRNewSpatialReference( NULL );
104 0 : OSRSetAxisMappingStrategy(hOutSourceSRS, OAMS_TRADITIONAL_GIS_ORDER);
105 :
106 0 : OSRImportFromUSGS( hOutSourceSRS, outsys, outzone, outparm, outdatum );
107 0 : hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84_LAT_LONG );
108 0 : OSRSetAxisMappingStrategy(hLatLong, OAMS_TRADITIONAL_GIS_ORDER);
109 :
110 0 : hForCT = OCTNewCoordinateTransformation( hLatLong, hOutSourceSRS );
111 :
112 0 : OSRDestroySpatialReference( hOutSourceSRS );
113 0 : OSRDestroySpatialReference( hLatLong );
114 :
115 0 : for_trans[outsys] = osr_for;
116 0 : }
117 :
118 : /******************************************************************************
119 : function for inverse gctp transformation
120 :
121 : gctp returns Longitude and Latitude values in radians
122 : ******************************************************************************/
123 :
124 0 : static int32 osr_inv(
125 : double x, /* (I) X projection coordinate */
126 : double y, /* (I) Y projection coordinate */
127 : double *lon, /* (O) Longitude */
128 : double *lat) /* (O) Latitude */
129 : {
130 :
131 0 : double dfX, dfY, dfZ = 0.0;
132 :
133 0 : dfX = x;
134 0 : dfY = y;
135 :
136 0 : OCTTransform( hInvCT, 1, &dfX, &dfY, &dfZ );
137 :
138 0 : *lon = dfX * RAD;
139 0 : *lat = dfY * RAD;
140 :
141 0 : return 0;
142 : }
143 :
144 : /******************************************************************************
145 : function to init a inverse gctp transformer
146 : ******************************************************************************/
147 :
148 0 : void inv_init(
149 : int32 insys, /* input system code */
150 : int32 inzone, /* input zone number */
151 : float64 *inparm, /* input array of projection parameters */
152 : int32 indatum, /* input datum code */
153 : CPL_UNUSED char *fn27, /* NAD 1927 parameter file */
154 : CPL_UNUSED char *fn83, /* NAD 1983 parameter file */
155 : int32 *iflg, /* status flag */
156 : int32 (*inv_trans[])(double, double, double*, double*))
157 : /* inverse function pointer */
158 : {
159 :
160 0 : OGRSpatialReferenceH hInSourceSRS, hLatLong = NULL;
161 0 : *iflg = 0;
162 :
163 0 : hInSourceSRS = OSRNewSpatialReference( NULL );
164 0 : OSRSetAxisMappingStrategy(hInSourceSRS, OAMS_TRADITIONAL_GIS_ORDER);
165 0 : OSRImportFromUSGS( hInSourceSRS, insys, inzone, inparm, indatum );
166 :
167 0 : hLatLong = OSRNewSpatialReference ( SRS_WKT_WGS84_LAT_LONG );
168 0 : OSRSetAxisMappingStrategy(hLatLong, OAMS_TRADITIONAL_GIS_ORDER);
169 :
170 0 : hInvCT = OCTNewCoordinateTransformation( hInSourceSRS, hLatLong );
171 :
172 0 : OSRDestroySpatialReference( hInSourceSRS );
173 0 : OSRDestroySpatialReference( hLatLong );
174 :
175 0 : inv_trans[insys] = osr_inv;
176 0 : }
177 :
178 : /******************************************************************************
179 : function to cleanup the transformers
180 :
181 : note: gctp does not have a function that does this
182 : ******************************************************************************/
183 : #ifndef GDAL_COMPILATION
184 : void gctp_destroy(void) {
185 : OCTDestroyCoordinateTransformation ( hForCT );
186 : OCTDestroyCoordinateTransformation ( hInvCT );
187 : }
188 : #endif
|
