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