1 /******************************************************************************
2 * Copyright (c) 2016, Howard Butler (howard@hobu.co)
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 Hobu, Inc. or Flaxen Geo Consulting nor the
17 * names of its contributors may be used to endorse or promote
18 * products derived from this software without specific prior
19 * 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 #pragma warning(push)
36 #pragma warning(disable: 4251)
37 #include <ogr_api.h>
38 #include <ogr_geometry.h>
39 #pragma warning(pop)
40
41 #include <pdal/Polygon.hpp>
42 #include <pdal/private/gdal/GDALUtils.hpp>
43
44 #include "../filters/private/pnp/GridPnp.hpp"
45
46 namespace pdal
47 {
48
49 struct Polygon::PrivateData
50 {
51 std::vector<GridPnp> m_grids;
52 };
53
54
Polygon()55 Polygon::Polygon()
56 {
57 init();
58 }
59
60
~Polygon()61 Polygon::~Polygon()
62 {}
63
64
Polygon(OGRGeometryH g)65 Polygon::Polygon(OGRGeometryH g) : Geometry(g)
66 {
67 init();
68 }
69
70
Polygon(OGRGeometryH g,const SpatialReference & srs)71 Polygon::Polygon(OGRGeometryH g, const SpatialReference& srs) : Geometry(g, srs)
72 {
73 init();
74 }
75
76
init()77 void Polygon::init()
78 {
79 m_pd.reset(new PrivateData());
80
81 // If the handle was null, we need to create an empty polygon.
82 if (!m_geom)
83 {
84 m_geom.reset(new OGRPolygon());
85 return;
86 }
87
88 OGRwkbGeometryType t = m_geom->getGeometryType();
89
90 if (!(t == wkbPolygon ||
91 t == wkbMultiPolygon ||
92 t == wkbPolygon25D ||
93 t == wkbMultiPolygon25D))
94 {
95 throw pdal::pdal_error("pdal::Polygon() cannot construct geometry "
96 "because OGR geometry is not Polygon or MultiPolygon.");
97 }
98 }
99
100
Polygon(const std::string & wkt_or_json,SpatialReference ref)101 Polygon::Polygon(const std::string& wkt_or_json, SpatialReference ref) :
102 Geometry(wkt_or_json, ref), m_pd(new PrivateData)
103 {}
104
105
Polygon(const BOX2D & box)106 Polygon::Polygon(const BOX2D& box) : m_pd(new PrivateData)
107 {
108 OGRPolygon *poly = new OGRPolygon();
109 m_geom.reset(poly);
110 OGRLinearRing *lr = new OGRLinearRing();
111 lr->addPoint(box.minx, box.miny);
112 lr->addPoint(box.maxx, box.miny);
113 lr->addPoint(box.maxx, box.maxy);
114 lr->addPoint(box.minx, box.maxy);
115 lr->addPoint(box.minx, box.miny);
116 poly->addRingDirectly(lr);
117 }
118
119
Polygon(const BOX3D & box)120 Polygon::Polygon(const BOX3D& box) : m_pd(new PrivateData)
121 {
122 OGRPolygon *poly = new OGRPolygon();
123 m_geom.reset(poly);
124 OGRLinearRing *lr = new OGRLinearRing();
125 lr->addPoint(box.minx, box.miny, box.minz);
126 lr->addPoint(box.minx, box.maxy, box.minz);
127 lr->addPoint(box.maxx, box.maxy, box.maxz);
128 lr->addPoint(box.maxx, box.miny, box.maxz);
129 lr->addPoint(box.minx, box.miny, box.minz);
130 poly->addRingDirectly(lr);
131 }
132
133
Polygon(const Polygon & poly)134 Polygon::Polygon(const Polygon& poly) : Geometry(poly)
135 {
136 init();
137 }
138
139
operator =(const Polygon & src)140 Polygon& Polygon::operator=(const Polygon& src)
141 {
142 ((Geometry *)this)->operator=((const Geometry&)src);
143 m_pd.reset(new PrivateData);
144 return *this;
145 }
146
147
modified()148 void Polygon::modified()
149 {
150 m_pd->m_grids.clear();
151 }
152
153
clear()154 void Polygon::clear()
155 {
156 m_geom.reset(new OGRPolygon());
157 modified();
158 }
159
160
simplify(double distance_tolerance,double area_tolerance,bool preserve_topology)161 void Polygon::simplify(double distance_tolerance, double area_tolerance,
162 bool preserve_topology)
163 {
164 throwNoGeos();
165
166 if (preserve_topology)
167 m_geom.reset(m_geom->SimplifyPreserveTopology(distance_tolerance));
168 else
169 m_geom.reset(m_geom->Simplify(distance_tolerance));
170
171 removeSmallRings(area_tolerance);
172 removeSmallHoles(area_tolerance);
173 modified();
174 }
175
176
removeSmallRings(double tolerance)177 void Polygon::removeSmallRings(double tolerance)
178 {
179 OGRwkbGeometryType t = m_geom->getGeometryType();
180 if (t == wkbPolygon || t == wkbPolygon25D)
181 {
182 if (area() < tolerance)
183 clear();
184 }
185 else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
186 {
187 OGRMultiPolygon *mPoly = static_cast<OGRMultiPolygon *>(m_geom.get());
188 for (int i = mPoly->getNumGeometries() - 1; i >= 0; --i)
189 {
190 OGRPolygon *poly =
191 static_cast<OGRPolygon *>(mPoly->getGeometryRef(i));
192 if (poly->get_Area() < tolerance)
193 mPoly->removeGeometry(i, true);
194 }
195 }
196 }
197
198
removeSmallHoles(double tolerance)199 void Polygon::removeSmallHoles(double tolerance)
200 {
201 auto remove = [tolerance](OGRPolygon *poly)
202 {
203 for (int i = poly->getNumInteriorRings() - 1; i >= 0; --i)
204 {
205 OGRLinearRing *lr = poly->getInteriorRing(i);
206 if (lr->get_Area() < tolerance)
207 OGR_G_RemoveGeometry(gdal::toHandle(poly), i + 1, true);
208 }
209 };
210
211 OGRwkbGeometryType t = m_geom->getGeometryType();
212 if (t == wkbPolygon || t == wkbPolygon25D)
213 remove(static_cast<OGRPolygon *>(m_geom.get()));
214 else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
215 {
216 OGRMultiPolygon *mPoly = static_cast<OGRMultiPolygon *>(m_geom.get());
217 for (int i = mPoly->getNumGeometries() - 1; i >= 0; --i)
218 remove(static_cast<OGRPolygon *>(mPoly->getGeometryRef(i)));
219 }
220 }
221
222
area() const223 double Polygon::area() const
224 {
225 if (!valid())
226 return 0;
227
228 throwNoGeos();
229
230 OGRwkbGeometryType t = m_geom->getGeometryType();
231 // Not until GDAL 2.3
232 /**
233 if (t == wkbPolygon || t == wkbPolygon25D)
234 return m_geom->toPolygon()->get_Area();
235 else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
236 return m_geom->toMultiPolygon()->get_Area();
237 **/
238 if (t == wkbPolygon || t == wkbPolygon25D)
239 {
240 OGRPolygon *p = static_cast<OGRPolygon *>(m_geom.get());
241 return p->get_Area();
242 }
243 else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
244 {
245 OGRMultiPolygon *p = static_cast<OGRMultiPolygon *>(m_geom.get());
246 return p->get_Area();
247 }
248 return 0;
249 }
250
251
covers(const PointRef & ref) const252 bool Polygon::covers(const PointRef& ref) const
253 {
254 throwNoGeos();
255
256 double x = ref.getFieldAs<double>(Dimension::Id::X);
257 double y = ref.getFieldAs<double>(Dimension::Id::Y);
258 double z = ref.getFieldAs<double>(Dimension::Id::Z);
259
260 OGRPoint p(x, y, z);
261 return m_geom->Contains(&p) || m_geom->Touches(&p);
262 }
263
264
overlaps(const Polygon & p) const265 bool Polygon::overlaps(const Polygon& p) const
266 {
267 throwNoGeos();
268
269 return m_geom->Overlaps(p.m_geom.get());
270 }
271
contains(const Polygon & p) const272 bool Polygon::contains(const Polygon& p) const
273 {
274 throwNoGeos();
275
276 return m_geom->Contains(p.m_geom.get());
277 }
278
disjoint(const Polygon & p) const279 bool Polygon::disjoint(const Polygon& p) const
280 {
281 throwNoGeos();
282
283 return m_geom->Disjoint(p.m_geom.get());
284 }
285
intersects(const Polygon & p) const286 bool Polygon::intersects(const Polygon& p) const
287 {
288 throwNoGeos();
289
290 return m_geom->Intersects(p.m_geom.get());
291 }
292
293 /// Determine whether this polygon contains a point.
294 /// \param x Point x coordinate.
295 /// \param y Point y coordinate.
296 /// \return Whether the polygon contains the point or not.
contains(double x,double y) const297 bool Polygon::contains(double x, double y) const
298 {
299 if (m_pd->m_grids.empty())
300 for (const Polygon& p : polygons())
301 m_pd->m_grids.emplace_back(p.exteriorRing(), p.interiorRings());
302 for (auto& g : m_pd->m_grids)
303 if (g.inside(x, y))
304 return true;
305 return false;
306 }
307
308
touches(const Polygon & p) const309 bool Polygon::touches(const Polygon& p) const
310 {
311 throwNoGeos();
312
313 return m_geom->Touches(p.m_geom.get());
314 }
315
within(const Polygon & p) const316 bool Polygon::within(const Polygon& p) const
317 {
318 throwNoGeos();
319
320 return m_geom->Within(p.m_geom.get());
321 }
322
crosses(const Polygon & p) const323 bool Polygon::crosses(const Polygon& p) const
324 {
325 throwNoGeos();
326
327 return m_geom->Crosses(p.m_geom.get());
328 }
329
polygons() const330 std::vector<Polygon> Polygon::polygons() const
331 {
332 std::vector<Polygon> polys;
333
334 OGRwkbGeometryType t = m_geom->getGeometryType();
335
336 if (t == wkbPolygon || t == wkbPolygon25D)
337 polys.emplace_back(*this);
338 else if (t == wkbMultiPolygon || t == wkbMultiPolygon25D)
339 {
340 // Not until GDAL 2.3
341 /**
342 OGRMultiPolygon *mPoly = m_geom->toMultiPolygon();
343 for (auto it = mPoly->begin(); it != mPoly->end(); ++it)
344 {
345 Polygon p;
346 p.m_geom.reset((*it)->clone());
347 polys.push_back(p);
348 }
349 **/
350 OGRMultiPolygon *mPoly = static_cast<OGRMultiPolygon *>(m_geom.get());
351 for (int i = 0; i < mPoly->getNumGeometries(); ++i)
352 {
353 Polygon p;
354 p.m_geom.reset(mPoly->getGeometryRef(i)->clone());
355 polys.push_back(p);
356 }
357 }
358 return polys;
359 }
360
361
exteriorRing() const362 Polygon::Ring Polygon::exteriorRing() const
363 {
364 Ring r;
365
366 OGRwkbGeometryType t = m_geom->getGeometryType();
367 if (t != wkbPolygon && t != wkbPolygon25D)
368 throw pdal_error("Request for exterior ring on non-polygon.");
369
370 // Not until GDAL 2.3
371 /**
372 OGRLinearRing *er = m_geom->toPolygon()->getExteriorRing();
373
374 // For some reason there's no operator -> on an iterator.
375 for (auto it = er->begin(); it != er->end(); ++it)
376 r.push_back({(*it).getX(), (*it).getY()});
377 **/
378 OGRLinearRing *er =
379 static_cast<OGRPolygon *>(m_geom.get())->getExteriorRing();
380 for (int i = 0; i < er->getNumPoints(); ++i)
381 r.push_back({er->getX(i), er->getY(i)});
382
383 return r;
384 }
385
386
interiorRings() const387 std::vector<Polygon::Ring> Polygon::interiorRings() const
388 {
389 std::vector<Ring> rings;
390
391 OGRwkbGeometryType t = m_geom->getGeometryType();
392 if (t != wkbPolygon && t != wkbPolygon25D)
393 throw pdal_error("Request for exterior ring on non-polygon.");
394
395 // OGRPolygon *poly = m_geom->toPolygon();
396 OGRPolygon *poly = static_cast<OGRPolygon *>(m_geom.get());
397 for (int i = 0; i < poly->getNumInteriorRings(); ++i)
398 {
399 OGRLinearRing *er = poly->getInteriorRing(i);
400
401 Ring r;
402 for (int j = 0; j < er->getNumPoints(); ++j)
403 r.push_back({er->getX(j), er->getY(j)});
404 rings.push_back(r);
405 }
406 return rings;
407 }
408
409 } // namespace pdal
410