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