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