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