1 /******************************************************************************
2 * Copyright (c) 2009, Howard Butler
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following
8 * conditions are met:
9 *
10 * * Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * * Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided
15 * with the distribution.
16 * * Neither the name of the Martin Isenburg or Iowa Department
17 * of Natural Resources nor the names of its contributors may be
18 * used to endorse or promote products derived from this software
19 * without specific prior written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
24 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
25 * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
26 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
27 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
28 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
29 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
31 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
32 * OF SUCH DAMAGE.
33 ****************************************************************************/
34
35 #include <memory>
36
37 #include <pdal/Metadata.hpp>
38 #include <pdal/PDALUtils.hpp>
39 #include <pdal/SpatialReference.hpp>
40 #include <pdal/private/SrsTransform.hpp>
41 #include <pdal/util/FileUtils.hpp>
42
43 // gdal
44 # pragma GCC diagnostic push
45 # pragma GCC diagnostic ignored "-Wfloat-equal"
46 #include <ogr_spatialref.h>
47 # pragma GCC diagnostic pop
48
49 #include <cpl_conv.h>
50
51 #include <pdal/util/Utils.hpp>
52
53 namespace
54 {
55
56 struct OGRDeleter
57 {
operator ()__anona8d4ed7f0111::OGRDeleter58 void operator()(OGRSpatialReference* o)
59 {
60 OSRDestroySpatialReference(o);
61 };
62 };
63
64 using OGRScopedSpatialReference =
65 std::unique_ptr<OGRSpatialReference, OGRDeleter>;
66
ogrCreateSrs(std::string s="")67 OGRScopedSpatialReference ogrCreateSrs(std::string s = "")
68 {
69 return OGRScopedSpatialReference(
70 static_cast<OGRSpatialReference*>(
71 OSRNewSpatialReference(s.size() ? s.c_str() : nullptr)));
72 }
73
74 }
75
76 namespace pdal
77 {
78
SpatialReference(const std::string & s)79 SpatialReference::SpatialReference(const std::string& s)
80 {
81 set(s);
82 }
83
84
85 //NOTE that this ctor allows a string constant to be used in places
86 // where a SpatialReference is extpected.
SpatialReference(const char * s)87 SpatialReference::SpatialReference(const char *s)
88 {
89 set(s);
90 }
91
92
empty() const93 bool SpatialReference::empty() const
94 {
95 return m_wkt.empty();
96 }
97
98
valid() const99 bool SpatialReference::valid() const
100 {
101 OGRSpatialReference current(m_wkt.data());
102
103 return OSRValidate(¤t) == OGRERR_NONE;
104 }
105
106
identifyHorizontalEPSG() const107 std::string SpatialReference::identifyHorizontalEPSG() const
108 {
109 OGRScopedSpatialReference srs(ogrCreateSrs(getHorizontal()));
110
111 if (!srs || srs->AutoIdentifyEPSG() != OGRERR_NONE)
112 return "";
113
114 if (const char* c = srs->GetAuthorityCode(nullptr))
115 return std::string(c);
116
117 return "";
118 }
119
120
identifyVerticalEPSG() const121 std::string SpatialReference::identifyVerticalEPSG() const
122 {
123 OGRScopedSpatialReference srs(ogrCreateSrs(getVertical()));
124
125 if (!srs || srs->AutoIdentifyEPSG() != OGRERR_NONE)
126 return "";
127
128 if (const char* c = srs->GetAuthorityCode(nullptr))
129 return std::string(c);
130
131 return "";
132 }
133
134
getWKT() const135 std::string SpatialReference::getWKT() const
136 {
137 return m_wkt;
138 }
139
140
parse(const std::string & s,std::string::size_type & pos)141 void SpatialReference::parse(const std::string& s, std::string::size_type& pos)
142 {
143 set(s.substr(pos));
144 }
145
146
set(std::string v)147 void SpatialReference::set(std::string v)
148 {
149 m_horizontalWkt.clear();
150 if (v.empty())
151 {
152 m_wkt.clear();
153 return;
154 }
155
156 if (isWKT(v))
157 {
158 m_wkt = v;
159 return;
160 }
161
162 OGRSpatialReference srs(NULL);
163
164 CPLErrorReset();
165 const char* input = v.c_str();
166 OGRErr err = srs.SetFromUserInput(const_cast<char *>(input));
167 if (err != OGRERR_NONE)
168 {
169 std::ostringstream oss;
170 std::string msg = CPLGetLastErrorMsg();
171 if (msg.empty())
172 msg = "(unknown reason)";
173 oss << "Could not import coordinate system '" << input << "': " <<
174 msg << ".";
175 throw pdal_error(oss.str());
176 }
177
178 char *poWKT = 0;
179 srs.exportToWkt(&poWKT);
180 m_wkt = poWKT;
181 CPLFree(poWKT);
182 }
183
184
getProj4() const185 std::string SpatialReference::getProj4() const
186 {
187 std::string tmp;
188
189 const char* poWKT = m_wkt.c_str();
190
191 OGRSpatialReference srs(NULL);
192 if (OGRERR_NONE == srs.importFromWkt(const_cast<char **>(&poWKT)))
193 {
194 char* proj4 = nullptr;
195 srs.exportToProj4(&proj4);
196 tmp = proj4;
197 CPLFree(proj4);
198 Utils::trim(tmp);
199 }
200
201 return tmp;
202 }
203
204
getVertical() const205 std::string SpatialReference::getVertical() const
206 {
207 std::string tmp;
208
209 OGRScopedSpatialReference poSRS = ogrCreateSrs(m_wkt);
210
211 // Above can fail if m_wkt is bad.
212 if (!poSRS)
213 return tmp;
214
215 char *pszWKT = NULL;
216
217 OGR_SRSNode* node = poSRS->GetAttrNode("VERT_CS");
218 if (node && poSRS)
219 {
220 node->exportToWkt(&pszWKT);
221 tmp = pszWKT;
222 CPLFree(pszWKT);
223 }
224
225 return tmp;
226 }
227
228
getVerticalUnits() const229 std::string SpatialReference::getVerticalUnits() const
230 {
231 std::string tmp;
232
233 OGRScopedSpatialReference poSRS = ogrCreateSrs(m_wkt);
234 if (poSRS)
235 {
236 OGR_SRSNode* node = poSRS->GetAttrNode("VERT_CS");
237 if (node)
238 {
239 char* units(nullptr);
240
241 // 'units' remains internal to the OGRSpatialReference
242 // and should not be freed, or modified. It may be invalidated
243 // on the next OGRSpatialReference call.
244 (void)poSRS->GetLinearUnits(&units);
245 tmp = units;
246 Utils::trim(tmp);
247 }
248 }
249
250 return tmp;
251 }
252
253
getHorizontal() const254 std::string SpatialReference::getHorizontal() const
255 {
256 if (m_horizontalWkt.empty())
257 {
258 OGRScopedSpatialReference poSRS = ogrCreateSrs(m_wkt);
259
260 if (poSRS)
261 {
262 char *pszWKT(nullptr);
263 poSRS->StripVertical();
264 poSRS->exportToWkt(&pszWKT);
265 m_horizontalWkt = pszWKT;
266 CPLFree(pszWKT);
267 }
268 }
269 return m_horizontalWkt;
270 }
271
272
getHorizontalUnits() const273 std::string SpatialReference::getHorizontalUnits() const
274 {
275 OGRScopedSpatialReference poSRS = ogrCreateSrs(m_wkt);
276
277 if (!poSRS)
278 return std::string();
279
280 char* units(nullptr);
281
282 // The returned value remains internal to the OGRSpatialReference
283 // and should not be freed, or modified. It may be invalidated on
284 // the next OGRSpatialReference call.
285 double u = poSRS->GetLinearUnits(&units);
286 std::string tmp(units);
287 Utils::trim(tmp);
288 return tmp;
289 }
290
291
equals(const SpatialReference & input) const292 bool SpatialReference::equals(const SpatialReference& input) const
293 {
294 if (getWKT() == input.getWKT())
295 return true;
296
297 OGRScopedSpatialReference current = ogrCreateSrs(getWKT());
298 OGRScopedSpatialReference other = ogrCreateSrs(input.getWKT());
299
300 if (!current || !other)
301 return false;
302
303 int output = OSRIsSame(current.get(), other.get());
304
305 return (output == 1);
306 }
307
308
operator ==(const SpatialReference & input) const309 bool SpatialReference::operator==(const SpatialReference& input) const
310 {
311 return this->equals(input);
312 }
313
314
operator !=(const SpatialReference & input) const315 bool SpatialReference::operator!=(const SpatialReference& input) const
316 {
317 return !(this->equals(input));
318 }
319
320
getName() const321 const std::string& SpatialReference::getName() const
322 {
323 static std::string name("pdal.spatialreference");
324 return name;
325 }
326
327
isGeographic() const328 bool SpatialReference::isGeographic() const
329 {
330 OGRScopedSpatialReference current = ogrCreateSrs(m_wkt);
331 if (!current)
332 return false;
333
334 bool output = OSRIsGeographic(current.get());
335 return output;
336 }
337
338
isGeocentric() const339 bool SpatialReference::isGeocentric() const
340 {
341 OGRScopedSpatialReference current = ogrCreateSrs(m_wkt);
342 if (!current)
343 return false;
344
345 bool output = OSRIsGeocentric(current.get());
346 return output;
347 }
348
349
isProjected() const350 bool SpatialReference::isProjected() const
351 {
352 OGRScopedSpatialReference current = ogrCreateSrs(m_wkt);
353 if (!current)
354 return false;
355
356 bool output = OSRIsProjected(current.get());
357 return output;
358 }
359
getAxisOrdering() const360 std::vector<int> SpatialReference::getAxisOrdering() const
361 {
362 std::vector<int> output;
363 OGRScopedSpatialReference current = ogrCreateSrs(m_wkt);
364 output = current.get()->GetDataAxisToSRSAxisMapping();
365 return output;
366 }
367
368
369
calculateZone(double lon,double lat)370 int SpatialReference::calculateZone(double lon, double lat)
371 {
372 int zone = 0;
373 lon = Utils::normalizeLongitude(lon);
374
375 // Special Norway processing.
376 if (lat >= 56.0 && lat < 64.0 && lon >= 3.0 && lon < 12.0 )
377 zone = 32;
378 // Special Svalbard processing.
379 else if (lat >= 72.0 && lat < 84.0)
380 {
381 if (lon >= 0.0 && lon < 9.0)
382 zone = 31;
383 else if (lon >= 9.0 && lon < 21.0)
384 zone = 33;
385 else if (lon >= 21.0 && lon < 33.0)
386 zone = 35;
387 else if (lon >= 33.0 && lon < 42.0)
388 zone = 37;
389 }
390 // Everywhere else.
391 else
392 {
393 zone = (int) floor((lon + 180.0) / 6) + 1;
394 if (lat < 0)
395 zone = -zone;
396 }
397
398 return zone;
399 }
400
401
402 /**
403 Create a spatial reference that represents a specific UTM zone.
404
405 \param zone Zone number. Must be non-zero and <= 60 and >= -60
406 \return A SpatialReference that represents the specified zone, or
407 an invalid SpatialReference on error.
408 */
wgs84FromZone(int zone)409 SpatialReference SpatialReference::wgs84FromZone(int zone)
410 {
411 uint32_t abszone(std::abs(zone));
412
413 if (abszone == 0 || abszone > 60)
414 return SpatialReference();
415
416 std::string code;
417 if (zone > 0)
418 code = "EPSG:326";
419 else
420 code = "EPSG:327";
421
422 code += ((abszone < 10) ? "0" : "") + Utils::toString(abszone);
423 return SpatialReference(code);
424 }
425
426
isWKT(const std::string & wkt)427 bool SpatialReference::isWKT(const std::string& wkt)
428 {
429 // List comes from GDAL. WKT includes FITTED_CS, but this isn't
430 // included in GDAL list. Not sure why.
431 StringList leaders { "PROJCS", "GEOGCS", "COMPD_CS", "GEOCCS",
432 "VERT_CS", "LOCAL_CS",
433 // New specification names.
434 "GEODCRS", "GEODETICCRS",
435 "PROJCRS", "PROJECTEDCRS",
436 "VERTCRS", "VERITCALCRS",
437 "ENGCRS", "ENGINEERINGCRS",
438 "IMAGECRS",
439 "PARAMETRICCRS",
440 "TIMECRS",
441 "COMPOUNDCRS"
442 };
443
444 for (const std::string& s : leaders)
445 if (wkt.compare(0, s.size(), s) == 0)
446 return true;
447 return false;
448 }
449
450
prettyWkt(const std::string & wkt)451 std::string SpatialReference::prettyWkt(const std::string& wkt)
452 {
453 std::string outWkt;
454
455 OGRScopedSpatialReference srs = ogrCreateSrs(wkt);
456 if (!srs)
457 return outWkt;
458
459 char *buf = nullptr;
460 srs->exportToPrettyWkt(&buf, FALSE);
461
462 outWkt = buf;
463 CPLFree(buf);
464 return outWkt;
465 }
466
getWKT1() const467 std::string SpatialReference::getWKT1() const
468 {
469 std::string wkt = getWKT();
470 if (wkt.empty())
471 return wkt;
472
473 OGRScopedSpatialReference srs = ogrCreateSrs(wkt);
474 std::string wkt1;
475 if (srs)
476 {
477 char *buf = nullptr;
478 const char* apszOptions[] =
479 { "FORMAT=WKT1_GDAL", "ALLOW_ELLIPSOIDAL_HEIGHT_AS_VERTICAL_CRS=YES", nullptr };
480
481 srs->exportToWkt(&buf, apszOptions);
482 if (buf)
483 {
484 wkt1 = buf;
485 CPLFree(buf);
486 }
487 }
488 if (wkt1.empty())
489 throw pdal_error("Couldn't convert spatial reference to WKT version 1.");
490 return wkt1;
491 }
492
493
getUTMZone() const494 int SpatialReference::getUTMZone() const
495 {
496 OGRScopedSpatialReference current = ogrCreateSrs(m_wkt);
497 if (!current)
498 throw pdal_error("Could not fetch current SRS");
499
500 int north(0);
501 int zone = OSRGetUTMZone(current.get(), &north);
502 return (north ? 1 : -1) * zone;
503 }
504
505
computeUTMZone(const BOX3D & cbox) const506 int SpatialReference::computeUTMZone(const BOX3D& cbox) const
507 {
508 SrsTransform transform(*this, SpatialReference("EPSG:4326"));
509
510 // We made the argument constant so copy so that we can modify.
511 BOX3D box(cbox);
512
513 transform.transform(box.minx, box.miny, box.minz);
514 transform.transform(box.maxx, box.maxy, box.maxz);
515
516 int minZone = calculateZone(box.minx, box.miny);
517 int maxZone = calculateZone(box.maxx, box.maxy);
518
519 if (minZone != maxZone)
520 {
521 std::ostringstream msg;
522 msg << "computeUTMZone failed due to zone crossing. Zones "
523 "are " << minZone << " and " << maxZone << ".";
524 throw pdal_error(msg.str());
525 }
526 return minZone;
527 }
528
529
toMetadata() const530 MetadataNode SpatialReference::toMetadata() const
531 {
532 MetadataNode root("srs");
533 root.add("horizontal", getHorizontal());
534 root.add("vertical", getVertical());
535 root.add("isgeographic", isGeographic());
536 root.add("isgeocentric", isGeocentric());
537 root.add("proj4", getProj4());
538 root.add("prettywkt", prettyWkt(getHorizontal()));
539 root.add("wkt", getHorizontal());
540 root.add("compoundwkt", getWKT());
541 root.add("prettycompoundwkt", prettyWkt(m_wkt));
542
543 MetadataNode units = root.add("units");
544 units.add("vertical", getVerticalUnits());
545 units.add("horizontal", getHorizontalUnits());
546
547 return root;
548 }
549
550
dump() const551 void SpatialReference::dump() const
552 {
553 std::cout << *this;
554 }
555
556
operator <<(std::ostream & ostr,const SpatialReference & srs)557 std::ostream& operator<<(std::ostream& ostr, const SpatialReference& srs)
558 {
559 ostr << SpatialReference::prettyWkt(srs.m_wkt);
560 return ostr;
561 }
562
563
operator >>(std::istream & istr,SpatialReference & srs)564 std::istream& operator>>(std::istream& istr, SpatialReference& srs)
565 {
566 std::ostringstream oss;
567 oss << istr.rdbuf();
568 srs.set(oss.str());
569
570 return istr;
571 }
572
573 } // namespace pdal
574