34#include "OCPNPlatform.h"
35#include "shapefile_basemap.h"
37#include "glChartCanvas.h"
46#define __CALL_CONVENTION
48#define __CALL_CONVENTION
52extern wxString gWorldShapefileLocation;
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;
75void __CALL_CONVENTION shpscombineCallback(GLdouble coords[3],
76 GLdouble *vertex_data[4],
81 vertex =
new GLvertexshp();
82 g_vertexesshp.push_back(vertex);
84 vertex->info.x = coords[0];
85 vertex->info.y = coords[1];
87 *dataOut = vertex->data;
90void __CALL_CONVENTION shpserrorCallback(GLenum errorCode) {
91 const GLubyte *estring;
92 estring = gluErrorString(errorCode);
96void __CALL_CONVENTION shpsbeginCallback(GLenum type) {
99 case GL_TRIANGLE_STRIP:
100 case GL_TRIANGLE_FAN:
104 printf(
"tess unhandled begin type: %d\n", type);
110void __CALL_CONVENTION shpsendCallback() {}
112void __CALL_CONVENTION shpsvertexCallback(GLvoid *arg) {
114 vertex = (GLvertexshp *)arg;
116 p.y = vertex->info.x;
117 p.x = vertex->info.y;
120 if (g_typeshp != GL_TRIANGLES) {
122 g_pvshp.push_back(g_p1shp);
123 g_pvshp.push_back(g_p2shp);
126 if (g_typeshp == GL_TRIANGLE_STRIP)
128 else if (g_posshp == 0)
133 g_pvshp.push_back(p);
138ShapeBaseChartSet::ShapeBaseChartSet() : _loaded(false) {
139 land_color = wxColor(170, 175, 80);
141void ShapeBaseChartSet::SetBasemapLandColor(wxColor color) {
144wxColor ShapeBaseChartSet::GetBasemapLandColor() {
return land_color; }
146wxPoint2DDouble ShapeBaseChartSet::GetDoublePixFromLL(
ViewPort &vp,
double lat,
149 p.m_x -= vp.rv_rect.x, p.m_y -= vp.rv_rect.y;
154 if (_basemap_map.find(Quality::crude) != _basemap_map.end()) {
155 return _basemap_map.at(Quality::crude);
157 if (_basemap_map.find(Quality::low) != _basemap_map.end()) {
158 return _basemap_map.at(Quality::low);
160 if (_basemap_map.find(Quality::medium) != _basemap_map.end()) {
161 return _basemap_map.at(Quality::medium);
163 if (_basemap_map.find(Quality::high) != _basemap_map.end()) {
164 return _basemap_map.at(Quality::high);
166 return _basemap_map.at(Quality::full);
170 if (_basemap_map.find(Quality::full) != _basemap_map.end()) {
171 return _basemap_map.at(Quality::full);
173 if (_basemap_map.find(Quality::high) != _basemap_map.end()) {
174 return _basemap_map.at(Quality::high);
176 if (_basemap_map.find(Quality::medium) != _basemap_map.end()) {
177 return _basemap_map.at(Quality::medium);
179 if (_basemap_map.find(Quality::low) != _basemap_map.end()) {
180 return _basemap_map.at(Quality::low);
182 return _basemap_map.at(Quality::crude);
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);
203 return LowestQualityBaseMap();
206void ShapeBaseChartSet::Reset() {
208 wxString basemap_dir;
209 if (gWorldShapefileLocation.empty()) {
210 basemap_dir = g_Platform->GetSharedDataDir();
211 basemap_dir.Append(
"basemap_shp");
213 basemap_dir = gWorldShapefileLocation;
216 LoadBasemaps(basemap_dir.ToStdString());
218void ShapeBaseChartSet::LoadBasemaps(
const std::string &dir) {
220 _basemap_map.clear();
222 if (fs::exists(ShapeBaseChart::ConstructPath(dir,
"crude_10x10"))) {
223 auto c =
ShapeBaseChart(ShapeBaseChart::ConstructPath(dir,
"crude_10x10"),
224 300000000, land_color);
226 _basemap_map.insert(std::make_pair(Quality::crude, c));
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)));
234 if (fs::exists(ShapeBaseChart::ConstructPath(dir,
"medium"))) {
235 _basemap_map.insert(std::make_pair(
237 ShapeBaseChart(ShapeBaseChart::ConstructPath(dir,
"medium"), 1000000,
240 if (fs::exists(ShapeBaseChart::ConstructPath(dir,
"high"))) {
241 _basemap_map.insert(std::make_pair(
243 ShapeBaseChart(ShapeBaseChart::ConstructPath(dir,
"high"), 300000,
246 if (fs::exists(ShapeBaseChart::ConstructPath(dir,
"full"))) {
247 _basemap_map.insert(std::make_pair(
258 if (!fs::exists(_filename)) {
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";
273 auto bounds = temp_reader->getBounds();
274 _is_usable = temp_reader->getCount() > 1 && bounds.getMaxX() <= 180 &&
275 bounds.getMinX() >= -180 && bounds.getMinY() >= -90 &&
283 _is_usable &= temp_reader->getGeometryType() == shp::GeometryType::Polygon;
286 for (
auto field : temp_reader->getFields()) {
287 if (field.getName() ==
"x") {
289 }
else if (field.getName() ==
"y") {
293 _is_tiled = (has_x && has_y);
294 if (_is_usable && _is_tiled) {
296 for (
auto const &feature : *temp_reader) {
302 auto f1 = feature.getAttributes();
305 _tiles[
LatLonKey(std::any_cast<int>(feature.getAttributes()[
"y"]),
306 std::any_cast<int>(feature.getAttributes()[
"x"]))]
312 _reader = temp_reader.release();
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()];
325 auto bbox = vp.GetBBox();
326 for (
auto &point : ring.getPoints()) {
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);
340 pnt.DrawPolygonTessellated(cnt, poly_pt, 0, 0);
346void ShapeBaseChart::AddPointToTessList(shp::Point &point,
ViewPort &vp,
347 GLUtesselator *tobj,
bool idl) {
351 if (glChartCanvas::HasNormalizedViewPort(vp)) {
352 q = ShapeBaseChartSet::GetDoublePixFromLL(vp, point.getY(), point.getX());
354 q.m_x = point.getY(), q.m_y = point.getX();
356 GLvertexshp *vertex =
new GLvertexshp();
357 g_vertexesshp.push_back(vertex);
358 if (vp.m_projection_type != PROJECTION_POLAR) {
362 if (idl && (point.getX() == 180)) {
363 if (vp.m_projection_type == PROJECTION_MERCATOR ||
364 vp.m_projection_type == PROJECTION_EQUIRECTANGULAR) {
372 vertex->info.x = q.m_x;
373 vertex->info.y = q.m_y;
375 gluTessVertex(tobj, (GLdouble *)vertex, (GLdouble *)vertex);
379void ShapeBaseChart::DoDrawPolygonFilledGL(
ocpnDC &pnt,
ViewPort &vp,
380 const shp::Feature &feature) {
384 vp.GetBBox().GetMinLon() <= -180 || vp.GetBBox().GetMaxLon() >= 180;
385 auto polygon =
static_cast<shp::Polygon *
>(feature.getGeometry());
386 for (
auto &ring : polygon->getRings()) {
388 GLUtesselator *tobj = gluNewTess();
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);
396 gluTessNormal(tobj, 0, 0, 1);
397 gluTessProperty(tobj, GLU_TESS_WINDING_RULE, GLU_TESS_WINDING_NONZERO);
399 gluTessBeginPolygon(tobj, NULL);
400 gluTessBeginContour(tobj);
401 for (
auto &point : ring.getPoints()) {
402 AddPointToTessList(point, vp, tobj, idl);
405 gluTessEndContour(tobj);
406 gluTessEndPolygon(tobj);
409 for (
auto ver : g_vertexesshp) delete ver;
410 g_vertexesshp.clear();
412 float_2Dpt *polyv =
new float_2Dpt[g_pvshp.size()];
414 for (
auto pt : g_pvshp) {
417 size_t polycnt = g_pvshp.size();
425 mat4x4_scale_aniso(mvp, m, 2.0 / (
float)vp.
pix_width,
429 float *pvt =
new float[2 * (polycnt)];
430 for (
size_t i = 0; i < polycnt; i++) {
431 float_2Dpt *pc = polyv + i;
436 pvt[(i * 2) + 1] = q.m_y;
439 GLShaderProgram *shader = pcolor_tri_shader_program[pnt.m_canvasIndex];
443 colorv[0] = _color.Red() / float(256);
444 colorv[1] = _color.Green() / float(256);
445 colorv[2] = _color.Blue() / float(256);
447 shader->SetUniform4fv(
"color", colorv);
449 shader->SetAttributePointerf(
"position", pvt);
451 glDrawArrays(GL_TRIANGLES, 0, polycnt);
456 glDeleteBuffers(1, &vbo);
461void ShapeBaseChart::DrawPolygonFilled(
ocpnDC &pnt,
ViewPort &vp) {
465 if (!_reader && !_loading) {
467 _loaded = std::async(std::launch::async, [&]() {
474 if (_loaded.wait_for(std::chrono::milliseconds(0)) ==
475 std::future_status::ready) {
476 _is_usable = _loaded.get();
483 LLBBox bbox = vp.GetBBox();
487 int lat_start = floor(bbox.GetMinLat());
489 lat_start = lat_start - (pmod + (lat_start % pmod));
491 lat_start = lat_start - (lat_start % pmod);
495 int lon_start = floor(bbox.GetMinLon());
497 lon_start = lon_start - (pmod + (lon_start % pmod));
499 lon_start = lon_start - (lon_start % pmod);
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) {
507 }
else if (j >= 180) {
510 for (
auto fid : _tiles[
LatLonKey(i, lon)]) {
511 auto const &feature = _reader->getFeature(fid);
513 DoDrawPolygonFilled(pnt, vp,
516 DoDrawPolygonFilledGL(pnt, vp,
523 for (
auto const &feature : *_reader) {
525 DoDrawPolygonFilled(pnt, vp,
528 DoDrawPolygonFilledGL(pnt, vp,
537 if (!_reader && !_loading) {
539 _loaded = std::async(std::launch::async, [&]() {
546 if (_loaded.wait_for(std::chrono::milliseconds(0)) ==
547 std::future_status::ready) {
548 _is_usable = _loaded.get();
555 double norm_lon1 = lon1;
556 double norm_lon2 = lon2;
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;
564 auto A = std::make_pair(lat1, norm_lon1);
565 auto B = std::make_pair(lat2, norm_lon2);
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);
584 for (
int i = latmin; i <= latmax; i +=
_dmod) {
585 for (
int j = lonmin; j < lonmax; j +=
_dmod) {
588 while (lon < -180) lon += 360;
589 while (lon >= 180) lon -= 360;
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)) {
605 for (
auto const &feature : *_reader) {
606 if (PolygonLineIntersect(feature, A, B)) {
618 if (_loaded.valid()) {
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;
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);
636 if (maxABx < minCDx || maxCDx < minABx)
return false;
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);
643 if (maxABy < minCDy || maxCDy < minABy)
return false;
646 double a1 = B.second - A.second;
647 double b1 = A.first - B.first;
648 double c1 = a1 * A.first + b1 * A.second;
651 double a2 = D.second - C.second;
652 double b2 = C.first - D.first;
653 double c2 = a2 * C.first + b2 * C.second;
655 double determinant = a1 * b2 - a2 * b1;
657 if (std::abs(determinant) < EPSILON) {
661 if (std::abs(a1 * C.first + b1 * C.second - c1) < EPSILON &&
662 std::abs(a1 * D.first + b1 * D.second - c1) < EPSILON) {
664 return !(maxABx < minCDx || maxCDx < minABx || maxABy < minCDy ||
671 double x = (b2 * c1 - b1 * c2) / determinant;
672 double y = (a1 * c2 - a2 * c1) / determinant;
675 return (minABx <= x && x <= maxABx && minABy <= y && y <= maxABy &&
676 minCDx <= x && x <= maxCDx && minCDy <= y && y <= maxCDy);
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();
688 auto polygon =
dynamic_cast<shp::Polygon *
>(geometry);
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);
699 for (
auto &ring : polygon->getRings()) {
700 const auto &points = ring.getPoints();
701 if (points.size() < 3)
continue;
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();
710 std::vector<std::pair<double, double>> ringPoints;
711 ringPoints.reserve(points.size());
713 for (
const auto &point : points) {
714 double lat = point.getY();
715 double lon = point.getX();
717 minRingLat = std::min(minRingLat, lat);
718 maxRingLat = std::max(maxRingLat, lat);
719 minRingLon = std::min(minRingLon, lon);
720 maxRingLon = std::max(maxRingLon, lon);
722 ringPoints.emplace_back(lat, lon);
726 if (maxLat < minRingLat || minLat > maxRingLat || maxLon < minRingLon ||
727 minLon > maxRingLon) {
732 size_t numPoints = ringPoints.size();
733 std::pair<double, double> previous = ringPoints.back();
735 for (
size_t i = 0; i < numPoints; i++) {
736 const auto ¤t = ringPoints[i];
738 if (LineLineIntersect(A, B, previous, current)) {
752 chart.SetColor(land_color);
753 chart.RenderViewOnDC(dc, vp);
General chart base definitions.
Wrapper class for OpenGL shader programs.
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.
int pix_height
Height of the viewport in physical pixels.
int pix_width
Width of the viewport in physical pixels.
wxPoint2DDouble GetDoublePixFromLL(double lat, double lon)
Convert latitude and longitude on the ViewPort to physical pixel coordinates with double precision.
double chart_scale
Chart scale denominator (e.g., 50000 for a 1:50000 scale).
Device context class that can use either wxDC or OpenGL for drawing.
Enhanced logging interface on top of wx/log.h.