1 /******************************************************************************
2 *
3 * Project: OpenGIS Simple Features Reference Implementation
4 * Purpose: OGRSpatialReference translation from OziExplorer
5 * georeferencing information.
6 * Author: Andrey Kiselev, dron@ak4719.spb.edu
7 *
8 ******************************************************************************
9 * Copyright (c) 2009, Andrey Kiselev <dron@ak4719.spb.edu>
10 * Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
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 "cpl_port.h"
32 #include "ogr_spatialref.h"
33
34 #include <cstdlib>
35 #include <cstring>
36
37 #include "cpl_conv.h"
38 #include "cpl_csv.h"
39 #include "cpl_error.h"
40 #include "cpl_string.h"
41 #include "ogr_core.h"
42 #include "ogr_srs_api.h"
43
44 CPL_CVSID("$Id: ogr_srs_ozi.cpp fa752ad6eabafaf630a704e1892a9d837d683cb3 2021-03-06 17:04:38 +0100 Even Rouault $")
45
46 /************************************************************************/
47 /* OSRImportFromOzi() */
48 /************************************************************************/
49
50 /**
51 * Import coordinate system from OziExplorer projection definition.
52 *
53 * This function will import projection definition in style, used by
54 * OziExplorer software.
55 *
56 * Note: another version of this function with a different signature existed
57 * in GDAL 1.X.
58 *
59 * @param hSRS spatial reference object.
60 * @param papszLines Map file lines. This is an array of strings containing
61 * the whole OziExplorer .MAP file. The array is terminated by a NULL pointer.
62 *
63 * @return OGRERR_NONE on success or an error code in case of failure.
64 *
65 * @since OGR 2.0
66 */
67
OSRImportFromOzi(OGRSpatialReferenceH hSRS,const char * const * papszLines)68 OGRErr OSRImportFromOzi( OGRSpatialReferenceH hSRS,
69 const char * const* papszLines )
70
71 {
72 VALIDATE_POINTER1( hSRS, "OSRImportFromOzi", OGRERR_FAILURE );
73
74 return OGRSpatialReference::FromHandle(hSRS)->importFromOzi( papszLines );
75 }
76
77 /************************************************************************/
78 /* importFromOzi() */
79 /************************************************************************/
80
81 /**
82 * Import coordinate system from OziExplorer projection definition.
83 *
84 * This method will import projection definition in style, used by
85 * OziExplorer software.
86 *
87 * @param papszLines Map file lines. This is an array of strings containing
88 * the whole OziExplorer .MAP file. The array is terminated by a NULL pointer.
89 *
90 * @return OGRERR_NONE on success or an error code in case of failure.
91 *
92 * @since OGR 1.10
93 */
94
importFromOzi(const char * const * papszLines)95 OGRErr OGRSpatialReference::importFromOzi( const char * const* papszLines )
96 {
97 const char *pszDatum;
98 const char *pszProj = nullptr;
99 const char *pszProjParams = nullptr;
100
101 Clear();
102
103 const int nLines = CSLCount(papszLines);
104 if( nLines < 5 )
105 return OGRERR_NOT_ENOUGH_DATA;
106
107 pszDatum = papszLines[4];
108
109 for( int iLine = 5; iLine < nLines; iLine++ )
110 {
111 if( STARTS_WITH_CI(papszLines[iLine], "Map Projection") )
112 {
113 pszProj = papszLines[iLine];
114 }
115 else if( STARTS_WITH_CI(papszLines[iLine], "Projection Setup") )
116 {
117 pszProjParams = papszLines[iLine];
118 }
119 }
120
121 if( !(pszDatum && pszProj && pszProjParams) )
122 return OGRERR_NOT_ENOUGH_DATA;
123
124 /* -------------------------------------------------------------------- */
125 /* Operate on the basis of the projection name. */
126 /* -------------------------------------------------------------------- */
127 char **papszProj = CSLTokenizeStringComplex( pszProj, ",", TRUE, TRUE );
128 char **papszProjParams = CSLTokenizeStringComplex( pszProjParams, ",",
129 TRUE, TRUE );
130 char **papszDatum = nullptr;
131
132 if( CSLCount(papszProj) < 2 )
133 {
134 goto not_enough_data;
135 }
136
137 if( STARTS_WITH_CI(papszProj[1], "Latitude/Longitude") )
138 {
139 // Do nothing.
140 }
141 else if( STARTS_WITH_CI(papszProj[1], "Mercator") )
142 {
143 if( CSLCount(papszProjParams) < 6 ) goto not_enough_data;
144 double dfScale = CPLAtof(papszProjParams[3]);
145 // If unset, default to scale = 1.
146 if( papszProjParams[3][0] == 0 ) dfScale = 1;
147 SetMercator( CPLAtof(papszProjParams[1]), CPLAtof(papszProjParams[2]),
148 dfScale,
149 CPLAtof(papszProjParams[4]), CPLAtof(papszProjParams[5]) );
150 }
151 else if( STARTS_WITH_CI(papszProj[1], "Transverse Mercator") )
152 {
153 if( CSLCount(papszProjParams) < 6 ) goto not_enough_data;
154 SetTM( CPLAtof(papszProjParams[1]), CPLAtof(papszProjParams[2]),
155 CPLAtof(papszProjParams[3]),
156 CPLAtof(papszProjParams[4]), CPLAtof(papszProjParams[5]) );
157 }
158 else if( STARTS_WITH_CI(papszProj[1], "Lambert Conformal Conic") )
159 {
160 if( CSLCount(papszProjParams) < 8 ) goto not_enough_data;
161 SetLCC( CPLAtof(papszProjParams[6]), CPLAtof(papszProjParams[7]),
162 CPLAtof(papszProjParams[1]), CPLAtof(papszProjParams[2]),
163 CPLAtof(papszProjParams[4]), CPLAtof(papszProjParams[5]) );
164 }
165 else if( STARTS_WITH_CI(papszProj[1], "Sinusoidal") )
166 {
167 if( CSLCount(papszProjParams) < 6 ) goto not_enough_data;
168 SetSinusoidal( CPLAtof(papszProjParams[2]),
169 CPLAtof(papszProjParams[4]), CPLAtof(papszProjParams[5]) );
170 }
171 else if( STARTS_WITH_CI(papszProj[1], "Albers Equal Area") )
172 {
173 if( CSLCount(papszProjParams) < 8 ) goto not_enough_data;
174 SetACEA( CPLAtof(papszProjParams[6]), CPLAtof(papszProjParams[7]),
175 CPLAtof(papszProjParams[1]), CPLAtof(papszProjParams[2]),
176 CPLAtof(papszProjParams[4]), CPLAtof(papszProjParams[5]) );
177 }
178 else if( STARTS_WITH_CI(
179 papszProj[1], "(UTM) Universal Transverse Mercator") &&
180 nLines > 5 )
181 {
182 // Look for the UTM zone in the calibration point data.
183 int iLine = 5; // Used after for.
184 for( ; iLine < nLines; iLine++ )
185 {
186 if( STARTS_WITH_CI(papszLines[iLine], "Point") )
187 {
188 char **papszTok =
189 CSLTokenizeString2(papszLines[iLine], ",",
190 CSLT_ALLOWEMPTYTOKENS
191 | CSLT_STRIPLEADSPACES
192 | CSLT_STRIPENDSPACES);
193 if( CSLCount(papszTok) < 17
194 || EQUAL(papszTok[2], "")
195 || EQUAL(papszTok[13], "")
196 || EQUAL(papszTok[14], "")
197 || EQUAL(papszTok[15], "")
198 || EQUAL(papszTok[16], "") )
199 {
200 CSLDestroy(papszTok);
201 continue;
202 }
203 SetUTM( atoi(papszTok[13]), EQUAL(papszTok[16], "N") );
204 CSLDestroy(papszTok);
205 break;
206 }
207 }
208 if( iLine == nLines ) // Try to guess the UTM zone.
209 {
210 float fMinLongitude = 1000.0f;
211 float fMaxLongitude = -1000.0f;
212 float fMinLatitude = 1000.0f;
213 float fMaxLatitude = -1000.0f;
214 bool bFoundMMPLL = false;
215 for( iLine = 5; iLine < nLines; iLine++ )
216 {
217 if( STARTS_WITH_CI(papszLines[iLine], "MMPLL") )
218 {
219 char **papszTok =
220 CSLTokenizeString2(papszLines[iLine], ",",
221 CSLT_ALLOWEMPTYTOKENS
222 | CSLT_STRIPLEADSPACES
223 | CSLT_STRIPENDSPACES);
224 if( CSLCount(papszTok) < 4 )
225 {
226 CSLDestroy(papszTok);
227 continue;
228 }
229 const float fLongitude =
230 static_cast<float>(CPLAtofM(papszTok[2]));
231 const float fLatitude =
232 static_cast<float>(CPLAtofM(papszTok[3]));
233 CSLDestroy(papszTok);
234
235 bFoundMMPLL = true;
236
237 if( fMinLongitude > fLongitude )
238 fMinLongitude = fLongitude;
239 if( fMaxLongitude < fLongitude )
240 fMaxLongitude = fLongitude;
241 if( fMinLatitude > fLatitude )
242 fMinLatitude = fLatitude;
243 if( fMaxLatitude < fLatitude )
244 fMaxLatitude = fLatitude;
245 }
246 }
247 const float fMedianLatitude = (fMinLatitude + fMaxLatitude) / 2;
248 const float fMedianLongitude = (fMinLongitude + fMaxLongitude) / 2;
249 if( bFoundMMPLL && fMaxLatitude <= 90 )
250 {
251 int nUtmZone = 0;
252 if( fMedianLatitude >= 56 && fMedianLatitude <= 64 &&
253 fMedianLongitude >= 3 && fMedianLongitude <= 12 )
254 nUtmZone = 32; // Norway exception.
255 else if( fMedianLatitude >= 72 && fMedianLatitude <= 84 &&
256 fMedianLongitude >= 0 && fMedianLongitude <= 42 )
257 // Svalbard exception.
258 nUtmZone =
259 static_cast<int>((fMedianLongitude + 3) / 12) * 2 + 31;
260 else
261 nUtmZone =
262 static_cast<int>((fMedianLongitude + 180 ) / 6) + 1;
263 SetUTM( nUtmZone, fMedianLatitude >= 0 );
264 }
265 else
266 {
267 CPLDebug( "OSR_Ozi", "UTM Zone not found");
268 }
269 }
270 }
271 else if( STARTS_WITH_CI(papszProj[1], "(I) France Zone I") )
272 {
273 SetLCC1SP( 49.5, 2.337229167, 0.99987734, 600000, 1200000 );
274 }
275 else if( STARTS_WITH_CI(papszProj[1], "(II) France Zone II") )
276 {
277 SetLCC1SP( 46.8, 2.337229167, 0.99987742, 600000, 2200000 );
278 }
279 else if( STARTS_WITH_CI(papszProj[1], "(III) France Zone III") )
280 {
281 SetLCC1SP( 44.1, 2.337229167, 0.99987750, 600000, 3200000 );
282 }
283 else if( STARTS_WITH_CI(papszProj[1], "(IV) France Zone IV") )
284 {
285 SetLCC1SP( 42.165, 2.337229167, 0.99994471, 234.358, 4185861.369 );
286 }
287
288 /*
289 * Note: The following projections have not been implemented yet
290 *
291 */
292
293 /*
294 else if( STARTS_WITH_CI(papszProj[1], "(BNG) British National Grid") )
295 {
296 }
297 else if( STARTS_WITH_CI(papszProj[1], "(IG) Irish Grid") )
298 {
299 }
300
301 else if( STARTS_WITH_CI(papszProj[1], "(NZG) New Zealand Grid") )
302 {
303 }
304 else if( STARTS_WITH_CI(papszProj[1], "(NZTM2) New Zealand TM 2000") )
305 {
306 }
307 else if( STARTS_WITH_CI(papszProj[1], "(SG) Swedish Grid") )
308 {
309 }
310 else if( STARTS_WITH_CI(papszProj[1], "(SUI) Swiss Grid") )
311 {
312 }
313 else if( STARTS_WITH_CI(papszProj[1], "(A)Lambert Azimuthual Equal Area") )
314 {
315 }
316 else if( STARTS_WITH_CI(papszProj[1], "(EQC) Equidistant Conic") )
317 {
318 }
319 else if( STARTS_WITH_CI(papszProj[1], "Polyconic (American)") )
320 {
321 }
322 else if( STARTS_WITH_CI(papszProj[1], "Van Der Grinten") )
323 {
324 }
325 else if( STARTS_WITH_CI(papszProj[1], "Vertical Near-Sided Perspective") )
326 {
327 }
328 else if( STARTS_WITH_CI(papszProj[1], "(WIV) Wagner IV") )
329 {
330 }
331 else if( STARTS_WITH_CI(papszProj[1], "Bonne") )
332 {
333 }
334 else if( STARTS_WITH_CI(papszProj[1],
335 "(MT0) Montana State Plane Zone 2500") )
336 {
337 }
338 else if( STARTS_WITH_CI(papszProj[1], "ITA1) Italy Grid Zone 1") )
339 {
340 }
341 else if( STARTS_WITH_CI(papszProj[1], "ITA2) Italy Grid Zone 2") )
342 {
343 }
344 else if( STARTS_WITH_CI(papszProj[1],
345 "(VICMAP-TM) Victoria Aust.(pseudo AMG)") )
346 {
347 }
348 else if( STARTS_WITH_CI(papszProj[1], "VICGRID) Victoria Australia") )
349 {
350 }
351 else if( STARTS_WITH_CI(papszProj[1],
352 "(VG94) VICGRID94 Victoria Australia") )
353 {
354 }
355 else if( STARTS_WITH_CI(papszProj[1], "Gnomonic") )
356 {
357 }
358 */
359 else
360 {
361 CPLDebug( "OSR_Ozi", "Unsupported projection: \"%s\"", papszProj[1] );
362 SetLocalCS( CPLString().Printf(R"("Ozi" projection "%s")",
363 papszProj[1]) );
364 }
365
366 /* -------------------------------------------------------------------- */
367 /* Try to translate the datum/spheroid. */
368 /* -------------------------------------------------------------------- */
369 papszDatum = CSLTokenizeString2( pszDatum, ",",
370 CSLT_ALLOWEMPTYTOKENS
371 | CSLT_STRIPLEADSPACES
372 | CSLT_STRIPENDSPACES );
373 if( papszDatum == nullptr )
374 goto not_enough_data;
375
376 if( !IsLocal() )
377 {
378 /* -------------------------------------------------------------------- */
379 /* Verify that we can find the CSV file containing the datums */
380 /* -------------------------------------------------------------------- */
381 if( CSVScanFileByName( CSVFilename( "ozi_datum.csv" ),
382 "EPSG_DATUM_CODE",
383 "4326", CC_Integer ) == nullptr )
384 {
385 CPLError(CE_Failure, CPLE_OpenFailed,
386 "Unable to open OZI support file %s. "
387 "Try setting the GDAL_DATA environment variable to point "
388 "to the directory containing OZI csv files.",
389 CSVFilename( "ozi_datum.csv" ));
390 goto other_error;
391 }
392
393 /* -------------------------------------------------------------------- */
394 /* Search for matching datum */
395 /* -------------------------------------------------------------------- */
396 const char *pszOziDatum = CSVFilename( "ozi_datum.csv" );
397 CPLString osDName = CSVGetField( pszOziDatum, "NAME", papszDatum[0],
398 CC_ApproxString, "NAME" );
399 if( osDName.empty() )
400 {
401 CPLError( CE_Failure, CPLE_AppDefined,
402 "Failed to find datum %s in ozi_datum.csv.",
403 papszDatum[0] );
404 goto other_error;
405 }
406
407 const int nDatumCode =
408 atoi( CSVGetField( pszOziDatum, "NAME", papszDatum[0],
409 CC_ApproxString, "EPSG_DATUM_CODE" ) );
410
411 if( nDatumCode > 0 ) // There is a matching EPSG code
412 {
413 OGRSpatialReference oGCS;
414 oGCS.importFromEPSG( nDatumCode );
415 CopyGeogCSFrom( &oGCS );
416 }
417 else // We use the parameters from the CSV files
418 {
419 CPLString osEllipseCode =
420 CSVGetField( pszOziDatum, "NAME", papszDatum[0],
421 CC_ApproxString, "ELLIPSOID_CODE" );
422 const double dfDeltaX =
423 CPLAtof(CSVGetField( pszOziDatum, "NAME", papszDatum[0],
424 CC_ApproxString, "DELTAX" ) );
425 const double dfDeltaY =
426 CPLAtof(CSVGetField( pszOziDatum, "NAME", papszDatum[0],
427 CC_ApproxString, "DELTAY" ) );
428 const double dfDeltaZ =
429 CPLAtof(CSVGetField( pszOziDatum, "NAME", papszDatum[0],
430 CC_ApproxString, "DELTAZ" ) );
431
432 /* -------------------------------------------------------------------- */
433 /* Verify that we can find the CSV file containing the ellipsoids. */
434 /* -------------------------------------------------------------------- */
435 if( CSVScanFileByName( CSVFilename( "ozi_ellips.csv" ),
436 "ELLIPSOID_CODE",
437 "20", CC_Integer ) == nullptr )
438 {
439 CPLError(
440 CE_Failure, CPLE_OpenFailed,
441 "Unable to open OZI support file %s. "
442 "Try setting the GDAL_DATA environment variable to point "
443 "to the directory containing OZI csv files.",
444 CSVFilename( "ozi_ellips.csv" ) );
445 goto other_error;
446 }
447
448 /* -------------------------------------------------------------------- */
449 /* Lookup the ellipse code. */
450 /* -------------------------------------------------------------------- */
451 const char *pszOziEllipse = CSVFilename( "ozi_ellips.csv" );
452
453 CPLString osEName =
454 CSVGetField( pszOziEllipse, "ELLIPSOID_CODE", osEllipseCode,
455 CC_ApproxString, "NAME" );
456 if( osEName.empty() )
457 {
458 CPLError( CE_Failure, CPLE_AppDefined,
459 "Failed to find ellipsoid %s in ozi_ellips.csv.",
460 osEllipseCode.c_str() );
461 goto other_error;
462 }
463
464 const double dfA =
465 CPLAtof(CSVGetField( pszOziEllipse, "ELLIPSOID_CODE",
466 osEllipseCode, CC_ApproxString, "A" ));
467 const double dfInvF =
468 CPLAtof(CSVGetField( pszOziEllipse, "ELLIPSOID_CODE",
469 osEllipseCode, CC_ApproxString, "INVF" ));
470
471 /* -------------------------------------------------------------------- */
472 /* Create geographic coordinate system. */
473 /* -------------------------------------------------------------------- */
474 SetGeogCS( osDName, osDName, osEName, dfA, dfInvF );
475 SetTOWGS84( dfDeltaX, dfDeltaY, dfDeltaZ );
476 }
477 }
478
479 /* -------------------------------------------------------------------- */
480 /* Grid units translation */
481 /* -------------------------------------------------------------------- */
482 if( IsLocal() || IsProjected() )
483 SetLinearUnits( SRS_UL_METER, 1.0 );
484
485 CSLDestroy(papszProj);
486 CSLDestroy(papszProjParams);
487 CSLDestroy(papszDatum);
488
489 return OGRERR_NONE;
490
491 not_enough_data:
492
493 CSLDestroy(papszProj);
494 CSLDestroy(papszProjParams);
495 CSLDestroy(papszDatum);
496
497 return OGRERR_NOT_ENOUGH_DATA;
498
499 other_error:
500
501 CSLDestroy(papszProj);
502 CSLDestroy(papszProjParams);
503 CSLDestroy(papszDatum);
504
505 return OGRERR_FAILURE;
506 }
507