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 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 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 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 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 200 void gctp_destroy(void) { 201 OCTDestroyCoordinateTransformation ( hForCT ); 202 OCTDestroyCoordinateTransformation ( hInvCT ); 203 } 204 #endif 205