Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: OGRSpatialReference translation to/from ESRI .prj definitions.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2000, Frank Warmerdam
9 : * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2013, Kyle Shannon <kyle at pobox dot com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "cpl_port.h"
16 : #include "ogr_spatialref.h"
17 : #include "ogr_srs_esri_names.h"
18 :
19 : #include <cmath>
20 : #include <climits>
21 : #include <cstddef>
22 : #include <cstdio>
23 : #include <cstdlib>
24 : #include <cstring>
25 : #include <algorithm>
26 : #include <limits>
27 : #include <string>
28 :
29 : #include "cpl_conv.h"
30 : #include "cpl_csv.h"
31 : #include "cpl_error.h"
32 : #include "cpl_multiproc.h"
33 : #include "cpl_string.h"
34 : #include "cpl_vsi.h"
35 : #include "ogr_core.h"
36 : #include "ogr_p.h"
37 : #include "ogr_srs_api.h"
38 :
39 : /* -------------------------------------------------------------------- */
40 : /* Table relating USGS and ESRI state plane zones. */
41 : /* -------------------------------------------------------------------- */
42 : constexpr int anUsgsEsriZones[] = {
43 : 101, 3101, 102, 3126, 201, 3151, 202, 3176, 203, 3201, 301, 3226,
44 : 302, 3251, 401, 3276, 402, 3301, 403, 3326, 404, 3351, 405, 3376,
45 : 406, 3401, 407, 3426, 501, 3451, 502, 3476, 503, 3501, 600, 3526,
46 : 700, 3551, 901, 3601, 902, 3626, 903, 3576, 1001, 3651, 1002, 3676,
47 : 1101, 3701, 1102, 3726, 1103, 3751, 1201, 3776, 1202, 3801, 1301, 3826,
48 : 1302, 3851, 1401, 3876, 1402, 3901, 1501, 3926, 1502, 3951, 1601, 3976,
49 : 1602, 4001, 1701, 4026, 1702, 4051, 1703, 6426, 1801, 4076, 1802, 4101,
50 : 1900, 4126, 2001, 4151, 2002, 4176, 2101, 4201, 2102, 4226, 2103, 4251,
51 : 2111, 6351, 2112, 6376, 2113, 6401, 2201, 4276, 2202, 4301, 2203, 4326,
52 : 2301, 4351, 2302, 4376, 2401, 4401, 2402, 4426, 2403, 4451, 2500, 0,
53 : 2501, 4476, 2502, 4501, 2503, 4526, 2600, 0, 2601, 4551, 2602, 4576,
54 : 2701, 4601, 2702, 4626, 2703, 4651, 2800, 4676, 2900, 4701, 3001, 4726,
55 : 3002, 4751, 3003, 4776, 3101, 4801, 3102, 4826, 3103, 4851, 3104, 4876,
56 : 3200, 4901, 3301, 4926, 3302, 4951, 3401, 4976, 3402, 5001, 3501, 5026,
57 : 3502, 5051, 3601, 5076, 3602, 5101, 3701, 5126, 3702, 5151, 3800, 5176,
58 : 3900, 0, 3901, 5201, 3902, 5226, 4001, 5251, 4002, 5276, 4100, 5301,
59 : 4201, 5326, 4202, 5351, 4203, 5376, 4204, 5401, 4205, 5426, 4301, 5451,
60 : 4302, 5476, 4303, 5501, 4400, 5526, 4501, 5551, 4502, 5576, 4601, 5601,
61 : 4602, 5626, 4701, 5651, 4702, 5676, 4801, 5701, 4802, 5726, 4803, 5751,
62 : 4901, 5776, 4902, 5801, 4903, 5826, 4904, 5851, 5001, 6101, 5002, 6126,
63 : 5003, 6151, 5004, 6176, 5005, 6201, 5006, 6226, 5007, 6251, 5008, 6276,
64 : 5009, 6301, 5010, 6326, 5101, 5876, 5102, 5901, 5103, 5926, 5104, 5951,
65 : 5105, 5976, 5201, 6001, 5200, 6026, 5200, 6076, 5201, 6051, 5202, 6051,
66 : 5300, 0, 5400, 0};
67 :
68 : /************************************************************************/
69 : /* ESRIToUSGSZone() */
70 : /* */
71 : /* Convert ESRI style state plane zones to USGS style state */
72 : /* plane zones. */
73 : /************************************************************************/
74 :
75 0 : static int ESRIToUSGSZone(int nESRIZone)
76 :
77 : {
78 : // anUsgsEsriZones is a series of ints where 2 consecutive integers
79 : // are used to map from USGS to ESRI state plane zones.
80 : // TODO(schwehr): Would be better as a std::map.
81 0 : const int nPairs = sizeof(anUsgsEsriZones) / (2 * sizeof(int));
82 :
83 0 : for (int i = 0; i < nPairs; i++)
84 : {
85 0 : if (anUsgsEsriZones[i * 2 + 1] == nESRIZone)
86 0 : return anUsgsEsriZones[i * 2];
87 : }
88 :
89 0 : return 0;
90 : }
91 :
92 : /************************************************************************/
93 : /* OSRImportFromESRI() */
94 : /************************************************************************/
95 :
96 : /**
97 : * \brief Import coordinate system from ESRI .prj format(s).
98 : *
99 : * This function is the same as the C++ method
100 : * OGRSpatialReference::importFromESRI().
101 : */
102 11 : OGRErr OSRImportFromESRI(OGRSpatialReferenceH hSRS, char **papszPrj)
103 :
104 : {
105 11 : VALIDATE_POINTER1(hSRS, "OSRImportFromESRI", OGRERR_FAILURE);
106 :
107 11 : return reinterpret_cast<OGRSpatialReference *>(hSRS)->importFromESRI(
108 11 : papszPrj);
109 : }
110 :
111 : /************************************************************************/
112 : /* OSR_GDV() */
113 : /* */
114 : /* Fetch a particular parameter out of the parameter list, or */
115 : /* the indicated default if it isn't available. This is a */
116 : /* helper function for importFromESRI(). */
117 : /************************************************************************/
118 :
119 134 : static double OSR_GDV(char **papszNV, const char *pszField,
120 : double dfDefaultValue)
121 :
122 : {
123 134 : if (papszNV == nullptr || papszNV[0] == nullptr)
124 0 : return dfDefaultValue;
125 :
126 134 : if (STARTS_WITH_CI(pszField, "PARAM_"))
127 : {
128 102 : int iLine = 0; // Used after for loop.
129 790 : for (; papszNV[iLine] != nullptr &&
130 790 : !STARTS_WITH_CI(papszNV[iLine], "Paramet");
131 : iLine++)
132 : {
133 : }
134 :
135 441 : for (int nOffset = atoi(pszField + 6);
136 441 : papszNV[iLine] != nullptr && nOffset > 0; iLine++)
137 : {
138 339 : if (strlen(papszNV[iLine]) > 0)
139 339 : nOffset--;
140 : }
141 :
142 102 : while (papszNV[iLine] != nullptr && strlen(papszNV[iLine]) == 0)
143 0 : iLine++;
144 :
145 102 : if (papszNV[iLine] != nullptr)
146 : {
147 102 : char *const pszLine = papszNV[iLine];
148 :
149 : // Trim comments.
150 7027 : for (int i = 0; pszLine[i] != '\0'; i++)
151 : {
152 6925 : if (pszLine[i] == '/' && pszLine[i + 1] == '*')
153 101 : pszLine[i] = '\0';
154 : }
155 :
156 102 : double dfValue = 0.0;
157 102 : char **papszTokens = CSLTokenizeString(papszNV[iLine]);
158 102 : if (CSLCount(papszTokens) == 3)
159 : {
160 : // http://agdcftp1.wr.usgs.gov/pub/projects/lcc/akcan_lcc/akcan.tar.gz
161 : // contains weird values for the second. Ignore it and
162 : // the result looks correct.
163 62 : double dfSecond = CPLAtof(papszTokens[2]);
164 62 : if (dfSecond < 0.0 || dfSecond >= 60.0)
165 0 : dfSecond = 0.0;
166 :
167 62 : dfValue = std::abs(CPLAtof(papszTokens[0])) +
168 62 : CPLAtof(papszTokens[1]) / 60.0 + dfSecond / 3600.0;
169 :
170 62 : if (CPLAtof(papszTokens[0]) < 0.0)
171 15 : dfValue *= -1;
172 : }
173 40 : else if (CSLCount(papszTokens) > 0)
174 : {
175 40 : dfValue = CPLAtof(papszTokens[0]);
176 : }
177 : else
178 : {
179 0 : dfValue = dfDefaultValue;
180 : }
181 :
182 102 : CSLDestroy(papszTokens);
183 :
184 102 : return dfValue;
185 : }
186 :
187 0 : return dfDefaultValue;
188 : }
189 :
190 32 : int iLine = 0; // Used after for loop.
191 156 : for (; papszNV[iLine] != nullptr &&
192 152 : !EQUALN(papszNV[iLine], pszField, strlen(pszField));
193 : iLine++)
194 : {
195 : }
196 :
197 32 : if (papszNV[iLine] == nullptr)
198 4 : return dfDefaultValue;
199 :
200 28 : return CPLAtof(papszNV[iLine] + strlen(pszField));
201 : }
202 :
203 : /************************************************************************/
204 : /* OSR_GDS() */
205 : /************************************************************************/
206 :
207 117 : static CPLString OSR_GDS(char **papszNV, const char *pszField,
208 : const char *pszDefaultValue)
209 :
210 : {
211 117 : if (papszNV == nullptr || papszNV[0] == nullptr)
212 0 : return pszDefaultValue;
213 :
214 117 : int iLine = 0; // Used after for loop.
215 370 : for (; papszNV[iLine] != nullptr &&
216 359 : !EQUALN(papszNV[iLine], pszField, strlen(pszField));
217 : iLine++)
218 : {
219 : }
220 :
221 117 : if (papszNV[iLine] == nullptr)
222 11 : return pszDefaultValue;
223 :
224 106 : char **papszTokens = CSLTokenizeString(papszNV[iLine]);
225 :
226 : CPLString osResult =
227 212 : CSLCount(papszTokens) > 1 ? papszTokens[1] : pszDefaultValue;
228 :
229 106 : CSLDestroy(papszTokens);
230 106 : return osResult;
231 : }
232 :
233 : /************************************************************************/
234 : /* importFromESRI() */
235 : /************************************************************************/
236 :
237 : /**
238 : * \brief Import coordinate system from ESRI .prj format(s).
239 : *
240 : * This function will read the text loaded from an ESRI .prj file, and
241 : * translate it into an OGRSpatialReference definition. This should support
242 : * many (but by no means all) old style (Arc/Info 7.x) .prj files, as well
243 : * as the newer pseudo-OGC WKT .prj files. Note that new style .prj files
244 : * are in OGC WKT format, but require some manipulation to correct datum
245 : * names, and units on some projection parameters. This is addressed within
246 : * importFromESRI() by an automatic call to morphFromESRI().
247 : *
248 : * Currently only GEOGRAPHIC, UTM, STATEPLANE, GREATBRITIAN_GRID, ALBERS,
249 : * EQUIDISTANT_CONIC, TRANSVERSE (mercator), POLAR, LAMBERT (Conic Conformal),
250 : * LAMBERT_AZIMUTHAL, MERCATOR and POLYCONIC
251 : * projections are supported from old style files.
252 : *
253 : * At this time there is no equivalent exportToESRI() method. Writing old
254 : * style .prj files is not supported by OGRSpatialReference. However the
255 : * morphToESRI() and exportToWkt() methods can be used to generate output
256 : * suitable to write to new style (Arc 8) .prj files.
257 : *
258 : * This function is the equivalent of the C function OSRImportFromESRI().
259 : *
260 : * @param papszPrj NULL terminated list of strings containing the definition.
261 : *
262 : * @return OGRERR_NONE on success or an error code in case of failure.
263 : */
264 :
265 943 : OGRErr OGRSpatialReference::importFromESRI(char **papszPrj)
266 :
267 : {
268 943 : if (papszPrj == nullptr || papszPrj[0] == nullptr)
269 29 : return OGRERR_CORRUPT_DATA;
270 :
271 : /* -------------------------------------------------------------------- */
272 : /* ArcGIS and related products now use a variant of Well Known */
273 : /* Text. Try to recognize this and ingest it. WKT is usually */
274 : /* all on one line, but we will accept multi-line formats and */
275 : /* concatenate. */
276 : /* -------------------------------------------------------------------- */
277 914 : if (STARTS_WITH_CI(papszPrj[0], "GEOGCS") ||
278 563 : STARTS_WITH_CI(papszPrj[0], "PROJCS") ||
279 43 : STARTS_WITH_CI(papszPrj[0], "LOCAL_CS")
280 : // Also accept COMPD_CS, even if it is unclear that it is valid
281 : // traditional ESRI WKT. But people might use such PRJ file
282 : // See https://github.com/OSGeo/gdal/issues/1881
283 43 : || STARTS_WITH_CI(papszPrj[0], "COMPD_CS"))
284 : {
285 1742 : std::string osWKT(papszPrj[0]);
286 873 : for (int i = 1; papszPrj[i] != nullptr; i++)
287 : {
288 2 : osWKT += papszPrj[i];
289 : }
290 871 : return importFromWkt(osWKT.c_str());
291 : }
292 :
293 : /* -------------------------------------------------------------------- */
294 : /* Operate on the basis of the projection name. */
295 : /* -------------------------------------------------------------------- */
296 86 : CPLString osProj = OSR_GDS(papszPrj, "Projection", "");
297 43 : bool bDatumApplied = false;
298 43 : bool bHasRadiusOfSphereOfReference = false;
299 43 : double dfRadiusOfSphereOfReference = 0.0;
300 :
301 43 : if (EQUAL(osProj, ""))
302 : {
303 5 : CPLDebug("OGR_ESRI", "Can't find Projection");
304 5 : return OGRERR_CORRUPT_DATA;
305 : }
306 38 : else if (EQUAL(osProj, "GEOGRAPHIC"))
307 : {
308 : // Nothing to do.
309 : }
310 35 : else if (EQUAL(osProj, "utm"))
311 : {
312 12 : const double dfOsrGdv = OSR_GDV(papszPrj, "zone", 0.0);
313 12 : if (dfOsrGdv > 0 && dfOsrGdv < 61)
314 : {
315 12 : const double dfYShift = OSR_GDV(papszPrj, "Yshift", 0.0);
316 :
317 12 : SetUTM(static_cast<int>(dfOsrGdv), dfYShift == 0.0);
318 : }
319 : else
320 : {
321 0 : const double dfCentralMeridian = OSR_GDV(papszPrj, "PARAM_1", 0.0);
322 0 : const double dfRefLat = OSR_GDV(papszPrj, "PARAM_2", 0.0);
323 0 : if (dfCentralMeridian >= -180.0 && dfCentralMeridian <= 180.0)
324 : {
325 0 : const int nZone = static_cast<int>(
326 0 : (dfCentralMeridian + 183.0) / 6.0 + 0.0000001);
327 0 : SetUTM(nZone, dfRefLat >= 0.0);
328 : }
329 : }
330 : }
331 23 : else if (EQUAL(osProj, "STATEPLANE"))
332 : {
333 4 : const double dfZone = OSR_GDV(papszPrj, "zone", 0.0);
334 :
335 4 : if (dfZone < std::numeric_limits<int>::min() ||
336 4 : dfZone > std::numeric_limits<int>::max() || std::isnan(dfZone))
337 : {
338 0 : CPLError(CE_Failure, CPLE_AppDefined, "zone out of range: %f",
339 : dfZone);
340 0 : return OGRERR_CORRUPT_DATA;
341 : }
342 :
343 4 : int nZone = static_cast<int>(dfZone);
344 :
345 4 : if (nZone != 0)
346 0 : nZone = ESRIToUSGSZone(nZone);
347 : else
348 : {
349 4 : const double dfFipszone = OSR_GDV(papszPrj, "fipszone", 0.0);
350 :
351 4 : if (dfFipszone < std::numeric_limits<int>::min() ||
352 8 : dfFipszone > std::numeric_limits<int>::max() ||
353 4 : std::isnan(dfFipszone))
354 : {
355 0 : CPLError(CE_Failure, CPLE_AppDefined,
356 : "fipszone out of range: %f", dfFipszone);
357 0 : return OGRERR_CORRUPT_DATA;
358 : }
359 :
360 4 : nZone = static_cast<int>(dfFipszone);
361 : }
362 :
363 4 : if (nZone != 0)
364 : {
365 4 : bDatumApplied = true;
366 4 : if (EQUAL(OSR_GDS(papszPrj, "Datum", "NAD83"), "NAD27"))
367 0 : SetStatePlane(nZone, FALSE);
368 : else
369 4 : SetStatePlane(nZone, TRUE);
370 : }
371 : }
372 38 : else if (EQUAL(osProj, "GREATBRITIAN_GRID") ||
373 19 : EQUAL(osProj, "GREATBRITAIN_GRID"))
374 : {
375 0 : const char *pszWkt =
376 : "PROJCS[\"OSGB 1936 / British National Grid\","
377 : "GEOGCS[\"OSGB 1936\",DATUM[\"OSGB_1936\","
378 : "SPHEROID[\"Airy 1830\",6377563.396,299.3249646]],"
379 : "PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],"
380 : "PROJECTION[\"Transverse_Mercator\"],"
381 : "PARAMETER[\"latitude_of_origin\",49],"
382 : "PARAMETER[\"central_meridian\",-2],"
383 : "PARAMETER[\"scale_factor\",0.999601272],"
384 : "PARAMETER[\"false_easting\",400000],"
385 : "PARAMETER[\"false_northing\",-100000],UNIT[\"metre\",1]]";
386 :
387 0 : bDatumApplied = true;
388 0 : importFromWkt(pszWkt);
389 : }
390 19 : else if (EQUAL(osProj, "ALBERS"))
391 : {
392 13 : SetACEA(OSR_GDV(papszPrj, "PARAM_1", 0.0),
393 : OSR_GDV(papszPrj, "PARAM_2", 0.0),
394 : OSR_GDV(papszPrj, "PARAM_4", 0.0),
395 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
396 : OSR_GDV(papszPrj, "PARAM_5", 0.0),
397 : OSR_GDV(papszPrj, "PARAM_6", 0.0));
398 : }
399 6 : else if (EQUAL(osProj, "LAMBERT"))
400 : {
401 0 : SetLCC(OSR_GDV(papszPrj, "PARAM_1", 0.0),
402 : OSR_GDV(papszPrj, "PARAM_2", 0.0),
403 : OSR_GDV(papszPrj, "PARAM_4", 0.0),
404 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
405 : OSR_GDV(papszPrj, "PARAM_5", 0.0),
406 : OSR_GDV(papszPrj, "PARAM_6", 0.0));
407 : }
408 6 : else if (EQUAL(osProj, "LAMBERT_AZIMUTHAL"))
409 : {
410 27 : for (int iLine = 0; papszPrj[iLine] != nullptr; iLine++)
411 : {
412 26 : if (strstr(papszPrj[iLine], "radius of the sphere of reference"))
413 : {
414 2 : bHasRadiusOfSphereOfReference = true;
415 2 : break;
416 : }
417 : }
418 3 : if (bHasRadiusOfSphereOfReference)
419 : {
420 : // Cf "Workstation" variation of
421 : // https://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Lambert_Azimuthal_Equal_Area
422 : // that is supposed to be applied to spheres only.
423 : // It is also documented at page 71 of
424 : // https://kartoweb.itc.nl/geometrics/Map%20projections/Understanding%20Map%20Projections.pdf
425 : // ("Understanding Map Projections", Melita Kennedy, ArcInfo 8)
426 : // We don't particularly enforce the restriction to spheres, so if
427 : // the "Parameters" line had non-spherical axis dimensions, they
428 : // would be used.
429 : // EPSG has a EPSG:1027 projection method "Lambert Azimuthal Equal Area (Spherical)"
430 : // that uses the radius of the authalic sphere, for non-spherical
431 : // ellipsoids, but it is not obvious that it would be appropriate here.
432 : // Examples:
433 : // - https://community.esri.com/t5/data-management-questions/problem-when-projecting-a-raster-file/td-p/400814/page/2
434 : // - https://lists.osgeo.org/pipermail/gdal-dev/2024-May/058990.html
435 2 : dfRadiusOfSphereOfReference = OSR_GDV(papszPrj, "PARAM_1", 0.0);
436 2 : SetLAEA(OSR_GDV(papszPrj, "PARAM_3", 0.0),
437 : OSR_GDV(papszPrj, "PARAM_2", 0.0),
438 : OSR_GDV(papszPrj, "PARAM_4", 0.0),
439 : OSR_GDV(papszPrj, "PARAM_5", 0.0));
440 : }
441 : else
442 : {
443 : // Example: https://trac.osgeo.org/gdal/ticket/4302
444 1 : SetLAEA(OSR_GDV(papszPrj, "PARAM_2", 0.0),
445 : OSR_GDV(papszPrj, "PARAM_1", 0.0),
446 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
447 : OSR_GDV(papszPrj, "PARAM_4", 0.0));
448 : }
449 : }
450 3 : else if (EQUAL(osProj, "EQUIDISTANT_CONIC"))
451 : {
452 1 : const double dfStdPCount = OSR_GDV(papszPrj, "PARAM_1", 0.0);
453 : // TODO(schwehr): What is a reasonable range for StdPCount?
454 1 : if (dfStdPCount < 0 || dfStdPCount > std::numeric_limits<int>::max() ||
455 0 : std::isnan(dfStdPCount))
456 : {
457 1 : CPLError(CE_Failure, CPLE_AppDefined, "StdPCount out of range: %lf",
458 : dfStdPCount);
459 1 : return OGRERR_CORRUPT_DATA;
460 : }
461 0 : const int nStdPCount = static_cast<int>(dfStdPCount);
462 :
463 0 : if (nStdPCount == 1)
464 : {
465 0 : SetEC(OSR_GDV(papszPrj, "PARAM_2", 0.0),
466 : OSR_GDV(papszPrj, "PARAM_2", 0.0),
467 : OSR_GDV(papszPrj, "PARAM_4", 0.0),
468 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
469 : OSR_GDV(papszPrj, "PARAM_5", 0.0),
470 : OSR_GDV(papszPrj, "PARAM_6", 0.0));
471 : }
472 : else
473 : {
474 0 : SetEC(OSR_GDV(papszPrj, "PARAM_2", 0.0),
475 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
476 : OSR_GDV(papszPrj, "PARAM_5", 0.0),
477 : OSR_GDV(papszPrj, "PARAM_4", 0.0),
478 : OSR_GDV(papszPrj, "PARAM_5", 0.0),
479 : OSR_GDV(papszPrj, "PARAM_7", 0.0));
480 : }
481 : }
482 2 : else if (EQUAL(osProj, "TRANSVERSE"))
483 : {
484 1 : SetTM(OSR_GDV(papszPrj, "PARAM_3", 0.0),
485 : OSR_GDV(papszPrj, "PARAM_2", 0.0),
486 : OSR_GDV(papszPrj, "PARAM_1", 0.0),
487 : OSR_GDV(papszPrj, "PARAM_4", 0.0),
488 : OSR_GDV(papszPrj, "PARAM_5", 0.0));
489 : }
490 1 : else if (EQUAL(osProj, "POLAR"))
491 : {
492 0 : SetPS(OSR_GDV(papszPrj, "PARAM_2", 0.0),
493 : OSR_GDV(papszPrj, "PARAM_1", 0.0), 1.0,
494 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
495 : OSR_GDV(papszPrj, "PARAM_4", 0.0));
496 : }
497 1 : else if (EQUAL(osProj, "MERCATOR"))
498 : {
499 1 : SetMercator2SP(OSR_GDV(papszPrj, "PARAM_2", 0.0), 0.0,
500 : OSR_GDV(papszPrj, "PARAM_1", 0.0),
501 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
502 : OSR_GDV(papszPrj, "PARAM_4", 0.0));
503 : }
504 0 : else if (EQUAL(osProj, SRS_PT_MERCATOR_AUXILIARY_SPHERE))
505 : {
506 : // This is EPSG:3875 Pseudo Mercator. We might as well import it from
507 : // the EPSG spec.
508 0 : importFromEPSG(3857);
509 0 : bDatumApplied = true;
510 : }
511 0 : else if (EQUAL(osProj, "POLYCONIC"))
512 : {
513 0 : SetPolyconic(OSR_GDV(papszPrj, "PARAM_2", 0.0),
514 : OSR_GDV(papszPrj, "PARAM_1", 0.0),
515 : OSR_GDV(papszPrj, "PARAM_3", 0.0),
516 : OSR_GDV(papszPrj, "PARAM_4", 0.0));
517 : }
518 : else
519 : {
520 0 : CPLDebug("OGR_ESRI", "Unsupported projection: %s", osProj.c_str());
521 0 : SetLocalCS(osProj);
522 : }
523 :
524 : /* -------------------------------------------------------------------- */
525 : /* Try to translate the datum/spheroid. */
526 : /* -------------------------------------------------------------------- */
527 37 : if (!IsLocal() && !bDatumApplied)
528 : {
529 66 : const CPLString osDatum = OSR_GDS(papszPrj, "Datum", "");
530 :
531 61 : if (EQUAL(osDatum, "NAD27") || EQUAL(osDatum, "NAD83") ||
532 61 : EQUAL(osDatum, "WGS84") || EQUAL(osDatum, "WGS72"))
533 : {
534 21 : SetWellKnownGeogCS(osDatum);
535 : }
536 12 : else if (EQUAL(osDatum, "EUR") || EQUAL(osDatum, "ED50"))
537 : {
538 0 : SetWellKnownGeogCS("EPSG:4230");
539 : }
540 12 : else if (EQUAL(osDatum, "GDA94"))
541 : {
542 9 : SetWellKnownGeogCS("EPSG:4283");
543 : }
544 : else
545 : {
546 6 : CPLString osSpheroid = OSR_GDS(papszPrj, "Spheroid", "");
547 :
548 6 : if (EQUAL(osSpheroid, "INT1909") ||
549 3 : EQUAL(osSpheroid, "INTERNATIONAL1909"))
550 : {
551 0 : OGRSpatialReference oGCS;
552 0 : oGCS.importFromEPSG(4022);
553 0 : CopyGeogCSFrom(&oGCS);
554 : }
555 3 : else if (EQUAL(osSpheroid, "AIRY"))
556 : {
557 0 : OGRSpatialReference oGCS;
558 0 : oGCS.importFromEPSG(4001);
559 0 : CopyGeogCSFrom(&oGCS);
560 : }
561 3 : else if (EQUAL(osSpheroid, "CLARKE1866"))
562 : {
563 0 : OGRSpatialReference oGCS;
564 0 : oGCS.importFromEPSG(4008);
565 0 : CopyGeogCSFrom(&oGCS);
566 : }
567 3 : else if (EQUAL(osSpheroid, "GRS80"))
568 : {
569 0 : OGRSpatialReference oGCS;
570 0 : oGCS.importFromEPSG(4019);
571 0 : CopyGeogCSFrom(&oGCS);
572 : }
573 3 : else if (EQUAL(osSpheroid, "KRASOVSKY") ||
574 6 : EQUAL(osSpheroid, "KRASSOVSKY") ||
575 3 : EQUAL(osSpheroid, "KRASSOWSKY"))
576 : {
577 0 : OGRSpatialReference oGCS;
578 0 : oGCS.importFromEPSG(4024);
579 0 : CopyGeogCSFrom(&oGCS);
580 : }
581 3 : else if (EQUAL(osSpheroid, "Bessel"))
582 : {
583 0 : OGRSpatialReference oGCS;
584 0 : oGCS.importFromEPSG(4004);
585 0 : CopyGeogCSFrom(&oGCS);
586 : }
587 : else
588 : {
589 3 : bool bFoundParameters = false;
590 14 : for (int iLine = 0; papszPrj[iLine] != nullptr; iLine++)
591 : {
592 14 : if (STARTS_WITH_CI(papszPrj[iLine], "Parameters"))
593 : {
594 6 : char **papszTokens = CSLTokenizeString(
595 3 : papszPrj[iLine] + strlen("Parameters"));
596 3 : if (CSLCount(papszTokens) == 2)
597 : {
598 2 : OGRSpatialReference oGCS;
599 2 : const double dfSemiMajor = CPLAtof(papszTokens[0]);
600 2 : const double dfSemiMinor = CPLAtof(papszTokens[1]);
601 : const double dfInvFlattening =
602 2 : OSRCalcInvFlattening(dfSemiMajor, dfSemiMinor);
603 2 : oGCS.SetGeogCS("unknown", "unknown", "unknown",
604 : dfSemiMajor, dfInvFlattening);
605 2 : CopyGeogCSFrom(&oGCS);
606 2 : bFoundParameters = true;
607 : }
608 3 : CSLDestroy(papszTokens);
609 3 : break;
610 : }
611 : }
612 3 : if (!bFoundParameters)
613 : {
614 1 : if (bHasRadiusOfSphereOfReference)
615 : {
616 2 : OGRSpatialReference oGCS;
617 1 : oGCS.SetGeogCS("unknown", "unknown", "unknown",
618 : dfRadiusOfSphereOfReference, 0);
619 1 : CopyGeogCSFrom(&oGCS);
620 : }
621 : else
622 : {
623 : // If unknown, default to WGS84 so there is something there.
624 0 : SetWellKnownGeogCS("WGS84");
625 : }
626 : }
627 : }
628 : }
629 : }
630 :
631 : /* -------------------------------------------------------------------- */
632 : /* Linear units translation */
633 : /* -------------------------------------------------------------------- */
634 37 : if (IsLocal() || IsProjected())
635 : {
636 34 : const double dfOldUnits = GetLinearUnits();
637 68 : const CPLString osValue = OSR_GDS(papszPrj, "Units", "");
638 68 : CPLString osOldAuth;
639 : {
640 34 : const char *pszOldAuth = GetAuthorityCode(nullptr);
641 34 : if (pszOldAuth)
642 4 : osOldAuth = pszOldAuth;
643 : }
644 :
645 34 : if (EQUAL(osValue, ""))
646 0 : SetLinearUnitsAndUpdateParameters(SRS_UL_METER, 1.0);
647 34 : else if (EQUAL(osValue, "FEET"))
648 2 : SetLinearUnitsAndUpdateParameters(SRS_UL_US_FOOT,
649 : CPLAtof(SRS_UL_US_FOOT_CONV));
650 32 : else if (CPLAtof(osValue) != 0.0)
651 1 : SetLinearUnitsAndUpdateParameters("user-defined",
652 1 : 1.0 / CPLAtof(osValue));
653 : else
654 31 : SetLinearUnitsAndUpdateParameters(osValue, 1.0);
655 :
656 : // Reinstall authority if linear units value has not changed (bug #1697)
657 34 : const double dfNewUnits = GetLinearUnits();
658 38 : if (IsProjected() && !osOldAuth.empty() && dfOldUnits != 0.0 &&
659 4 : std::abs(dfNewUnits / dfOldUnits - 1) < 1e-8)
660 : {
661 1 : SetAuthority("PROJCS", "EPSG", atoi(osOldAuth));
662 : }
663 : }
664 :
665 37 : return OGRERR_NONE;
666 : }
667 :
668 : /************************************************************************/
669 : /* FindCodeFromDict() */
670 : /* */
671 : /* Find the code from a dict file. */
672 : /************************************************************************/
673 0 : static int FindCodeFromDict(const char *pszDictFile, const char *CSName,
674 : char *code)
675 : {
676 : /* -------------------------------------------------------------------- */
677 : /* Find and open file. */
678 : /* -------------------------------------------------------------------- */
679 0 : const char *pszFilename = CPLFindFile("gdal", pszDictFile);
680 0 : if (pszFilename == nullptr)
681 0 : return OGRERR_UNSUPPORTED_SRS;
682 :
683 0 : VSILFILE *fp = VSIFOpenL(pszFilename, "rb");
684 0 : if (fp == nullptr)
685 0 : return OGRERR_UNSUPPORTED_SRS;
686 :
687 : /* -------------------------------------------------------------------- */
688 : /* Process lines. */
689 : /* -------------------------------------------------------------------- */
690 0 : OGRErr eErr = OGRERR_UNSUPPORTED_SRS;
691 0 : const char *pszLine = nullptr;
692 :
693 0 : while ((pszLine = CPLReadLineL(fp)) != nullptr)
694 : {
695 0 : if (pszLine[0] == '#')
696 0 : continue;
697 :
698 0 : if (strstr(pszLine, CSName))
699 : {
700 0 : const char *pComma = strchr(pszLine, ',');
701 0 : if (pComma)
702 : {
703 0 : strncpy(code, pszLine, pComma - pszLine);
704 0 : code[pComma - pszLine] = '\0';
705 0 : eErr = OGRERR_NONE;
706 : }
707 0 : break;
708 : }
709 : }
710 :
711 : /* -------------------------------------------------------------------- */
712 : /* Cleanup */
713 : /* -------------------------------------------------------------------- */
714 0 : VSIFCloseL(fp);
715 :
716 0 : return eErr;
717 : }
718 :
719 : /************************************************************************/
720 : /* ImportFromESRIStatePlaneWKT() */
721 : /* */
722 : /* Search a ESRI State Plane WKT and import it. */
723 : /************************************************************************/
724 :
725 3 : OGRErr OGRSpatialReference::ImportFromESRIStatePlaneWKT(int code,
726 : const char *datumName,
727 : const char *unitsName,
728 : int pcsCode,
729 : const char *csName)
730 : {
731 : // If the CS name is known.
732 3 : if (code == 0 && !datumName && !unitsName && pcsCode == 32767 && csName)
733 : {
734 0 : char codeS[10] = {};
735 0 : if (FindCodeFromDict("esri_StatePlane_extra.wkt", csName, codeS) !=
736 : OGRERR_NONE)
737 0 : return OGRERR_FAILURE;
738 0 : return importFromDict("esri_StatePlane_extra.wkt", codeS);
739 : }
740 :
741 3 : int searchCode = -1;
742 3 : if (unitsName == nullptr)
743 0 : unitsName = "";
744 :
745 : // Find state plane prj str by pcs code only.
746 3 : if (code == 0 && !datumName && pcsCode != 32767)
747 : {
748 0 : int unitCode = 1;
749 0 : if (EQUAL(unitsName, "international_feet"))
750 0 : unitCode = 3;
751 0 : else if (strstr(unitsName, "feet") || strstr(unitsName, "foot"))
752 0 : unitCode = 2;
753 :
754 0 : for (int i = 0; statePlanePcsCodeToZoneCode[i] != 0; i += 2)
755 : {
756 0 : if (pcsCode == statePlanePcsCodeToZoneCode[i])
757 : {
758 0 : searchCode = statePlanePcsCodeToZoneCode[i + 1];
759 0 : const int unitIndex = searchCode % 10;
760 0 : if ((unitCode == 1 && !(unitIndex == 0 || unitIndex == 1)) ||
761 0 : (unitCode == 2 &&
762 0 : !(unitIndex == 2 || unitIndex == 3 || unitIndex == 4)) ||
763 0 : (unitCode == 3 && !(unitIndex == 5 || unitIndex == 6)))
764 : {
765 0 : searchCode -= unitIndex;
766 0 : switch (unitIndex)
767 : {
768 0 : case 0:
769 : case 3:
770 : case 5:
771 0 : if (unitCode == 2)
772 0 : searchCode += 3;
773 0 : else if (unitCode == 3)
774 0 : searchCode += 5;
775 0 : break;
776 0 : case 1:
777 : case 2:
778 : case 6:
779 0 : if (unitCode == 1)
780 0 : searchCode += 1;
781 0 : if (unitCode == 2)
782 0 : searchCode += 2;
783 0 : else if (unitCode == 3)
784 0 : searchCode += 6;
785 0 : break;
786 0 : case 4:
787 : // FIXME? The following cond is not possible:
788 : // if( unitCode == 2 )
789 : // searchCode += 4;
790 0 : break;
791 : }
792 : }
793 0 : break;
794 : }
795 0 : }
796 : }
797 : else // Find state plane prj str by all inputs.
798 : {
799 3 : if (code < 0 || code > INT_MAX / 10)
800 0 : return OGRERR_FAILURE;
801 :
802 : // Need to have a special EPSG-ESRI zone code mapping first.
803 360 : for (int i = 0; statePlaneZoneMapping[i] != 0; i += 3)
804 : {
805 357 : if (code == statePlaneZoneMapping[i] &&
806 0 : (statePlaneZoneMapping[i + 1] == -1 ||
807 0 : pcsCode == statePlaneZoneMapping[i + 1]))
808 : {
809 0 : code = statePlaneZoneMapping[i + 2];
810 0 : break;
811 : }
812 : }
813 3 : searchCode = code * 10;
814 3 : if (!datumName)
815 : {
816 0 : CPLError(CE_Failure, CPLE_AppDefined, "datumName is NULL.");
817 0 : return OGRERR_FAILURE;
818 : }
819 3 : if (EQUAL(datumName, "HARN"))
820 : {
821 0 : if (EQUAL(unitsName, "international_feet"))
822 0 : searchCode += 5;
823 0 : else if (strstr(unitsName, "feet") || strstr(unitsName, "foot"))
824 0 : searchCode += 3;
825 : }
826 3 : else if (strstr(datumName, "NAD") && strstr(datumName, "83"))
827 : {
828 3 : if (EQUAL(unitsName, "meters"))
829 0 : searchCode += 1;
830 3 : else if (EQUAL(unitsName, "international_feet"))
831 0 : searchCode += 6;
832 3 : else if (strstr(unitsName, "feet") || strstr(unitsName, "foot"))
833 3 : searchCode += 2;
834 : }
835 0 : else if (strstr(datumName, "NAD") && strstr(datumName, "27") &&
836 0 : !EQUAL(unitsName, "meters"))
837 : {
838 0 : searchCode += 4;
839 : }
840 : else
841 0 : searchCode = -1;
842 : }
843 3 : if (searchCode > 0)
844 : {
845 3 : char codeS[20] = {};
846 3 : snprintf(codeS, sizeof(codeS), "%d", static_cast<int>(searchCode));
847 3 : return importFromDict("esri_StatePlane_extra.wkt", codeS);
848 : }
849 0 : return OGRERR_FAILURE;
850 : }
|