Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: OGRSpatialReference translation to/from PCI georeferencing
5 : * information.
6 : * Author: Andrey Kiselev, dron@ak4719.spb.edu
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2003, Andrey Kiselev <dron@ak4719.spb.edu>
10 : * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "cpl_port.h"
16 : #include "ogr_srs_api.h"
17 :
18 : #include <cctype>
19 : #include <cstdio>
20 : #include <cstdlib>
21 : #include <cstring>
22 :
23 : #include "cpl_conv.h"
24 : #include "cpl_csv.h"
25 : #include "cpl_error.h"
26 : #include "cpl_string.h"
27 : #include "cpl_vsi.h"
28 : #include "ogr_core.h"
29 : #include "ogr_p.h"
30 : #include "ogr_spatialref.h"
31 :
32 : // PCI uses a 16-character string for coordinate system and datum/ellipsoid.
33 : constexpr int knProjSize = 16;
34 :
35 : struct PCIDatums
36 : {
37 : const char *pszPCIDatum;
38 : int nEPSGCode;
39 : };
40 :
41 : static const PCIDatums asDatums[] = {
42 : {"D-01", 4267}, // NAD27 (USA, NADCON)
43 : {"D-03", 4267}, // NAD27 (Canada, NTv1)
44 : {"D-02", 4269}, // NAD83 (USA, NADCON)
45 : {"D-04", 4269}, // NAD83 (Canada, NTv1)
46 : {"D000", 4326}, // WGS 1984
47 : {"D001", 4322}, // WGS 1972
48 : {"D008", 4296}, // Sudan
49 : {"D013", 4601}, // Antigua Island Astro 1943
50 : {"D029", 4202}, // Australian Geodetic 1966
51 : {"D030", 4203}, // Australian Geodetic 1984
52 : {"D033", 4216}, // Bermuda 1957
53 : {"D034", 4165}, // Bissau
54 : {"D036", 4219}, // Bukit Rimpah
55 : {"D038", 4221}, // Campo Inchauspe
56 : {"D040", 4222}, // Cape
57 : {"D042", 4223}, // Carthage
58 : {"D044", 4224}, // Chua Astro
59 : {"D045", 4225}, // Corrego Alegre
60 : {"D046", 4155}, // Dabola (Guinea)
61 : {"D066", 4272}, // Geodetic Datum 1949 (New Zealand)
62 : {"D071", 4255}, // Herat North (Afghanistan)
63 : {"D077", 4239}, // Indian 1954 (Thailand, Vietnam)
64 : {"D078", 4240}, // Indian 1975 (Thailand)
65 : {"D083", 4244}, // Kandawala (Sri Lanka)
66 : {"D085", 4245}, // Kertau 1948 (West Malaysia & Singapore)
67 : {"D088", 4250}, // Leigon (Ghana)
68 : {"D089", 4251}, // Liberia 1964 (Liberia)
69 : {"D092", 4256}, // Mahe 1971 (Mahe Island)
70 : {"D093", 4262}, // Massawa (Ethiopia (Eritrea))
71 : {"D094", 4261}, // Merchich (Morocco)
72 : {"D098",
73 : 4604}, // Montserrat Island Astro 1958 (Montserrat (Leeward Islands))
74 : {"D110", 4267}, // NAD27 / Alaska
75 : {"D139", 4282}, // Pointe Noire 1948 (Congo)
76 : {"D140", 4615}, // Porto Santo 1936 (Porto Santo, Madeira Islands)
77 : {"D151", 4139}, // Puerto Rico (Puerto Rico, Virgin Islands)
78 : {"D153", 4287}, // Qornoq (Greenland (South))
79 : {"D158", 4292}, // Sapper Hill 1943 (East Falkland Island)
80 : {"D159", 4293}, // Schwarzeck (Namibia)
81 : {"D160", 4616}, // Selvagem Grande 1938 (Salvage Islands)
82 : {"D176", 4297}, // Tananarive Observatory 1925 (Madagascar)
83 : {"D177", 4298}, // Timbalai 1948 (Brunei, East Malaysia (Sabah, Sarawak))
84 : {"D187", 4309}, // Yacare (Uruguay)
85 : {"D188", 4311}, // Zanderij (Suriname)
86 : {"D401", 4124}, // RT90 (Sweden)
87 : {"D501", 4312}, // MGI (Hermannskogel, Austria)
88 : {nullptr, 0}};
89 :
90 : static const PCIDatums asEllips[] = {
91 : {"E000", 7008}, // Clarke, 1866 (NAD1927)
92 : {"E001", 7034}, // Clarke, 1880
93 : {"E002", 7004}, // Bessel, 1841
94 : {"E004", 7022}, // International, 1924 (Hayford, 1909)
95 : {"E005", 7043}, // WGS, 1972
96 : {"E006", 7042}, // Everest, 1830
97 : {"E008", 7019}, // GRS, 1980 (NAD1983)
98 : {"E009", 7001}, // Airy, 1830
99 : {"E010", 7018}, // Modified Everest
100 : {"E011", 7002}, // Modified Airy
101 : {"E012", 7030}, // WGS, 1984 (GPS)
102 : {"E014", 7003}, // Australian National, 1965
103 : {"E015", 7024}, // Krassovsky, 1940
104 : {"E016", 7053}, // Hough
105 : {"E019", 7052}, // normal sphere
106 : {"E333", 7046}, // Bessel 1841 (Japan By Law)
107 : {"E900", 7006}, // Bessel, 1841 (Namibia)
108 : {"E901", 7044}, // Everest, 1956
109 : {"E902", 7056}, // Everest, 1969
110 : {"E903", 7016}, // Everest (Sabah & Sarawak)
111 : {"E904", 7020}, // Helmert, 1906
112 : {"E907", 7036}, // South American, 1969
113 : {"E910", 7041}, // ATS77
114 : {nullptr, 0}};
115 :
116 : /************************************************************************/
117 : /* OSRImportFromPCI() */
118 : /************************************************************************/
119 :
120 : /**
121 : * \brief Import coordinate system from PCI projection definition.
122 : *
123 : * This function is the same as OGRSpatialReference::importFromPCI().
124 : */
125 :
126 8 : OGRErr OSRImportFromPCI(OGRSpatialReferenceH hSRS, const char *pszProj,
127 : const char *pszUnits, double *padfPrjParams)
128 :
129 : {
130 8 : VALIDATE_POINTER1(hSRS, "OSRImportFromPCI", OGRERR_FAILURE);
131 :
132 8 : return OGRSpatialReference::FromHandle(hSRS)->importFromPCI(
133 8 : pszProj, pszUnits, padfPrjParams);
134 : }
135 :
136 : /************************************************************************/
137 : /* importFromPCI() */
138 : /************************************************************************/
139 :
140 : /**
141 : * \brief Import coordinate system from PCI projection definition.
142 : *
143 : * PCI software uses 16-character string to specify coordinate system
144 : * and datum/ellipsoid. You should supply at least this string to the
145 : * importFromPCI() function.
146 : *
147 : * This function is the equivalent of the C function OSRImportFromPCI().
148 : *
149 : * @param pszProj NULL terminated string containing the definition. Looks
150 : * like "pppppppppppp Ennn" or "pppppppppppp Dnnn", where "pppppppppppp" is
151 : * a projection code, "Ennn" is an ellipsoid code, "Dnnn" a datum code.
152 : *
153 : * @param pszUnits Grid units code ("DEGREE" or "METRE"). If NULL "METRE" will
154 : * be used.
155 : *
156 : * @param padfPrjParams Array of 17 coordinate system parameters:
157 : *
158 : * [0] Spheroid semi major axis
159 : * [1] Spheroid semi minor axis
160 : * [2] Reference Longitude
161 : * [3] Reference Latitude
162 : * [4] First Standard Parallel
163 : * [5] Second Standard Parallel
164 : * [6] False Easting
165 : * [7] False Northing
166 : * [8] Scale Factor
167 : * [9] Height above sphere surface
168 : * [10] Longitude of 1st point on center line
169 : * [11] Latitude of 1st point on center line
170 : * [12] Longitude of 2nd point on center line
171 : * [13] Latitude of 2nd point on center line
172 : * [14] Azimuth east of north for center line
173 : * [15] Landsat satellite number
174 : * [16] Landsat path number
175 : *
176 : * Particular projection uses different parameters, unused ones may be set to
177 : * zero. If NULL is supplied instead of an array pointer, default values will be
178 : * used (i.e., zeroes).
179 : *
180 : * @return OGRERR_NONE on success or an error code in case of failure.
181 : */
182 :
183 1205 : OGRErr OGRSpatialReference::importFromPCI(const char *pszProj,
184 : const char *pszUnits,
185 : const double *padfPrjParams)
186 :
187 : {
188 1205 : Clear();
189 :
190 2410 : if (pszProj == nullptr ||
191 1205 : CPLStrnlen(pszProj, knProjSize) < static_cast<size_t>(knProjSize))
192 : {
193 3 : return OGRERR_CORRUPT_DATA;
194 : }
195 :
196 1202 : CPLDebug("OSR_PCI", "Trying to import projection \"%s\"", pszProj);
197 :
198 : /* -------------------------------------------------------------------- */
199 : /* Use safe defaults if projection parameters are not supplied. */
200 : /* -------------------------------------------------------------------- */
201 : static const double adfZeroedPrjParams[17] = {0};
202 1202 : if (padfPrjParams == nullptr)
203 : {
204 0 : padfPrjParams = adfZeroedPrjParams;
205 : }
206 :
207 : /* -------------------------------------------------------------------- */
208 : /* Extract and "normalize" the earthmodel to look like E001, */
209 : /* D-02 or D109. */
210 : /* -------------------------------------------------------------------- */
211 1202 : char szEarthModel[5] = {};
212 :
213 1202 : strcpy(szEarthModel, "");
214 1202 : const char *pszEM = pszProj + strlen(pszProj) - 1;
215 14108 : while (pszEM != pszProj)
216 : {
217 14108 : if (*pszEM == 'e' || *pszEM == 'E' || *pszEM == 'd' || *pszEM == 'D')
218 : {
219 1202 : int nCode = atoi(pszEM + 1);
220 :
221 1202 : if (nCode >= -99 && nCode <= 999)
222 1202 : snprintf(szEarthModel, sizeof(szEarthModel), "%c%03d",
223 1202 : CPLToupper(static_cast<unsigned char>(*pszEM)), nCode);
224 :
225 1202 : break;
226 : }
227 :
228 12906 : pszEM--;
229 : }
230 :
231 1196 : const bool bIsNAD27 = EQUAL(pszEM, "E000") || EQUAL(pszEM, "D-01") ||
232 1196 : EQUAL(pszEM, "D-03") || EQUAL(pszEM, "D-07") ||
233 1196 : EQUAL(pszEM, "D-09") || EQUAL(pszEM, "D-11") ||
234 2398 : EQUAL(pszEM, "D-13") || EQUAL(pszEM, "D-17");
235 :
236 : /* -------------------------------------------------------------------- */
237 : /* Operate on the basis of the projection name. */
238 : /* -------------------------------------------------------------------- */
239 1202 : if (STARTS_WITH_CI(pszProj, "LONG/LAT"))
240 : {
241 : // TODO(schwehr): A NOP is okay?
242 : }
243 1193 : else if (STARTS_WITH_CI(pszProj, "METER") ||
244 1186 : STARTS_WITH_CI(pszProj, "METRE"))
245 : {
246 1171 : SetLocalCS("METER");
247 1171 : SetLinearUnits("METER", 1.0);
248 : }
249 22 : else if (STARTS_WITH_CI(pszProj, "FEET") || STARTS_WITH_CI(pszProj, "FOOT"))
250 : {
251 0 : SetLocalCS("FEET");
252 0 : SetLinearUnits("FEET", CPLAtof(SRS_UL_FOOT_CONV));
253 : }
254 22 : else if (STARTS_WITH_CI(pszProj, "ACEA"))
255 : {
256 0 : SetACEA(padfPrjParams[4], padfPrjParams[5], padfPrjParams[3],
257 0 : padfPrjParams[2], padfPrjParams[6], padfPrjParams[7]);
258 : }
259 22 : else if (STARTS_WITH_CI(pszProj, "AE"))
260 : {
261 0 : SetAE(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
262 0 : padfPrjParams[7]);
263 : }
264 22 : else if (STARTS_WITH_CI(pszProj, "CASS "))
265 : {
266 0 : SetCS(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
267 0 : padfPrjParams[7]);
268 : }
269 22 : else if (STARTS_WITH_CI(pszProj, "EC"))
270 : {
271 2 : SetEC(padfPrjParams[4], padfPrjParams[5], padfPrjParams[3],
272 2 : padfPrjParams[2], padfPrjParams[6], padfPrjParams[7]);
273 : }
274 20 : else if (STARTS_WITH_CI(pszProj, "ER"))
275 : {
276 : // PCI and GCTP don't support natural origin lat.
277 0 : SetEquirectangular2(0.0, padfPrjParams[2], padfPrjParams[3],
278 0 : padfPrjParams[6], padfPrjParams[7]);
279 : }
280 20 : else if (STARTS_WITH_CI(pszProj, "GNO"))
281 : {
282 0 : SetGnomonic(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
283 0 : padfPrjParams[7]);
284 : }
285 : // FIXME: GVNP --- General Vertical Near- Side Perspective skipped.
286 : // FIXME: GOOD -- Our Goode's is not the interrupted version from PCI.
287 20 : else if (STARTS_WITH_CI(pszProj, "LAEA"))
288 : {
289 0 : SetLAEA(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
290 0 : padfPrjParams[7]);
291 : }
292 20 : else if (STARTS_WITH_CI(pszProj, "LCC "))
293 : {
294 4 : SetLCC(padfPrjParams[4], padfPrjParams[5], padfPrjParams[3],
295 4 : padfPrjParams[2], padfPrjParams[6], padfPrjParams[7]);
296 : }
297 16 : else if (STARTS_WITH_CI(pszProj, "LCC_1SP "))
298 : {
299 0 : SetLCC1SP(padfPrjParams[3], padfPrjParams[2], padfPrjParams[8],
300 0 : padfPrjParams[6], padfPrjParams[7]);
301 : }
302 16 : else if (STARTS_WITH_CI(pszProj, "MC"))
303 : {
304 0 : SetMC(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
305 0 : padfPrjParams[7]);
306 : }
307 16 : else if (STARTS_WITH_CI(pszProj, "MER"))
308 : {
309 1 : if (EQUAL(pszEM, "D894") && padfPrjParams[3] == 0 &&
310 1 : padfPrjParams[2] == 0 && padfPrjParams[8] == 1.0 &&
311 1 : padfPrjParams[6] == 0 && padfPrjParams[7] == 0)
312 : {
313 : // Special case for Web Mercator
314 1 : return importFromEPSG(3857);
315 : }
316 0 : SetMercator(padfPrjParams[3], padfPrjParams[2],
317 0 : (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
318 0 : padfPrjParams[6], padfPrjParams[7]);
319 : }
320 15 : else if (STARTS_WITH_CI(pszProj, "OG"))
321 : {
322 0 : SetOrthographic(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
323 0 : padfPrjParams[7]);
324 : }
325 15 : else if (STARTS_WITH_CI(pszProj, "OM "))
326 : {
327 0 : if (padfPrjParams[10] == 0.0 && padfPrjParams[11] == 0.0 &&
328 0 : padfPrjParams[12] == 0.0 && padfPrjParams[13] == 0.0)
329 : {
330 0 : SetHOM(padfPrjParams[3], padfPrjParams[2], padfPrjParams[14],
331 0 : padfPrjParams[14], // Use azimuth for grid angle.
332 0 : padfPrjParams[8], padfPrjParams[6], padfPrjParams[7]);
333 : }
334 : else
335 : {
336 0 : SetHOM2PNO(padfPrjParams[3], padfPrjParams[11], padfPrjParams[10],
337 0 : padfPrjParams[13], padfPrjParams[12], padfPrjParams[8],
338 0 : padfPrjParams[6], padfPrjParams[7]);
339 : }
340 : }
341 15 : else if (STARTS_WITH_CI(pszProj, "PC"))
342 : {
343 0 : SetPolyconic(padfPrjParams[3], padfPrjParams[2], padfPrjParams[6],
344 0 : padfPrjParams[7]);
345 : }
346 15 : else if (STARTS_WITH_CI(pszProj, "PS"))
347 : {
348 0 : SetPS(padfPrjParams[3], padfPrjParams[2],
349 0 : (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
350 0 : padfPrjParams[6], padfPrjParams[7]);
351 : }
352 15 : else if (STARTS_WITH_CI(pszProj, "ROB"))
353 : {
354 0 : SetRobinson(padfPrjParams[2], padfPrjParams[6], padfPrjParams[7]);
355 : }
356 :
357 15 : else if (STARTS_WITH_CI(pszProj, "SGDO"))
358 : {
359 0 : SetOS(padfPrjParams[3], padfPrjParams[2],
360 0 : (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
361 0 : padfPrjParams[6], padfPrjParams[7]);
362 : }
363 15 : else if (STARTS_WITH_CI(pszProj, "SG"))
364 : {
365 0 : SetStereographic(padfPrjParams[3], padfPrjParams[2],
366 0 : (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
367 0 : padfPrjParams[6], padfPrjParams[7]);
368 : }
369 15 : else if (STARTS_WITH_CI(pszProj, "SIN"))
370 : {
371 0 : SetSinusoidal(padfPrjParams[2], padfPrjParams[6], padfPrjParams[7]);
372 : }
373 : // FIXME: SOM --- Space Oblique Mercator skipped.
374 15 : else if (STARTS_WITH_CI(pszProj, "SPCS"))
375 : {
376 : const int iZone =
377 0 : static_cast<int>(CPLScanLong(const_cast<char *>(pszProj) + 5, 4));
378 :
379 0 : SetStatePlane(iZone, !bIsNAD27);
380 0 : SetLinearUnitsAndUpdateParameters(SRS_UL_METER, 1.0);
381 : }
382 15 : else if (STARTS_WITH_CI(pszProj, "SPIF"))
383 : {
384 : const int iZone =
385 0 : static_cast<int>(CPLScanLong(const_cast<char *>(pszProj) + 5, 4));
386 :
387 0 : SetStatePlane(iZone, !bIsNAD27);
388 0 : SetLinearUnitsAndUpdateParameters(SRS_UL_FOOT,
389 : CPLAtof(SRS_UL_FOOT_CONV));
390 : }
391 15 : else if (STARTS_WITH_CI(pszProj, "SPAF"))
392 : {
393 : const int iZone =
394 0 : static_cast<int>(CPLScanLong(const_cast<char *>(pszProj) + 5, 4));
395 :
396 0 : SetStatePlane(iZone, !bIsNAD27);
397 0 : SetLinearUnitsAndUpdateParameters(SRS_UL_US_FOOT,
398 : CPLAtof(SRS_UL_US_FOOT_CONV));
399 : }
400 15 : else if (STARTS_WITH_CI(pszProj, "TM"))
401 : {
402 0 : SetTM(padfPrjParams[3], padfPrjParams[2],
403 0 : (padfPrjParams[8] != 0.0) ? padfPrjParams[8] : 1.0,
404 0 : padfPrjParams[6], padfPrjParams[7]);
405 : }
406 :
407 15 : else if (STARTS_WITH_CI(pszProj, "UTM"))
408 : {
409 15 : bool bNorth = true;
410 :
411 : int iZone =
412 15 : static_cast<int>(CPLScanLong(const_cast<char *>(pszProj) + 4, 5));
413 15 : if (iZone < 0)
414 : {
415 0 : iZone = -iZone;
416 0 : bNorth = false;
417 : }
418 :
419 : // Check for a zone letter. PCI uses, accidentally, MGRS
420 : // type row lettering in its UTM projection.
421 15 : char byZoneID = 0;
422 :
423 15 : if (strlen(pszProj) > 10 && pszProj[10] != ' ')
424 4 : byZoneID = pszProj[10];
425 :
426 : // Determine if the MGRS zone falls above or below the equator.
427 15 : if (byZoneID != 0)
428 : {
429 4 : CPLDebug("OSR_PCI", "Found MGRS zone in UTM projection string: %c",
430 : byZoneID);
431 :
432 4 : if (byZoneID >= 'N' && byZoneID <= 'X')
433 : {
434 3 : bNorth = true;
435 : }
436 1 : else if (byZoneID >= 'C' && byZoneID <= 'M')
437 : {
438 1 : bNorth = false;
439 : }
440 : else
441 : {
442 : // Yikes. Most likely we got something that was not really
443 : // an MGRS zone code so we ignore it.
444 : }
445 : }
446 :
447 15 : SetUTM(iZone, bNorth);
448 : }
449 0 : else if (STARTS_WITH_CI(pszProj, "VDG"))
450 : {
451 0 : SetVDG(padfPrjParams[2], padfPrjParams[6], padfPrjParams[7]);
452 : }
453 : else
454 : {
455 0 : CPLDebug("OSR_PCI", "Unsupported projection: %s", pszProj);
456 0 : SetLocalCS(pszProj);
457 : }
458 :
459 : /* ==================================================================== */
460 : /* Translate the datum/spheroid. */
461 : /* ==================================================================== */
462 :
463 : /* -------------------------------------------------------------------- */
464 : /* We have an earthmodel string, look it up in the datum list. */
465 : /* -------------------------------------------------------------------- */
466 2402 : if (strlen(szEarthModel) > 0 &&
467 1201 : (GetRoot() == nullptr || IsProjected() || IsGeographic()))
468 : {
469 30 : const PCIDatums *pasDatum = asDatums;
470 :
471 : // Search for matching datum.
472 729 : while (pasDatum->pszPCIDatum)
473 : {
474 715 : if (EQUALN(szEarthModel, pasDatum->pszPCIDatum, 4))
475 : {
476 32 : OGRSpatialReference oGCS;
477 16 : oGCS.importFromEPSG(pasDatum->nEPSGCode);
478 16 : CopyGeogCSFrom(&oGCS);
479 16 : break;
480 : }
481 699 : pasDatum++;
482 : }
483 :
484 : /* --------------------------------------------------------------------
485 : */
486 : /* If we did not find a datum definition in our in-code EPSG */
487 : /* lookup table, then try fetching from the pci_datum.txt */
488 : /* file. */
489 : /* --------------------------------------------------------------------
490 : */
491 30 : char **papszDatumDefn = nullptr;
492 :
493 30 : if (!pasDatum->pszPCIDatum && szEarthModel[0] == 'D')
494 : {
495 1 : const char *pszDatumCSV = CSVFilename("pci_datum.txt");
496 1 : VSILFILE *fp = pszDatumCSV ? VSIFOpenL(pszDatumCSV, "r") : nullptr;
497 :
498 1 : if (fp != nullptr)
499 : {
500 1 : char **papszLineItems = nullptr;
501 :
502 315 : while ((papszLineItems = CSVReadParseLineL(fp)) != nullptr)
503 : {
504 607 : if (CSLCount(papszLineItems) > 3 &&
505 292 : EQUALN(papszLineItems[0], szEarthModel,
506 : sizeof(szEarthModel) - 1))
507 : {
508 1 : papszDatumDefn = papszLineItems;
509 1 : strncpy(szEarthModel, papszLineItems[2],
510 : sizeof(szEarthModel) - 1);
511 1 : break;
512 : }
513 314 : CSLDestroy(papszLineItems);
514 : }
515 :
516 1 : VSIFCloseL(fp);
517 : }
518 : }
519 :
520 : /* --------------------------------------------------------------------
521 : */
522 : /* If not, look in the ellipsoid/EPSG matching list. */
523 : /* --------------------------------------------------------------------
524 : */
525 30 : if (!pasDatum->pszPCIDatum) // No matching; search for ellipsoids.
526 : {
527 14 : char *pszName = nullptr;
528 14 : double dfSemiMajor = 0.0;
529 14 : double dfInvFlattening = 0.0;
530 14 : int nEPSGCode = 0;
531 :
532 14 : pasDatum = asEllips;
533 :
534 63 : while (pasDatum->pszPCIDatum)
535 : {
536 62 : if (EQUALN(szEarthModel, pasDatum->pszPCIDatum, 4))
537 : {
538 13 : nEPSGCode = pasDatum->nEPSGCode;
539 13 : CPL_IGNORE_RET_VAL(
540 13 : OSRGetEllipsoidInfo(pasDatum->nEPSGCode, &pszName,
541 : &dfSemiMajor, &dfInvFlattening));
542 13 : break;
543 : }
544 49 : pasDatum++;
545 : }
546 :
547 : /* --------------------------------------------------------------------
548 : */
549 : /* If we don't find it in that list, do a lookup in the */
550 : /* pci_ellips.txt file. */
551 : /* --------------------------------------------------------------------
552 : */
553 14 : if (!pasDatum->pszPCIDatum && szEarthModel[0] == 'E')
554 : {
555 1 : const char *pszCSV = CSVFilename("pci_ellips.txt");
556 1 : VSILFILE *fp = pszCSV ? VSIFOpenL(pszCSV, "r") : nullptr;
557 :
558 1 : if (fp != nullptr)
559 : {
560 1 : char **papszLineItems = nullptr;
561 :
562 86 : while ((papszLineItems = CSVReadParseLineL(fp)) != nullptr)
563 : {
564 138 : if (CSLCount(papszLineItems) > 3 &&
565 52 : EQUALN(papszLineItems[0], szEarthModel, 4))
566 : {
567 1 : dfSemiMajor = CPLAtof(papszLineItems[2]);
568 : const double dfSemiMinor =
569 1 : CPLAtof(papszLineItems[3]);
570 1 : dfInvFlattening =
571 1 : OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor);
572 1 : break;
573 : }
574 85 : CSLDestroy(papszLineItems);
575 : }
576 1 : CSLDestroy(papszLineItems);
577 :
578 1 : VSIFCloseL(fp);
579 : }
580 : }
581 :
582 : /* --------------------------------------------------------------------
583 : */
584 : /* Custom spheroid? */
585 : /* --------------------------------------------------------------------
586 : */
587 14 : if (dfSemiMajor == 0.0 && STARTS_WITH_CI(szEarthModel, "E999") &&
588 0 : padfPrjParams[0] != 0.0)
589 : {
590 0 : dfSemiMajor = padfPrjParams[0];
591 0 : double dfSemiMinor = padfPrjParams[1];
592 0 : dfInvFlattening =
593 0 : OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor);
594 : }
595 :
596 : /* --------------------------------------------------------------------
597 : */
598 : /* If nothing else, fall back to WGS84 parameters. */
599 : /* --------------------------------------------------------------------
600 : */
601 14 : if (dfSemiMajor == 0.0)
602 : {
603 0 : dfSemiMajor = SRS_WGS84_SEMIMAJOR;
604 0 : dfInvFlattening = SRS_WGS84_INVFLATTENING;
605 : }
606 :
607 : /* --------------------------------------------------------------------
608 : */
609 : /* Now try to put this all together into a GEOGCS definition.
610 : */
611 : /* --------------------------------------------------------------------
612 : */
613 :
614 28 : CPLString osEllipseName;
615 14 : if (pszName)
616 13 : osEllipseName = pszName;
617 : else
618 1 : osEllipseName.Printf("Unknown - PCI %s", szEarthModel);
619 14 : CPLFree(pszName);
620 :
621 28 : CPLString osDatumName;
622 14 : if (papszDatumDefn)
623 1 : osDatumName = papszDatumDefn[1];
624 : else
625 13 : osDatumName.Printf("Unknown - PCI %s", szEarthModel);
626 :
627 28 : const CPLString osGCSName = osDatumName;
628 :
629 14 : SetGeogCS(osGCSName, osDatumName, osEllipseName, dfSemiMajor,
630 : dfInvFlattening);
631 :
632 : // Do we have an ellipsoid EPSG code?
633 14 : if (nEPSGCode != 0)
634 13 : SetAuthority("SPHEROID", "EPSG", nEPSGCode);
635 :
636 : // Do we have 7 datum shift parameters?
637 15 : if (papszDatumDefn != nullptr && CSLCount(papszDatumDefn) >= 15 &&
638 1 : CPLAtof(papszDatumDefn[14]) != 0.0)
639 : {
640 1 : double dfScale = CPLAtof(papszDatumDefn[14]);
641 :
642 : // we want scale in parts per million off 1.0
643 : // but pci uses a mix of forms.
644 1 : if (dfScale >= 0.999 && dfScale <= 1.001)
645 1 : dfScale = (dfScale - 1.0) * 1000000.0;
646 :
647 1 : SetTOWGS84(
648 1 : CPLAtof(papszDatumDefn[3]), CPLAtof(papszDatumDefn[4]),
649 1 : CPLAtof(papszDatumDefn[5]), CPLAtof(papszDatumDefn[11]),
650 1 : CPLAtof(papszDatumDefn[12]), CPLAtof(papszDatumDefn[13]),
651 : dfScale);
652 : }
653 :
654 : // Do we have 7 datum shift parameters?
655 13 : else if (papszDatumDefn != nullptr &&
656 13 : CSLCount(papszDatumDefn) == 11 &&
657 0 : (CPLAtof(papszDatumDefn[3]) != 0.0 ||
658 0 : CPLAtof(papszDatumDefn[4]) != 0.0 ||
659 0 : CPLAtof(papszDatumDefn[5]) != 0.0))
660 : {
661 0 : SetTOWGS84(CPLAtof(papszDatumDefn[3]),
662 0 : CPLAtof(papszDatumDefn[4]),
663 0 : CPLAtof(papszDatumDefn[5]));
664 : }
665 : }
666 :
667 30 : CSLDestroy(papszDatumDefn);
668 : }
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Grid units translation */
672 : /* -------------------------------------------------------------------- */
673 1201 : if ((IsLocal() || IsProjected()) && pszUnits)
674 : {
675 5 : if (EQUAL(pszUnits, "METRE"))
676 5 : SetLinearUnits(SRS_UL_METER, 1.0);
677 0 : else if (EQUAL(pszUnits, "DEGREE"))
678 0 : SetAngularUnits(SRS_UA_DEGREE, CPLAtof(SRS_UA_DEGREE_CONV));
679 : else
680 0 : SetLinearUnits(SRS_UL_METER, 1.0);
681 : }
682 :
683 1201 : return OGRERR_NONE;
684 : }
685 :
686 : /************************************************************************/
687 : /* OSRExportToPCI() */
688 : /************************************************************************/
689 : /**
690 : * \brief Export coordinate system in PCI projection definition.
691 : *
692 : * This function is the same as OGRSpatialReference::exportToPCI().
693 : */
694 7 : OGRErr OSRExportToPCI(OGRSpatialReferenceH hSRS, char **ppszProj,
695 : char **ppszUnits, double **ppadfPrjParams)
696 :
697 : {
698 7 : VALIDATE_POINTER1(hSRS, "OSRExportToPCI", OGRERR_FAILURE);
699 :
700 7 : *ppszProj = nullptr;
701 7 : *ppszUnits = nullptr;
702 7 : *ppadfPrjParams = nullptr;
703 :
704 7 : return OGRSpatialReference::FromHandle(hSRS)->exportToPCI(
705 7 : ppszProj, ppszUnits, ppadfPrjParams);
706 : }
707 :
708 : /************************************************************************/
709 : /* exportToPCI() */
710 : /************************************************************************/
711 :
712 : /**
713 : * \brief Export coordinate system in PCI projection definition.
714 : *
715 : * Converts the loaded coordinate reference system into PCI projection
716 : * definition to the extent possible. The strings returned in ppszProj,
717 : * ppszUnits and ppadfPrjParams array should be deallocated by the caller
718 : * with CPLFree() when no longer needed.
719 : *
720 : * LOCAL_CS coordinate systems are not translatable. An empty string
721 : * will be returned along with OGRERR_NONE.
722 : *
723 : * This method is the equivalent of the C function OSRExportToPCI().
724 : *
725 : * @param ppszProj pointer to which dynamically allocated PCI projection
726 : * definition will be assigned.
727 : *
728 : * @param ppszUnits pointer to which dynamically allocated units definition
729 : * will be assigned.
730 : *
731 : * @param ppadfPrjParams pointer to which dynamically allocated array of
732 : * 17 projection parameters will be assigned. See importFromPCI() for the list
733 : * of parameters.
734 : *
735 : * @return OGRERR_NONE on success or an error code on failure.
736 : */
737 :
738 69 : OGRErr OGRSpatialReference::exportToPCI(char **ppszProj, char **ppszUnits,
739 : double **ppadfPrjParams) const
740 :
741 : {
742 69 : const char *pszProjection = GetAttrValue("PROJECTION");
743 :
744 : /* -------------------------------------------------------------------- */
745 : /* Fill all projection parameters with zero. */
746 : /* -------------------------------------------------------------------- */
747 69 : *ppadfPrjParams = static_cast<double *>(CPLMalloc(17 * sizeof(double)));
748 1242 : for (int i = 0; i < 17; i++)
749 1173 : (*ppadfPrjParams)[i] = 0.0;
750 :
751 : /* -------------------------------------------------------------------- */
752 : /* Get the prime meridian info. */
753 : /* -------------------------------------------------------------------- */
754 : #if 0
755 : const OGR_SRSNode *poPRIMEM = GetAttrNode( "PRIMEM" );
756 : double dfFromGreenwich = 0.0;
757 :
758 : if( poPRIMEM != NULL && poPRIMEM->GetChildCount() >= 2
759 : && CPLAtof(poPRIMEM->GetChild(1)->GetValue()) != 0.0 )
760 : {
761 : dfFromGreenwich = CPLAtof(poPRIMEM->GetChild(1)->GetValue());
762 : }
763 : #endif
764 :
765 : /* ==================================================================== */
766 : /* Handle the projection definition. */
767 : /* ==================================================================== */
768 69 : char szProj[knProjSize + 1] = {};
769 :
770 69 : if (IsLocal())
771 : {
772 2 : if (GetLinearUnits() > 0.30479999 && GetLinearUnits() < 0.3048010)
773 0 : CPLPrintStringFill(szProj, "FEET", knProjSize);
774 : else
775 2 : CPLPrintStringFill(szProj, "METER", knProjSize);
776 : }
777 67 : else if (pszProjection == nullptr)
778 : {
779 54 : CPLPrintStringFill(szProj, "LONG/LAT", knProjSize);
780 : }
781 13 : else if (EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
782 : {
783 0 : CPLPrintStringFill(szProj, "ACEA", knProjSize);
784 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
785 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
786 0 : (*ppadfPrjParams)[4] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
787 0 : (*ppadfPrjParams)[5] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0);
788 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
789 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
790 : }
791 13 : else if (EQUAL(pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT))
792 : {
793 0 : CPLPrintStringFill(szProj, "AE", knProjSize);
794 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
795 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
796 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
797 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
798 : }
799 13 : else if (EQUAL(pszProjection, SRS_PT_CASSINI_SOLDNER))
800 : {
801 0 : CPLPrintStringFill(szProj, "CASS", knProjSize);
802 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
803 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
804 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
805 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
806 : }
807 13 : else if (EQUAL(pszProjection, SRS_PT_EQUIDISTANT_CONIC))
808 : {
809 1 : CPLPrintStringFill(szProj, "EC", knProjSize);
810 1 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
811 1 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
812 1 : (*ppadfPrjParams)[4] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
813 1 : (*ppadfPrjParams)[5] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0);
814 1 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
815 1 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
816 : }
817 12 : else if (EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR))
818 : {
819 0 : CPLPrintStringFill(szProj, "ER", knProjSize);
820 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
821 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
822 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
823 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
824 : }
825 12 : else if (EQUAL(pszProjection, SRS_PT_GNOMONIC))
826 : {
827 0 : CPLPrintStringFill(szProj, "GNO", knProjSize);
828 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
829 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
830 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
831 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
832 : }
833 12 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
834 : {
835 0 : CPLPrintStringFill(szProj, "LAEA", knProjSize);
836 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
837 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
838 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
839 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
840 : }
841 12 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
842 : {
843 2 : CPLPrintStringFill(szProj, "LCC", knProjSize);
844 2 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
845 2 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
846 2 : (*ppadfPrjParams)[4] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0);
847 2 : (*ppadfPrjParams)[5] = GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0);
848 2 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
849 2 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
850 : }
851 10 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
852 : {
853 0 : CPLPrintStringFill(szProj, "LCC_1SP", knProjSize);
854 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
855 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
856 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
857 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
858 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
859 : }
860 10 : else if (EQUAL(pszProjection, SRS_PT_MILLER_CYLINDRICAL))
861 : {
862 0 : CPLPrintStringFill(szProj, "MC", knProjSize);
863 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
864 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
865 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
866 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
867 : }
868 10 : else if (EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
869 : {
870 1 : CPLPrintStringFill(szProj, "MER", knProjSize);
871 1 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
872 1 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
873 1 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
874 1 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
875 1 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
876 : }
877 9 : else if (EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC))
878 : {
879 0 : CPLPrintStringFill(szProj, "OG", knProjSize);
880 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
881 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
882 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
883 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
884 : }
885 9 : else if (EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
886 : {
887 0 : CPLPrintStringFill(szProj, "OM", knProjSize);
888 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
889 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
890 0 : (*ppadfPrjParams)[14] = GetNormProjParm(SRS_PP_AZIMUTH, 0.0);
891 : // Note: Ignoring rectified_grid_angle which has no PCI analog.
892 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0);
893 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
894 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
895 : }
896 9 : else if (EQUAL(pszProjection,
897 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN))
898 : {
899 0 : CPLPrintStringFill(szProj, "OM", knProjSize);
900 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_CENTER, 0.0);
901 0 : (*ppadfPrjParams)[11] =
902 0 : GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_1, 0.0);
903 0 : (*ppadfPrjParams)[10] =
904 0 : GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_1, 0.0);
905 0 : (*ppadfPrjParams)[13] =
906 0 : GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_2, 0.0);
907 0 : (*ppadfPrjParams)[12] =
908 0 : GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_2, 0.0);
909 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0);
910 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
911 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
912 : }
913 9 : else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
914 : {
915 0 : CPLPrintStringFill(szProj, "PC", knProjSize);
916 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
917 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
918 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
919 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
920 : }
921 9 : else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
922 : {
923 0 : CPLPrintStringFill(szProj, "PS", knProjSize);
924 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
925 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
926 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
927 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
928 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
929 : }
930 9 : else if (EQUAL(pszProjection, SRS_PT_ROBINSON))
931 : {
932 0 : CPLPrintStringFill(szProj, "ROB", knProjSize);
933 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
934 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
935 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
936 : }
937 9 : else if (EQUAL(pszProjection, SRS_PT_OBLIQUE_STEREOGRAPHIC))
938 : {
939 0 : CPLPrintStringFill(szProj, "SGDO", knProjSize);
940 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
941 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
942 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
943 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
944 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
945 : }
946 9 : else if (EQUAL(pszProjection, SRS_PT_STEREOGRAPHIC))
947 : {
948 0 : CPLPrintStringFill(szProj, "SG", knProjSize);
949 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
950 0 : (*ppadfPrjParams)[3] = GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
951 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
952 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
953 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
954 : }
955 9 : else if (EQUAL(pszProjection, SRS_PT_SINUSOIDAL))
956 : {
957 0 : CPLPrintStringFill(szProj, "SIN", knProjSize);
958 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER, 0.0);
959 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
960 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
961 : }
962 9 : else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
963 : {
964 9 : int bNorth = FALSE;
965 9 : int nZone = GetUTMZone(&bNorth);
966 :
967 9 : if (nZone != 0)
968 : {
969 9 : CPLPrintStringFill(szProj, "UTM", knProjSize);
970 9 : if (bNorth)
971 9 : CPLPrintInt32(szProj + 5, nZone, 4);
972 : else
973 0 : CPLPrintInt32(szProj + 5, -nZone, 4);
974 : }
975 : else
976 : {
977 0 : CPLPrintStringFill(szProj, "TM", knProjSize);
978 0 : (*ppadfPrjParams)[2] =
979 0 : GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
980 0 : (*ppadfPrjParams)[3] =
981 0 : GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
982 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
983 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
984 0 : (*ppadfPrjParams)[8] = GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0);
985 : }
986 : }
987 0 : else if (EQUAL(pszProjection, SRS_PT_VANDERGRINTEN))
988 : {
989 0 : CPLPrintStringFill(szProj, "VDG", knProjSize);
990 0 : (*ppadfPrjParams)[2] = GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
991 0 : (*ppadfPrjParams)[6] = GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
992 0 : (*ppadfPrjParams)[7] = GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
993 : }
994 : // Projection unsupported by PCI.
995 : else
996 : {
997 0 : CPLDebug("OSR_PCI",
998 : "Projection \"%s\" unsupported by PCI. "
999 : "PIXEL value will be used.",
1000 : pszProjection);
1001 0 : CPLPrintStringFill(szProj, "PIXEL", knProjSize);
1002 : }
1003 :
1004 : /* ==================================================================== */
1005 : /* Translate the earth model. */
1006 : /* ==================================================================== */
1007 :
1008 : /* -------------------------------------------------------------------- */
1009 : /* Is this a well known datum? */
1010 : /* -------------------------------------------------------------------- */
1011 69 : const char *pszDatum = GetAttrValue("DATUM");
1012 69 : char szEarthModel[5] = {};
1013 :
1014 69 : if (pszDatum == nullptr || strlen(pszDatum) == 0)
1015 : {
1016 : // Do nothing.
1017 : }
1018 67 : else if (EQUAL(pszDatum, SRS_DN_NAD27))
1019 : {
1020 6 : CPLPrintStringFill(szEarthModel, "D-01", 4);
1021 : }
1022 61 : else if (EQUAL(pszDatum, SRS_DN_NAD83))
1023 : {
1024 0 : CPLPrintStringFill(szEarthModel, "D-02", 4);
1025 : }
1026 61 : else if (EQUAL(pszDatum, SRS_DN_WGS84))
1027 : {
1028 53 : CPLPrintStringFill(szEarthModel, "D000", 4);
1029 53 : if (pszProjection && EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
1030 : {
1031 : // Special case for Web Mercator
1032 1 : const char *method_name = "";
1033 1 : GetWKT2ProjectionMethod(&method_name, nullptr, nullptr);
1034 1 : if (EQUAL(method_name, "Popular Visualisation Pseudo Mercator"))
1035 : {
1036 1 : CPLPrintStringFill(szEarthModel, "D894", 4);
1037 : }
1038 : }
1039 : }
1040 :
1041 : /* -------------------------------------------------------------------- */
1042 : /* If not a very well known datum, try for an EPSG based */
1043 : /* translation. */
1044 : /* -------------------------------------------------------------------- */
1045 69 : if (szEarthModel[0] == '\0')
1046 : {
1047 10 : const char *pszAuthority = GetAuthorityName("GEOGCS");
1048 :
1049 10 : if (pszAuthority && EQUAL(pszAuthority, "EPSG"))
1050 : {
1051 1 : const int nGCS_EPSG = atoi(GetAuthorityCode("GEOGCS"));
1052 :
1053 10 : for (int i = 0; asDatums[i].nEPSGCode != 0; i++)
1054 : {
1055 10 : if (asDatums[i].nEPSGCode == nGCS_EPSG)
1056 : {
1057 1 : snprintf(szEarthModel, sizeof(szEarthModel), "%s",
1058 1 : asDatums[i].pszPCIDatum);
1059 1 : break;
1060 : }
1061 : }
1062 : }
1063 : }
1064 :
1065 : /* -------------------------------------------------------------------- */
1066 : /* If we haven't found something yet, try translating the */
1067 : /* ellipsoid. */
1068 : /* -------------------------------------------------------------------- */
1069 69 : const double dfSemiMajor = GetSemiMajor();
1070 69 : const double dfInvFlattening = GetInvFlattening();
1071 69 : if (szEarthModel[0] == '\0')
1072 : {
1073 9 : const PCIDatums *pasDatum = asEllips;
1074 :
1075 137 : while (pasDatum->pszPCIDatum)
1076 : {
1077 133 : double dfSM = 0.0;
1078 133 : double dfIF = 0.0;
1079 :
1080 133 : if (OSRGetEllipsoidInfo(pasDatum->nEPSGCode, nullptr, &dfSM,
1081 133 : &dfIF) == OGRERR_NONE &&
1082 143 : CPLIsEqual(dfSemiMajor, dfSM) &&
1083 10 : CPLIsEqual(dfInvFlattening, dfIF))
1084 : {
1085 5 : CPLPrintStringFill(szEarthModel, pasDatum->pszPCIDatum, 4);
1086 5 : break;
1087 : }
1088 :
1089 128 : pasDatum++;
1090 : }
1091 : }
1092 :
1093 : // Try to find in pci_ellips.txt.
1094 69 : if (szEarthModel[0] == '\0')
1095 : {
1096 4 : const char *pszCSV = CSVFilename("pci_ellips.txt");
1097 : const double dfSemiMinor =
1098 4 : OSRCalcSemiMinorFromInvFlattening(dfSemiMajor, dfInvFlattening);
1099 :
1100 4 : VSILFILE *fp = pszCSV ? VSIFOpenL(pszCSV, "r") : nullptr;
1101 :
1102 4 : if (fp != nullptr)
1103 : {
1104 4 : char **papszLineItems = nullptr;
1105 :
1106 194 : while ((papszLineItems = CSVReadParseLineL(fp)) != nullptr)
1107 : {
1108 194 : if (CSLCount(papszLineItems) >= 4 &&
1109 198 : CPLIsEqual(dfSemiMajor, CPLAtof(papszLineItems[2])) &&
1110 4 : CPLIsEqual(dfSemiMinor, CPLAtof(papszLineItems[3])))
1111 : {
1112 4 : snprintf(szEarthModel, sizeof(szEarthModel), "%s",
1113 : papszLineItems[0]);
1114 4 : break;
1115 : }
1116 :
1117 190 : CSLDestroy(papszLineItems);
1118 : }
1119 :
1120 4 : CSLDestroy(papszLineItems);
1121 4 : VSIFCloseL(fp);
1122 : }
1123 : }
1124 :
1125 : // Custom ellipsoid parameters.
1126 69 : if (szEarthModel[0] == '\0')
1127 : {
1128 0 : CPLPrintStringFill(szEarthModel, "E999", 4);
1129 0 : (*ppadfPrjParams)[0] = dfSemiMajor;
1130 0 : (*ppadfPrjParams)[1] =
1131 0 : OSRCalcSemiMinorFromInvFlattening(dfSemiMajor, dfInvFlattening);
1132 : }
1133 :
1134 : /* -------------------------------------------------------------------- */
1135 : /* If we have a non-parameteric ellipsoid, scan the */
1136 : /* pci_datum.txt for a match. */
1137 : /* -------------------------------------------------------------------- */
1138 69 : if (szEarthModel[0] == 'E' && !EQUAL(szEarthModel, "E999") &&
1139 : pszDatum != nullptr)
1140 : {
1141 7 : const char *pszDatumCSV = CSVFilename("pci_datum.txt");
1142 7 : double adfTOWGS84[7] = {};
1143 7 : const bool bHaveTOWGS84 = GetTOWGS84(adfTOWGS84, 7) == OGRERR_NONE;
1144 :
1145 7 : VSILFILE *fp = pszDatumCSV ? VSIFOpenL(pszDatumCSV, "r") : nullptr;
1146 :
1147 7 : if (fp != nullptr)
1148 : {
1149 7 : char **papszLineItems = nullptr;
1150 :
1151 3285 : while ((papszLineItems = CSVReadParseLineL(fp)) != nullptr)
1152 : {
1153 : // Compare based on datum name. This is mostly for
1154 : // PCI round-tripping. We won't usually get exact matches
1155 : // from other sources.
1156 3280 : if (CSLCount(papszLineItems) > 3 &&
1157 3281 : EQUAL(papszLineItems[1], pszDatum) &&
1158 1 : EQUAL(papszLineItems[2], szEarthModel))
1159 : {
1160 1 : snprintf(szEarthModel, sizeof(szEarthModel), "%s",
1161 : papszLineItems[0]);
1162 1 : break;
1163 : }
1164 :
1165 3279 : bool bTOWGS84Match = bHaveTOWGS84;
1166 :
1167 3279 : if (CSLCount(papszLineItems) < 11)
1168 848 : bTOWGS84Match = false;
1169 :
1170 3780 : if (bTOWGS84Match &&
1171 501 : (!CPLIsEqual(adfTOWGS84[0], CPLAtof(papszLineItems[3])) ||
1172 1 : !CPLIsEqual(adfTOWGS84[1], CPLAtof(papszLineItems[4])) ||
1173 1 : !CPLIsEqual(adfTOWGS84[2], CPLAtof(papszLineItems[5]))))
1174 500 : bTOWGS84Match = false;
1175 :
1176 3280 : if (bTOWGS84Match && CSLCount(papszLineItems) >= 15 &&
1177 1 : (!CPLIsEqual(adfTOWGS84[3], CPLAtof(papszLineItems[11])) ||
1178 1 : !CPLIsEqual(adfTOWGS84[4], CPLAtof(papszLineItems[12])) ||
1179 1 : !CPLIsEqual(adfTOWGS84[5], CPLAtof(papszLineItems[13]))))
1180 0 : bTOWGS84Match = false;
1181 :
1182 3279 : if (bTOWGS84Match && CSLCount(papszLineItems) >= 15)
1183 : {
1184 1 : double dfScale = CPLAtof(papszLineItems[14]);
1185 :
1186 : // Convert to parts per million if is a 1 based scaling.
1187 1 : if (dfScale >= 0.999 && dfScale <= 1.001)
1188 1 : dfScale = (dfScale - 1.0) * 1000000.0;
1189 :
1190 1 : if (fabs(adfTOWGS84[6] - dfScale) >
1191 1 : 1e-11 * fabs(adfTOWGS84[6]))
1192 0 : bTOWGS84Match = false;
1193 : }
1194 :
1195 3279 : if (bTOWGS84Match && CSLCount(papszLineItems) < 15 &&
1196 0 : (!CPLIsEqual(adfTOWGS84[3], 0.0) ||
1197 0 : !CPLIsEqual(adfTOWGS84[4], 0.0) ||
1198 0 : !CPLIsEqual(adfTOWGS84[5], 0.0) ||
1199 0 : !CPLIsEqual(adfTOWGS84[6], 0.0)))
1200 0 : bTOWGS84Match = false;
1201 :
1202 3279 : if (bTOWGS84Match)
1203 : {
1204 1 : snprintf(szEarthModel, sizeof(szEarthModel), "%s",
1205 : papszLineItems[0]);
1206 1 : break;
1207 : }
1208 :
1209 3278 : CSLDestroy(papszLineItems);
1210 : }
1211 :
1212 7 : CSLDestroy(papszLineItems);
1213 7 : VSIFCloseL(fp);
1214 : }
1215 : }
1216 :
1217 69 : CPLPrintStringFill(szProj + 12, szEarthModel, 4);
1218 :
1219 69 : CPLDebug("OSR_PCI", "Translated as '%s'", szProj);
1220 :
1221 : /* -------------------------------------------------------------------- */
1222 : /* Translate the linear units. */
1223 : /* -------------------------------------------------------------------- */
1224 69 : const char *pszUnits =
1225 69 : STARTS_WITH_CI(szProj, "LONG/LAT") ? "DEGREE" : "METRE";
1226 :
1227 : /* -------------------------------------------------------------------- */
1228 : /* Report results. */
1229 : /* -------------------------------------------------------------------- */
1230 69 : szProj[knProjSize] = '\0';
1231 69 : *ppszProj = CPLStrdup(szProj);
1232 :
1233 69 : *ppszUnits = CPLStrdup(pszUnits);
1234 :
1235 69 : return OGRERR_NONE;
1236 : }
|
