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