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