OpenCPN Partial API docs
Loading...
Searching...
No Matches
shapefile_basemap.cpp
1/******************************************************************************
2 *
3 * Project: OpenCPN
4 * Purpose: Shapefile basemap
5 *
6 ***************************************************************************
7 * Copyright (C) 2012-2023 by David S. Register *
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 * This program is distributed in the hope that it will be useful, *
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
17 * GNU General Public License for more details. *
18 * *
19 * You should have received a copy of the GNU General Public License *
20 * along with this program; if not, write to the *
21 * Free Software Foundation, Inc., *
22 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. *
23 ***************************************************************************
24 *
25 *
26 */
27
28#include <list>
29
30// Include OCPNPlatform.h before shapefile_basemap.h to prevent obscure syntax
31// error when compiling with VS2022
32#include <list>
33
34#include "OCPNPlatform.h"
35#include "shapefile_basemap.h"
36#include "chartbase.h"
37#include "glChartCanvas.h"
38
39#include "model/logger.h"
40
41#ifdef ocpnUSE_GL
42#include "shaders.h"
43#endif
44
45#ifdef __WXMSW__
46#define __CALL_CONVENTION //__stdcall
47#else
48#define __CALL_CONVENTION
49#endif
50
51extern OCPNPlatform *g_Platform;
52extern wxString gWorldShapefileLocation;
53
54ShapeBaseChartSet gShapeBasemap;
55
56#ifdef ocpnUSE_GL
57
58typedef union {
59 GLdouble data[6];
60 struct sGLvertex {
61 GLdouble x;
62 GLdouble y;
63 GLdouble z;
64 GLdouble r;
65 GLdouble g;
66 GLdouble b;
67 } info;
68} GLvertexshp;
69
70static std::list<float_2Dpt> g_pvshp;
71static std::list<GLvertexshp *> g_vertexesshp;
72static int g_typeshp, g_posshp;
73static float_2Dpt g_p1shp, g_p2shp;
74
75void __CALL_CONVENTION shpscombineCallback(GLdouble coords[3],
76 GLdouble *vertex_data[4],
77 GLfloat weight[4],
78 GLdouble **dataOut) {
79 GLvertexshp *vertex;
80
81 vertex = new GLvertexshp();
82 g_vertexesshp.push_back(vertex);
83
84 vertex->info.x = coords[0];
85 vertex->info.y = coords[1];
86
87 *dataOut = vertex->data;
88}
89
90void __CALL_CONVENTION shpserrorCallback(GLenum errorCode) {
91 const GLubyte *estring;
92 estring = gluErrorString(errorCode);
93 // wxLogMessage( "OpenGL Tessellation Error: %s", estring );
94}
95
96void __CALL_CONVENTION shpsbeginCallback(GLenum type) {
97 switch (type) {
98 case GL_TRIANGLES:
99 case GL_TRIANGLE_STRIP:
100 case GL_TRIANGLE_FAN:
101 g_typeshp = type;
102 break;
103 default:
104 printf("tess unhandled begin type: %d\n", type);
105 }
106
107 g_posshp = 0;
108}
109
110void __CALL_CONVENTION shpsendCallback() {}
111
112void __CALL_CONVENTION shpsvertexCallback(GLvoid *arg) {
113 GLvertexshp *vertex;
114 vertex = (GLvertexshp *)arg;
115 float_2Dpt p;
116 p.y = vertex->info.x;
117 p.x = vertex->info.y;
118
119 // convert strips and fans into triangles
120 if (g_typeshp != GL_TRIANGLES) {
121 if (g_posshp > 2) {
122 g_pvshp.push_back(g_p1shp);
123 g_pvshp.push_back(g_p2shp);
124 }
125
126 if (g_typeshp == GL_TRIANGLE_STRIP)
127 g_p1shp = g_p2shp;
128 else if (g_posshp == 0)
129 g_p1shp = p;
130 g_p2shp = p;
131 }
132
133 g_pvshp.push_back(p);
134 g_posshp++;
135}
136#endif
137
138ShapeBaseChartSet::ShapeBaseChartSet() : _loaded(false) {
139 land_color = wxColor(170, 175, 80);
140}
141void ShapeBaseChartSet::SetBasemapLandColor(wxColor color) {
142 land_color = color;
143}
144wxColor ShapeBaseChartSet::GetBasemapLandColor() { return land_color; }
145
146wxPoint2DDouble ShapeBaseChartSet::GetDoublePixFromLL(ViewPort &vp, double lat,
147 double lon) {
148 wxPoint2DDouble p = vp.GetDoublePixFromLL(lat, lon);
149 p.m_x -= vp.rv_rect.x, p.m_y -= vp.rv_rect.y;
150 return p;
151}
152
153ShapeBaseChart &ShapeBaseChartSet::LowestQualityBaseMap() {
154 if (_basemap_map.find(Quality::crude) != _basemap_map.end()) {
155 return _basemap_map.at(Quality::crude);
156 }
157 if (_basemap_map.find(Quality::low) != _basemap_map.end()) {
158 return _basemap_map.at(Quality::low);
159 }
160 if (_basemap_map.find(Quality::medium) != _basemap_map.end()) {
161 return _basemap_map.at(Quality::medium);
162 }
163 if (_basemap_map.find(Quality::high) != _basemap_map.end()) {
164 return _basemap_map.at(Quality::high);
165 }
166 return _basemap_map.at(Quality::full);
167}
168
169ShapeBaseChart &ShapeBaseChartSet::HighestQualityBaseMap() {
170 if (_basemap_map.find(Quality::full) != _basemap_map.end()) {
171 return _basemap_map.at(Quality::full);
172 }
173 if (_basemap_map.find(Quality::high) != _basemap_map.end()) {
174 return _basemap_map.at(Quality::high);
175 }
176 if (_basemap_map.find(Quality::medium) != _basemap_map.end()) {
177 return _basemap_map.at(Quality::medium);
178 }
179 if (_basemap_map.find(Quality::low) != _basemap_map.end()) {
180 return _basemap_map.at(Quality::low);
181 }
182 return _basemap_map.at(Quality::crude);
183}
184
185ShapeBaseChart &ShapeBaseChartSet::SelectBaseMap(const size_t &scale) {
186 if (_basemap_map.find(Quality::full) != _basemap_map.end() &&
187 _basemap_map.at(Quality::full).IsUsable() &&
188 scale <= _basemap_map.at(Quality::full).MinScale()) {
189 return _basemap_map.at(Quality::full);
190 } else if (_basemap_map.find(Quality::high) != _basemap_map.end() &&
191 _basemap_map.at(Quality::high).IsUsable() &&
192 scale <= _basemap_map.at(Quality::high).MinScale()) {
193 return _basemap_map.at(Quality::high);
194 } else if (_basemap_map.find(Quality::medium) != _basemap_map.end() &&
195 _basemap_map.at(Quality::medium).IsUsable() &&
196 scale <= _basemap_map.at(Quality::medium).MinScale()) {
197 return _basemap_map.at(Quality::medium);
198 } else if (_basemap_map.find(Quality::low) != _basemap_map.end() &&
199 _basemap_map.at(Quality::low).IsUsable() &&
200 scale <= _basemap_map.at(Quality::low).MinScale()) {
201 return _basemap_map.at(Quality::low);
202 }
203 return LowestQualityBaseMap();
204}
205
206void ShapeBaseChartSet::Reset() {
207 // Establish the Shapefile basemap location
208 wxString basemap_dir;
209 if (gWorldShapefileLocation.empty()) {
210 basemap_dir = g_Platform->GetSharedDataDir();
211 basemap_dir.Append("basemap_shp");
212 } else {
213 basemap_dir = gWorldShapefileLocation;
214 }
215
216 LoadBasemaps(basemap_dir.ToStdString());
217}
218void ShapeBaseChartSet::LoadBasemaps(const std::string &dir) {
219 _loaded = false;
220 _basemap_map.clear();
221
222 if (fs::exists(ShapeBaseChart::ConstructPath(dir, "crude_10x10"))) {
223 auto c = ShapeBaseChart(ShapeBaseChart::ConstructPath(dir, "crude_10x10"),
224 300000000, land_color);
225 c._dmod = 10;
226 _basemap_map.insert(std::make_pair(Quality::crude, c));
227 }
228
229 if (fs::exists(ShapeBaseChart::ConstructPath(dir, "low"))) {
230 _basemap_map.insert(std::make_pair(
231 Quality::low, ShapeBaseChart(ShapeBaseChart::ConstructPath(dir, "low"),
232 15000000, land_color)));
233 }
234 if (fs::exists(ShapeBaseChart::ConstructPath(dir, "medium"))) {
235 _basemap_map.insert(std::make_pair(
236 Quality::medium,
237 ShapeBaseChart(ShapeBaseChart::ConstructPath(dir, "medium"), 1000000,
238 land_color)));
239 }
240 if (fs::exists(ShapeBaseChart::ConstructPath(dir, "high"))) {
241 _basemap_map.insert(std::make_pair(
242 Quality::high,
243 ShapeBaseChart(ShapeBaseChart::ConstructPath(dir, "high"), 300000,
244 land_color)));
245 }
246 if (fs::exists(ShapeBaseChart::ConstructPath(dir, "full"))) {
247 _basemap_map.insert(std::make_pair(
248 Quality::full,
249 ShapeBaseChart(ShapeBaseChart::ConstructPath(dir, "full"), 10000,
250 land_color)));
251 }
252 _loaded = true;
253 // if(_basemap_map.size())
254 // LowestQualityBaseMap().LoadSHP();
255}
256
258 if (!fs::exists(_filename)) {
259 _is_usable = false;
260 return false;
261 }
262 std::unique_ptr<shp::ShapefileReader> temp_reader(
263 new shp::ShapefileReader(_filename));
264 if (!temp_reader->isOpen()) {
265 MESSAGE_LOG << "Shapefile " << _filename << " is not opened";
266 _is_usable = false;
267 return false;
268 }
269 if (!_loading) {
270 // Check if loading was cancelled
271 return false;
272 }
273 auto bounds = temp_reader->getBounds();
274 _is_usable = temp_reader->getCount() > 1 && bounds.getMaxX() <= 180 &&
275 bounds.getMinX() >= -180 && bounds.getMinY() >= -90 &&
276 bounds.getMaxY() <=
277 90; // TODO - Do we care whether the planet is covered?
278 if (!_loading) {
279 // Check if loading was cancelled
280 _is_usable = false;
281 return false;
282 }
283 _is_usable &= temp_reader->getGeometryType() == shp::GeometryType::Polygon;
284 bool has_x = false;
285 bool has_y = false;
286 for (auto field : temp_reader->getFields()) {
287 if (field.getName() == "x") {
288 has_x = true;
289 } else if (field.getName() == "y") {
290 has_y = true;
291 }
292 }
293 _is_tiled = (has_x && has_y);
294 if (_is_usable && _is_tiled) {
295 size_t feat{0};
296 for (auto const &feature : *temp_reader) {
297 if (!_loading) {
298 // Check if loading was cancelled
299 _is_usable = false;
300 return false;
301 }
302 auto f1 = feature.getAttributes();
303 // Create a LatLonKey using the 'y' (latitude) and 'x' (longitude)
304 // attributes These values represent the top-left corner of the tiles
305 _tiles[LatLonKey(std::any_cast<int>(feature.getAttributes()["y"]),
306 std::any_cast<int>(feature.getAttributes()["x"]))]
307 .push_back(feat);
308 feat++;
309 }
310 }
311 if (_loading) { // Only set reader if loading wasn't cancelled
312 _reader = temp_reader.release();
313 }
314 return _is_usable;
315}
316
317void ShapeBaseChart::DoDrawPolygonFilled(ocpnDC &pnt, ViewPort &vp,
318 const shp::Feature &feature) {
319 double old_x = -9999999.0, old_y = -9999999.0;
320 auto polygon = static_cast<shp::Polygon *>(feature.getGeometry());
321 pnt.SetBrush(_color);
322 for (auto &ring : polygon->getRings()) {
323 wxPoint *poly_pt = new wxPoint[ring.getPoints().size()];
324 size_t cnt{0};
325 auto bbox = vp.GetBBox();
326 for (auto &point : ring.getPoints()) {
327 // if (bbox.ContainsMarge(point.getY(), point.getX(), 0.05)) {
328 wxPoint2DDouble q =
329 ShapeBaseChartSet::GetDoublePixFromLL(vp, point.getY(), point.getX());
330 if (round(q.m_x) != round(old_x) || round(q.m_y) != round(old_y)) {
331 poly_pt[cnt].x = round(q.m_x);
332 poly_pt[cnt].y = round(q.m_y);
333 cnt++;
334 }
335 old_x = q.m_x;
336 old_y = q.m_y;
337 //}
338 }
339 if (cnt > 1) {
340 pnt.DrawPolygonTessellated(cnt, poly_pt, 0, 0);
341 }
342 delete[] poly_pt;
343 }
344}
345
346void ShapeBaseChart::AddPointToTessList(shp::Point &point, ViewPort &vp,
347 GLUtesselator *tobj, bool idl) {
348#ifdef ocpnUSE_GL
349
350 wxPoint2DDouble q;
351 if (glChartCanvas::HasNormalizedViewPort(vp)) {
352 q = ShapeBaseChartSet::GetDoublePixFromLL(vp, point.getY(), point.getX());
353 } else { // tesselation directly from lat/lon
354 q.m_x = point.getY(), q.m_y = point.getX();
355 }
356 GLvertexshp *vertex = new GLvertexshp();
357 g_vertexesshp.push_back(vertex);
358 if (vp.m_projection_type != PROJECTION_POLAR) {
359 // need to correctly pick +180 or -180 longitude for projections
360 // that have a discontiguous date line
361
362 if (idl && (point.getX() == 180)) {
363 if (vp.m_projection_type == PROJECTION_MERCATOR ||
364 vp.m_projection_type == PROJECTION_EQUIRECTANGULAR) {
365 // q.m_x -= 40058986 * 4096.0; // 360 degrees in normalized
366 // viewport
367 } else {
368 q.m_x -= 360; // lat/lon coordinates
369 }
370 }
371 }
372 vertex->info.x = q.m_x;
373 vertex->info.y = q.m_y;
374
375 gluTessVertex(tobj, (GLdouble *)vertex, (GLdouble *)vertex);
376#endif
377}
378
379void ShapeBaseChart::DoDrawPolygonFilledGL(ocpnDC &pnt, ViewPort &vp,
380 const shp::Feature &feature) {
381#ifdef ocpnUSE_GL
382
383 bool idl =
384 vp.GetBBox().GetMinLon() <= -180 || vp.GetBBox().GetMaxLon() >= 180;
385 auto polygon = static_cast<shp::Polygon *>(feature.getGeometry());
386 for (auto &ring : polygon->getRings()) {
387 size_t cnt{0};
388 GLUtesselator *tobj = gluNewTess();
389
390 gluTessCallback(tobj, GLU_TESS_VERTEX, (_GLUfuncptr)&shpsvertexCallback);
391 gluTessCallback(tobj, GLU_TESS_BEGIN, (_GLUfuncptr)&shpsbeginCallback);
392 gluTessCallback(tobj, GLU_TESS_END, (_GLUfuncptr)&shpsendCallback);
393 gluTessCallback(tobj, GLU_TESS_COMBINE, (_GLUfuncptr)&shpscombineCallback);
394 gluTessCallback(tobj, GLU_TESS_ERROR, (_GLUfuncptr)&shpserrorCallback);
395
396 gluTessNormal(tobj, 0, 0, 1);
397 gluTessProperty(tobj, GLU_TESS_WINDING_RULE, GLU_TESS_WINDING_NONZERO);
398
399 gluTessBeginPolygon(tobj, NULL);
400 gluTessBeginContour(tobj);
401 for (auto &point : ring.getPoints()) {
402 AddPointToTessList(point, vp, tobj, idl);
403 cnt++;
404 }
405 gluTessEndContour(tobj);
406 gluTessEndPolygon(tobj);
407 gluDeleteTess(tobj);
408
409 for (auto ver : g_vertexesshp) delete ver;
410 g_vertexesshp.clear();
411 }
412 float_2Dpt *polyv = new float_2Dpt[g_pvshp.size()];
413 int cnt = 0;
414 for (auto pt : g_pvshp) {
415 polyv[cnt++] = pt;
416 }
417 size_t polycnt = g_pvshp.size();
418 g_pvshp.clear();
419
420 GLuint vbo = 0;
421
422 // Build the shader viewport transform matrix
423 mat4x4 m, mvp;
424 mat4x4_identity(m);
425 mat4x4_scale_aniso(mvp, m, 2.0 / (float)vp.pix_width,
426 2.0 / (float)vp.pix_height, 1.0);
427 mat4x4_translate_in_place(mvp, -vp.pix_width / 2, vp.pix_height / 2, 0);
428
429 float *pvt = new float[2 * (polycnt)];
430 for (size_t i = 0; i < polycnt; i++) {
431 float_2Dpt *pc = polyv + i;
432 // wxPoint2DDouble q(pc->y, pc->x);// = vp.GetDoublePixFromLL(pc->y, pc->x);
433 wxPoint2DDouble q = vp.GetDoublePixFromLL(pc->y, pc->x);
434
435 pvt[i * 2] = q.m_x;
436 pvt[(i * 2) + 1] = q.m_y;
437 }
438
439 GLShaderProgram *shader = pcolor_tri_shader_program[pnt.m_canvasIndex];
440 shader->Bind();
441
442 float colorv[4];
443 colorv[0] = _color.Red() / float(256);
444 colorv[1] = _color.Green() / float(256);
445 colorv[2] = _color.Blue() / float(256);
446 colorv[3] = 1.0;
447 shader->SetUniform4fv("color", colorv);
448
449 shader->SetAttributePointerf("position", pvt);
450
451 glDrawArrays(GL_TRIANGLES, 0, polycnt);
452
453 delete[] pvt;
454 delete[] polyv;
455
456 glDeleteBuffers(1, &vbo);
457 shader->UnBind();
458#endif
459}
460
461void ShapeBaseChart::DrawPolygonFilled(ocpnDC &pnt, ViewPort &vp) {
462 if (!_is_usable) {
463 return;
464 }
465 if (!_reader && !_loading) {
466 _loading = true;
467 _loaded = std::async(std::launch::async, [&]() {
468 bool ret = LoadSHP();
469 _loading = false;
470 return ret;
471 });
472 }
473 if (_loading) {
474 if (_loaded.wait_for(std::chrono::milliseconds(0)) ==
475 std::future_status::ready) {
476 _is_usable = _loaded.get();
477 } else {
478 return; // not yet loaded
479 }
480 }
481
482 int pmod = _dmod;
483 LLBBox bbox = vp.GetBBox();
484
485 // Calculate the starting latitude for the tiles grid
486 // Using floor() aligns with the top-left corner tile model
487 int lat_start = floor(bbox.GetMinLat());
488 if (lat_start < 0)
489 lat_start = lat_start - (pmod + (lat_start % pmod));
490 else
491 lat_start = lat_start - (lat_start % pmod);
492
493 // Calculate the starting longitude for the tiles grid
494 // Using floor() aligns with the top-left corner tile model
495 int lon_start = floor(bbox.GetMinLon());
496 if (lon_start < 0)
497 lon_start = lon_start - (pmod + (lon_start % pmod));
498 else
499 lon_start = lon_start - (lon_start % pmod);
500
501 if (_is_tiled) {
502 for (int i = lat_start; i < ceil(bbox.GetMaxLat()) + pmod; i += pmod) {
503 for (int j = lon_start; j < ceil(bbox.GetMaxLon()) + pmod; j += pmod) {
504 int lon{j};
505 if (j < -180) {
506 lon = j + 360;
507 } else if (j >= 180) {
508 lon = j - 360;
509 }
510 for (auto fid : _tiles[LatLonKey(i, lon)]) {
511 auto const &feature = _reader->getFeature(fid);
512 if (pnt.GetDC()) {
513 DoDrawPolygonFilled(pnt, vp,
514 feature); // Parallelize using std::async?
515 } else {
516 DoDrawPolygonFilledGL(pnt, vp,
517 feature); // Parallelize using std::async?
518 }
519 }
520 }
521 }
522 } else {
523 for (auto const &feature : *_reader) {
524 if (pnt.GetDC()) {
525 DoDrawPolygonFilled(pnt, vp,
526 feature); // Parallelize using std::async?
527 } else {
528 DoDrawPolygonFilledGL(pnt, vp,
529 feature); // Parallelize using std::async?
530 }
531 }
532 }
533}
534
535bool ShapeBaseChart::CrossesLand(double &lat1, double &lon1, double &lat2,
536 double &lon2) {
537 if (!_reader && !_loading) {
538 _loading = true;
539 _loaded = std::async(std::launch::async, [&]() {
540 bool ret = LoadSHP();
541 _loading = false;
542 return ret;
543 });
544 }
545 if (_loading) {
546 if (_loaded.wait_for(std::chrono::milliseconds(0)) ==
547 std::future_status::ready) {
548 _is_usable = _loaded.get();
549 } else {
550 // Chart not yet loaded. Assume no land crossing.
551 return false;
552 }
553 }
554
555 double norm_lon1 = lon1;
556 double norm_lon2 = lon2;
557
558 while (norm_lon1 < -180) norm_lon1 += 360;
559 while (norm_lon1 > 180) norm_lon1 -= 360;
560 while (norm_lon2 < -180) norm_lon2 += 360;
561 while (norm_lon2 > 180) norm_lon2 -= 360;
562
563 // Create normalized coordinate pairs
564 auto A = std::make_pair(lat1, norm_lon1);
565 auto B = std::make_pair(lat2, norm_lon2);
566
567 if (_is_tiled) {
568 // Calculate grid-aligned min/max coordinates based on _dmod
569 double minLat = std::min(lat1, lat2);
570 double maxLat = std::max(lat1, lat2);
571 double minLon = std::min(norm_lon1, norm_lon2);
572 double maxLon = std::max(norm_lon1, norm_lon2);
573
574 // Calculate grid-aligned indices based on the top-left corner model
575 // For latitudes: ceil(lat) gives the top edge of the tile containing that
576 // latitude For longitudes: floor(lon) gives the left edge of the tile
577 // containing that longitude
578 int latmax = ceil(maxLat / _dmod) * _dmod;
579 int latmin = floor(minLat / _dmod) * _dmod;
580 int lonmin = floor(minLon / _dmod) * _dmod;
581 int lonmax = ceil(maxLon / _dmod) * _dmod;
582
583 // Adjust the loop to iterate through all relevant tiles
584 for (int i = latmin; i <= latmax; i += _dmod) {
585 for (int j = lonmin; j < lonmax; j += _dmod) {
586 int lon = j;
587 // Normalize tile longitude to valid range
588 while (lon < -180) lon += 360;
589 while (lon >= 180) lon -= 360;
590
591 // Check if the tile exists in our map
592 auto tileIter = _tiles.find(LatLonKey(i, lon));
593 if (tileIter != _tiles.end()) {
594 for (auto fid : tileIter->second) {
595 auto const &feature = _reader->getFeature(fid);
596 if (PolygonLineIntersect(feature, A, B)) {
597 return true;
598 }
599 }
600 }
601 }
602 }
603 } else {
604 // Non-tiled case: check all features
605 for (auto const &feature : *_reader) {
606 if (PolygonLineIntersect(feature, A, B)) {
607 return true;
608 }
609 }
610 }
611
612 return false;
613}
614
616 if (_loading) {
617 _loading = false;
618 if (_loaded.valid()) {
619 _loaded.wait(); // Wait for async operation to complete
620 }
621 }
622}
623
624bool ShapeBaseChart::LineLineIntersect(const std::pair<double, double> &A,
625 const std::pair<double, double> &B,
626 const std::pair<double, double> &C,
627 const std::pair<double, double> &D) {
628 constexpr double EPSILON = 1e-9; // Tolerance for floating-point comparisons
629
630 // Fast bounding box check - early rejection for improved performance
631 double minABx = std::min(A.first, B.first);
632 double maxABx = std::max(A.first, B.first);
633 double minCDx = std::min(C.first, D.first);
634 double maxCDx = std::max(C.first, D.first);
635
636 if (maxABx < minCDx || maxCDx < minABx) return false;
637
638 double minABy = std::min(A.second, B.second);
639 double maxABy = std::max(A.second, B.second);
640 double minCDy = std::min(C.second, D.second);
641 double maxCDy = std::max(C.second, D.second);
642
643 if (maxABy < minCDy || maxCDy < minABy) return false;
644
645 // Line AB represented as a1x + b1y = c1
646 double a1 = B.second - A.second;
647 double b1 = A.first - B.first;
648 double c1 = a1 * A.first + b1 * A.second;
649
650 // Line CD represented as a2x + b2y = c2
651 double a2 = D.second - C.second;
652 double b2 = C.first - D.first;
653 double c2 = a2 * C.first + b2 * C.second;
654
655 double determinant = a1 * b2 - a2 * b1;
656
657 if (std::abs(determinant) < EPSILON) {
658 // Lines are parallel or collinear.
659 // Check if collinear by testing if points from one segment lie on the other
660 // line.
661 if (std::abs(a1 * C.first + b1 * C.second - c1) < EPSILON &&
662 std::abs(a1 * D.first + b1 * D.second - c1) < EPSILON) {
663 // Segments are collinear, check for overlap
664 return !(maxABx < minCDx || maxCDx < minABx || maxABy < minCDy ||
665 maxCDy < minABy);
666 }
667 return false; // Parallel but not collinear
668 }
669
670 // Calculate intersection point
671 double x = (b2 * c1 - b1 * c2) / determinant;
672 double y = (a1 * c2 - a2 * c1) / determinant;
673
674 // Check if intersection point lies on both segments
675 return (minABx <= x && x <= maxABx && minABy <= y && y <= maxABy &&
676 minCDx <= x && x <= maxCDx && minCDy <= y && y <= maxCDy);
677}
678
679bool ShapeBaseChart::PolygonLineIntersect(const shp::Feature &feature,
680 const std::pair<double, double> &A,
681 const std::pair<double, double> &B) {
682 auto geometry = feature.getGeometry();
683 if (!geometry) {
684 return false;
685 }
686
687 // Try to cast to polygon - if this fails, the geometry isn't a polygon
688 auto polygon = dynamic_cast<shp::Polygon *>(geometry);
689 if (!polygon) {
690 return false;
691 }
692
693 // Calculate line segment bounding box (lat, lon)
694 double minLat = std::min(A.first, B.first);
695 double maxLat = std::max(A.first, B.first);
696 double minLon = std::min(A.second, B.second);
697 double maxLon = std::max(A.second, B.second);
698
699 for (auto &ring : polygon->getRings()) {
700 const auto &points = ring.getPoints();
701 if (points.size() < 3) continue;
702
703 // Quick bounding box check for ring
704 double minRingLat = std::numeric_limits<double>::max();
705 double maxRingLat = std::numeric_limits<double>::lowest();
706 double minRingLon = std::numeric_limits<double>::max();
707 double maxRingLon = std::numeric_limits<double>::lowest();
708
709 // Pre-compute all point pairs to avoid repeated getX/getY calls
710 std::vector<std::pair<double, double>> ringPoints;
711 ringPoints.reserve(points.size());
712
713 for (const auto &point : points) {
714 double lat = point.getY();
715 double lon = point.getX();
716
717 minRingLat = std::min(minRingLat, lat);
718 maxRingLat = std::max(maxRingLat, lat);
719 minRingLon = std::min(minRingLon, lon);
720 maxRingLon = std::max(maxRingLon, lon);
721
722 ringPoints.emplace_back(lat, lon);
723 }
724
725 // Early reject if bounding boxes don't overlap
726 if (maxLat < minRingLat || minLat > maxRingLat || maxLon < minRingLon ||
727 minLon > maxRingLon) {
728 continue;
729 }
730
731 // Check each edge for intersection
732 size_t numPoints = ringPoints.size();
733 std::pair<double, double> previous = ringPoints.back();
734
735 for (size_t i = 0; i < numPoints; i++) {
736 const auto &current = ringPoints[i];
737
738 if (LineLineIntersect(A, B, previous, current)) {
739 return true;
740 }
741
742 previous = current;
743 }
744 }
745
746 return false;
747}
748
749void ShapeBaseChartSet::RenderViewOnDC(ocpnDC &dc, ViewPort &vp) {
750 if (IsUsable()) {
751 ShapeBaseChart &chart = SelectBaseMap(vp.chart_scale);
752 chart.SetColor(land_color);
753 chart.RenderViewOnDC(dc, vp);
754 }
755}
General chart base definitions.
Wrapper class for OpenGL shader programs.
Definition shaders.h:57
A latitude/longitude key for 1x1 or 10x10 degree grid tiles.
Provides platform-specific support utilities for OpenCPN.
Manages a set of ShapeBaseChart objects at different resolutions.
bool IsUsable()
Checks if the chart set contains at least one usable chart.
Represents a basemap chart based on shapefile data.
bool CrossesLand(double &lat1, double &lon1, double &lat2, double &lon2)
Determines if a line segment between two geographical points intersects any land mass represented in ...
void CancelLoading()
Cancel the chart loading operation.
int _dmod
Tile size in degrees.
bool LoadSHP()
Loads the shapefile data into memory.
ViewPort - Core geographic projection and coordinate transformation engine.
Definition viewport.h:84
int pix_height
Height of the viewport in physical pixels.
Definition viewport.h:261
int pix_width
Width of the viewport in physical pixels.
Definition viewport.h:259
wxPoint2DDouble GetDoublePixFromLL(double lat, double lon)
Convert latitude and longitude on the ViewPort to physical pixel coordinates with double precision.
Definition viewport.cpp:147
double chart_scale
Chart scale denominator (e.g., 50000 for a 1:50000 scale).
Definition viewport.h:247
Device context class that can use either wxDC or OpenGL for drawing.
Definition ocpndc.h:60
Enhanced logging interface on top of wx/log.h.