Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: Generate an OGRSpatialReference object based on an EPSG
5 : * PROJCS, or GEOGCS code.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2000, Frank Warmerdam
10 : * Copyright (c) 2008-2013, 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 :
36 : #include <cmath>
37 : #include <cstddef>
38 : #include <cstdio>
39 : #include <cstring>
40 : #include <cstdlib>
41 : #include <map>
42 : #include <memory>
43 : #include <string>
44 : #include <vector>
45 : #include <limits>
46 :
47 : #include "cpl_conv.h"
48 : #include "cpl_csv.h"
49 : #include "cpl_error.h"
50 : #include "cpl_multiproc.h"
51 : #include "cpl_string.h"
52 : #include "ogr_core.h"
53 : #include "ogr_p.h"
54 : #include "ogr_proj_p.h"
55 : #include "ogr_spatialref.h"
56 :
57 : #include "proj.h"
58 :
59 : extern void OGRsnPrintDouble(char *pszStrBuf, size_t size, double dfValue);
60 :
61 : /************************************************************************/
62 : /* OSRGetEllipsoidInfo() */
63 : /************************************************************************/
64 :
65 : /**
66 : * Fetch info about an ellipsoid.
67 : *
68 : * This helper function will return ellipsoid parameters corresponding to EPSG
69 : * code provided. Axes are always returned in meters. Semi major computed
70 : * based on inverse flattening where that is provided.
71 : *
72 : * @param nCode EPSG code of the requested ellipsoid
73 : *
74 : * @param ppszName pointer to string where ellipsoid name will be returned. It
75 : * is caller responsibility to free this string after using with CPLFree().
76 : *
77 : * @param pdfSemiMajor pointer to variable where semi major axis will be
78 : * returned.
79 : *
80 : * @param pdfInvFlattening pointer to variable where inverse flattening will
81 : * be returned.
82 : *
83 : * @return OGRERR_NONE on success or an error code in case of failure.
84 : **/
85 :
86 245 : OGRErr OSRGetEllipsoidInfo(int nCode, char **ppszName, double *pdfSemiMajor,
87 : double *pdfInvFlattening)
88 :
89 : {
90 490 : CPLString osCode;
91 245 : osCode.Printf("%d", nCode);
92 245 : auto ellipsoid = proj_create_from_database(
93 : OSRGetProjTLSContext(), "EPSG", osCode.c_str(), PJ_CATEGORY_ELLIPSOID,
94 : false, nullptr);
95 245 : if (!ellipsoid)
96 : {
97 0 : return OGRERR_UNSUPPORTED_SRS;
98 : }
99 :
100 245 : if (ppszName)
101 : {
102 62 : *ppszName = CPLStrdup(proj_get_name(ellipsoid));
103 : }
104 245 : proj_ellipsoid_get_parameters(OSRGetProjTLSContext(), ellipsoid,
105 : pdfSemiMajor, nullptr, nullptr,
106 : pdfInvFlattening);
107 245 : proj_destroy(ellipsoid);
108 :
109 245 : return OGRERR_NONE;
110 : }
111 :
112 : /************************************************************************/
113 : /* SetStatePlane() */
114 : /************************************************************************/
115 :
116 : /**
117 : * \brief Set State Plane projection definition.
118 : *
119 : * This will attempt to generate a complete definition of a state plane
120 : * zone based on generating the entire SRS from the EPSG tables. If the
121 : * EPSG tables are unavailable, it will produce a stubbed LOCAL_CS definition
122 : * and return OGRERR_FAILURE.
123 : *
124 : * This method is the same as the C function OSRSetStatePlaneWithUnits().
125 : *
126 : * @param nZone State plane zone number, in the USGS numbering scheme (as
127 : * distinct from the Arc/Info and Erdas numbering scheme.
128 : *
129 : * @param bNAD83 TRUE if the NAD83 zone definition should be used or FALSE
130 : * if the NAD27 zone definition should be used.
131 : *
132 : * @param pszOverrideUnitName Linear unit name to apply overriding the
133 : * legal definition for this zone.
134 : *
135 : * @param dfOverrideUnit Linear unit conversion factor to apply overriding
136 : * the legal definition for this zone.
137 : *
138 : * @return OGRERR_NONE on success, or OGRERR_FAILURE on failure, mostly likely
139 : * due to the EPSG tables not being accessible.
140 : */
141 :
142 10 : OGRErr OGRSpatialReference::SetStatePlane(int nZone, int bNAD83,
143 : const char *pszOverrideUnitName,
144 : double dfOverrideUnit)
145 :
146 : {
147 :
148 : /* -------------------------------------------------------------------- */
149 : /* Get the index id from stateplane.csv. */
150 : /* -------------------------------------------------------------------- */
151 :
152 10 : if (!bNAD83 && nZone > INT_MAX - 10000)
153 0 : return OGRERR_FAILURE;
154 :
155 10 : const int nAdjustedId = bNAD83 ? nZone : nZone + 10000;
156 :
157 : /* -------------------------------------------------------------------- */
158 : /* Turn this into a PCS code. We assume there will only be one */
159 : /* PCS corresponding to each Proj_ code since the proj code */
160 : /* already effectively indicates NAD27 or NAD83. */
161 : /* -------------------------------------------------------------------- */
162 10 : char szID[32] = {};
163 10 : snprintf(szID, sizeof(szID), "%d", nAdjustedId);
164 10 : const int nPCSCode = atoi(CSVGetField(CSVFilename("stateplane.csv"), "ID",
165 : szID, CC_Integer, "EPSG_PCS_CODE"));
166 10 : if (nPCSCode < 1)
167 : {
168 : static bool bFailureReported = false;
169 :
170 0 : if (!bFailureReported)
171 : {
172 0 : bFailureReported = true;
173 0 : CPLError(CE_Warning, CPLE_OpenFailed,
174 : "Unable to find state plane zone in stateplane.csv, "
175 : "likely because the GDAL data files cannot be found. "
176 : "Using incomplete definition of state plane zone.");
177 : }
178 :
179 0 : Clear();
180 0 : if (bNAD83)
181 : {
182 0 : char szName[128] = {};
183 0 : snprintf(szName, sizeof(szName), "State Plane Zone %d / NAD83",
184 : nZone);
185 0 : SetLocalCS(szName);
186 0 : SetLinearUnits(SRS_UL_METER, 1.0);
187 : }
188 : else
189 : {
190 0 : char szName[128] = {};
191 0 : snprintf(szName, sizeof(szName), "State Plane Zone %d / NAD27",
192 : nZone);
193 0 : SetLocalCS(szName);
194 0 : SetLinearUnits(SRS_UL_US_FOOT, CPLAtof(SRS_UL_US_FOOT_CONV));
195 : }
196 :
197 0 : return OGRERR_FAILURE;
198 : }
199 :
200 : /* -------------------------------------------------------------------- */
201 : /* Define based on a full EPSG definition of the zone. */
202 : /* -------------------------------------------------------------------- */
203 10 : OGRErr eErr = importFromEPSG(nPCSCode);
204 :
205 10 : if (eErr != OGRERR_NONE)
206 0 : return eErr;
207 :
208 : /* -------------------------------------------------------------------- */
209 : /* Apply units override if required. */
210 : /* */
211 : /* We will need to adjust the linear projection parameter to */
212 : /* match the provided units, and clear the authority code. */
213 : /* -------------------------------------------------------------------- */
214 14 : if (pszOverrideUnitName != nullptr && dfOverrideUnit != 0.0 &&
215 4 : fabs(dfOverrideUnit - GetLinearUnits()) > 0.0000000001)
216 : {
217 2 : const double dfFalseEasting = GetNormProjParm(SRS_PP_FALSE_EASTING);
218 2 : const double dfFalseNorthing = GetNormProjParm(SRS_PP_FALSE_NORTHING);
219 :
220 2 : SetLinearUnits(pszOverrideUnitName, dfOverrideUnit);
221 :
222 2 : SetNormProjParm(SRS_PP_FALSE_EASTING, dfFalseEasting);
223 2 : SetNormProjParm(SRS_PP_FALSE_NORTHING, dfFalseNorthing);
224 :
225 2 : OGR_SRSNode *const poPROJCS = GetAttrNode("PROJCS");
226 2 : if (poPROJCS != nullptr && poPROJCS->FindChild("AUTHORITY") != -1)
227 : {
228 0 : poPROJCS->DestroyChild(poPROJCS->FindChild("AUTHORITY"));
229 : }
230 : }
231 :
232 10 : return OGRERR_NONE;
233 : }
234 :
235 : /************************************************************************/
236 : /* OSRSetStatePlane() */
237 : /************************************************************************/
238 :
239 : /**
240 : * \brief Set State Plane projection definition.
241 : *
242 : * This function is the same as OGRSpatialReference::SetStatePlane().
243 : */
244 :
245 1 : OGRErr OSRSetStatePlane(OGRSpatialReferenceH hSRS, int nZone, int bNAD83)
246 :
247 : {
248 1 : VALIDATE_POINTER1(hSRS, "OSRSetStatePlane", OGRERR_FAILURE);
249 :
250 1 : return reinterpret_cast<OGRSpatialReference *>(hSRS)->SetStatePlane(nZone,
251 1 : bNAD83);
252 : }
253 :
254 : /************************************************************************/
255 : /* OSRSetStatePlaneWithUnits() */
256 : /************************************************************************/
257 :
258 : /**
259 : * \brief Set State Plane projection definition.
260 : *
261 : * This function is the same as OGRSpatialReference::SetStatePlane().
262 : */
263 :
264 3 : OGRErr OSRSetStatePlaneWithUnits(OGRSpatialReferenceH hSRS, int nZone,
265 : int bNAD83, const char *pszOverrideUnitName,
266 : double dfOverrideUnit)
267 :
268 : {
269 3 : VALIDATE_POINTER1(hSRS, "OSRSetStatePlaneWithUnits", OGRERR_FAILURE);
270 :
271 3 : return reinterpret_cast<OGRSpatialReference *>(hSRS)->SetStatePlane(
272 3 : nZone, bNAD83, pszOverrideUnitName, dfOverrideUnit);
273 : }
274 :
275 : /************************************************************************/
276 : /* AutoIdentifyEPSG() */
277 : /************************************************************************/
278 :
279 : /**
280 : * \brief Set EPSG authority info if possible.
281 : *
282 : * This method inspects a WKT definition, and adds EPSG authority nodes
283 : * where an aspect of the coordinate system can be easily and safely
284 : * corresponded with an EPSG identifier. In practice, this method will
285 : * evolve over time. In theory it can add authority nodes for any object
286 : * (i.e. spheroid, datum, GEOGCS, units, and PROJCS) that could have an
287 : * authority node. Mostly this is useful to inserting appropriate
288 : * PROJCS codes for common formulations (like UTM n WGS84).
289 : *
290 : * If it success the OGRSpatialReference is updated in place, and the
291 : * method return OGRERR_NONE. If the method fails to identify the
292 : * general coordinate system OGRERR_UNSUPPORTED_SRS is returned but no
293 : * error message is posted via CPLError().
294 : *
295 : * This method is the same as the C function OSRAutoIdentifyEPSG().
296 : *
297 : * Since GDAL 2.3, the FindMatches() method can also be used for improved
298 : * matching by researching the EPSG catalog.
299 : *
300 : * @return OGRERR_NONE or OGRERR_UNSUPPORTED_SRS.
301 : *
302 : * @see OGRSpatialReference::FindBestMatch()
303 : */
304 :
305 312 : OGRErr OGRSpatialReference::AutoIdentifyEPSG()
306 :
307 : {
308 : /* -------------------------------------------------------------------- */
309 : /* Do we have a GEOGCS node, but no authority? If so, try */
310 : /* guessing it. */
311 : /* -------------------------------------------------------------------- */
312 417 : if ((IsProjected() || IsGeographic()) &&
313 571 : GetAuthorityCode("GEOGCS") == nullptr &&
314 154 : GetAttrValue("PROJCRS|BASEGEOGCRS|ID") == nullptr)
315 : {
316 152 : const int nGCS = GetEPSGGeogCS();
317 152 : if (nGCS != -1)
318 46 : SetAuthority("GEOGCS", "EPSG", nGCS);
319 : }
320 :
321 312 : if (IsProjected() && GetAuthorityCode("PROJCS") == nullptr)
322 : {
323 186 : const char *pszProjection = GetAttrValue("PROJECTION");
324 :
325 : /* --------------------------------------------------------------------
326 : */
327 : /* Is this a UTM coordinate system with a common GEOGCS? */
328 : /* --------------------------------------------------------------------
329 : */
330 186 : int nZone = 0;
331 186 : int bNorth = FALSE;
332 186 : if ((nZone = GetUTMZone(&bNorth)) != 0)
333 : {
334 106 : const char *pszAuthName = GetAuthorityName("PROJCS|GEOGCS");
335 106 : const char *pszAuthCode = GetAuthorityCode("PROJCS|GEOGCS");
336 :
337 106 : if (pszAuthName == nullptr || pszAuthCode == nullptr)
338 : {
339 : // Don't exactly recognise datum.
340 : }
341 91 : else if (EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4326)
342 : {
343 : // WGS84
344 0 : if (bNorth)
345 0 : SetAuthority("PROJCS", "EPSG", 32600 + nZone);
346 : else
347 0 : SetAuthority("PROJCS", "EPSG", 32700 + nZone);
348 : }
349 91 : else if (EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4267 &&
350 88 : nZone >= 3 && nZone <= 22 && bNorth)
351 : {
352 88 : SetAuthority("PROJCS", "EPSG", 26700 + nZone); // NAD27
353 : }
354 3 : else if (EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4269 &&
355 1 : nZone >= 3 && nZone <= 23 && bNorth)
356 : {
357 1 : SetAuthority("PROJCS", "EPSG", 26900 + nZone); // NAD83
358 : }
359 2 : else if (EQUAL(pszAuthName, "EPSG") && atoi(pszAuthCode) == 4322)
360 : { // WGS72
361 0 : if (bNorth)
362 0 : SetAuthority("PROJCS", "EPSG", 32200 + nZone);
363 : else
364 0 : SetAuthority("PROJCS", "EPSG", 32300 + nZone);
365 : }
366 : }
367 :
368 : /* --------------------------------------------------------------------
369 : */
370 : /* Is this a Polar Stereographic system on WGS 84 ? */
371 : /* --------------------------------------------------------------------
372 : */
373 80 : else if (pszProjection != nullptr &&
374 80 : EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
375 : {
376 4 : const char *pszAuthName = GetAuthorityName("PROJCS|GEOGCS");
377 4 : const char *pszAuthCode = GetAuthorityCode("PROJCS|GEOGCS");
378 : const double dfLatOrigin =
379 4 : GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
380 :
381 4 : if (pszAuthName != nullptr && EQUAL(pszAuthName, "EPSG") &&
382 4 : pszAuthCode != nullptr && atoi(pszAuthCode) == 4326 &&
383 4 : fabs(fabs(dfLatOrigin) - 71.0) < 1e-15 &&
384 3 : fabs(GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)) < 1e-15 &&
385 3 : fabs(GetProjParm(SRS_PP_SCALE_FACTOR, 1.0) - 1.0) < 1e-15 &&
386 3 : fabs(GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0)) < 1e-15 &&
387 11 : fabs(GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0)) < 1e-15 &&
388 3 : fabs(GetLinearUnits() - 1.0) < 1e-15)
389 : {
390 3 : if (dfLatOrigin > 0)
391 : // Arctic Polar Stereographic
392 1 : SetAuthority("PROJCS", "EPSG", 3995);
393 : else
394 : // Antarctic Polar Stereographic
395 2 : SetAuthority("PROJCS", "EPSG", 3031);
396 : }
397 : }
398 : }
399 :
400 : /* -------------------------------------------------------------------- */
401 : /* Return. */
402 : /* -------------------------------------------------------------------- */
403 312 : if (IsProjected() && GetAuthorityCode("PROJCS") != nullptr)
404 113 : return OGRERR_NONE;
405 :
406 199 : if (IsGeographic() && GetAuthorityCode("GEOGCS") != nullptr)
407 30 : return OGRERR_NONE;
408 :
409 169 : return OGRERR_UNSUPPORTED_SRS;
410 : }
411 :
412 : /************************************************************************/
413 : /* OSRAutoIdentifyEPSG() */
414 : /************************************************************************/
415 :
416 : /**
417 : * \brief Set EPSG authority info if possible.
418 : *
419 : * This function is the same as OGRSpatialReference::AutoIdentifyEPSG().
420 : *
421 : * Since GDAL 2.3, the OSRFindMatches() function can also be used for improved
422 : * matching by researching the EPSG catalog.
423 : *
424 : */
425 :
426 5 : OGRErr OSRAutoIdentifyEPSG(OGRSpatialReferenceH hSRS)
427 :
428 : {
429 5 : VALIDATE_POINTER1(hSRS, "OSRAutoIdentifyEPSG", OGRERR_FAILURE);
430 :
431 5 : return reinterpret_cast<OGRSpatialReference *>(hSRS)->AutoIdentifyEPSG();
432 : }
|
