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