diff --git a/MapDisplay.pro b/MapDisplay.pro index 1ae3afd6..ad85b646 100644 --- a/MapDisplay.pro +++ b/MapDisplay.pro @@ -18,6 +18,7 @@ SOURCES += \ geocomputation.cpp \ geofeatureoperator.cpp \ geofileparser.cpp \ + geospatialanalysis.cpp \ icons.cpp \ importscenedatadialog.cpp \ layeroperator.cpp \ @@ -57,6 +58,7 @@ HEADERS += \ geocomputation.h \ geofeatureoperator.h \ geofileparser.h \ + geospatialanalysis.h \ icons.h \ importscenedatadialog.h \ layeroperator.h \ @@ -114,6 +116,7 @@ DEFINES += _HAS_STD_BYTE=0 INCLUDEPATH += "../../include" INCLUDEPATH += "../../include/private" INCLUDEPATH += "../../sample/extensions4Qt" +INCLUDEPATH += "./include" #附加库(目录以L标记,附加库名以l标记) win32{ @@ -125,6 +128,7 @@ win32{ contains(QMAKE_HOST.arch, x86_64){ DESTDIR = "../../sample/debug/x64" LIBS += -L"../../sample/debug/x64" -lExtensions4Qt + LIBS += -L"$$PWD/lib/libd_x64" -lerkird LIBPATH = "../../lib/libd_x64" } LIBS += -lSuToolkitd \ @@ -163,6 +167,7 @@ win32{ -lSuSymbolMarker3Dd\ -lSuLayer3DFiled\ -lSuLayer3DTreed\ + -lSuGridAnalystd\ # -lSuFileParserGLTFd\ }else:CONFIG(release, debug|release){ @@ -170,10 +175,12 @@ win32{ contains(QMAKE_HOST.arch, x86_64){ DESTDIR = "../../sample/release/x64" LIBS += -L"../../sample/release/x64" -lExtensions4Qt + LIBS += -L"$$PWD/lib/lib_x64" -lerkir LIBPATH = "../../lib/lib_x64" }else{ DESTDIR = "../release/x86" LIBS += -L"../release/x86" -lExtensions4Qt + LIBS += -L"$$PWD/lib/lib_x64" -lerkir LIBPATH = "../../lib/lib" } LIBS += -lSuToolkit \ @@ -211,6 +218,7 @@ win32{ -lSuSymbolMarker3D\ -lSuLayer3DFile\ -lSuLayer3DTree\ + -lSuGridAnalyst\ # -lSuFileParserGLTF\ } } diff --git a/MapDisplay.pro.user b/MapDisplay.pro.user index a7f7efeb..b4847440 100644 --- a/MapDisplay.pro.user +++ b/MapDisplay.pro.user @@ -1,6 +1,6 @@ - + EnvironmentId @@ -174,7 +174,7 @@ false - Path=D:\Program Files\Microsoft Visual Studio\2022\Professional\VC\Tools\MSVC\14.36.32532\bin\HostX64\x64;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\VC\VCPackages;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\TestWindow;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\TeamFoundation\Team Explorer;D:\Program Files\Microsoft Visual Studio\2022\Professional\MSBuild\Current\bin\Roslyn;D:\Program Files\Microsoft Visual Studio\2022\Professional\Team Tools\Performance Tools\x64;D:\Program Files\Microsoft Visual Studio\2022\Professional\Team Tools\Performance Tools;C:\Program Files (x86)\Windows Kits\10\bin\10.0.22000.0\\x64;C:\Program Files (x86)\Windows Kits\10\bin\\x64;D:\Program Files\Microsoft Visual Studio\2022\Professional\\MSBuild\Current\Bin\amd64;C:\Windows\Microsoft.NET\Framework64\v4.0.30319;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\Tools\;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\bin;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\libnvvp;D:\Program Files (x86)\VMware\VMware Workstation\bin\;C:\Program Files (x86)\Common Files\Oracle\Java\javapath;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Windows\System32\OpenSSH\;D:\Java\jdk1.8.0_361\bin;D:\Java\jdk1.8.0_361\jre\bin;D:\Program Files\nodejs\node_global;D:\Qt\5.15.2\msvc2019_64\bin;D:\Qt\5.15.2\mingw81_64\bin;D:\Qt\Tools\mingw810_64\bin;D:\Anaconda3;D:\Anaconda3\Scripts;D:\Anaconda3\Library\bin;C:\Program Files (x86)\NVIDIA Corporation\PhysX\Common;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\include;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\lib;D:\Program Files\CMake\bin;D:\Program Files\nodejs\;D:\Program Files\Git\cmd;C:\Users\wuche\AppData\Local\Microsoft\WindowsApps;D:\Program Files\Microsoft VS Code\bin;C:\Users\wuche\AppData\Roaming\npm;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\CMake\CMake\bin;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\CMake\Ninja;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\VC\Linux\bin\ConnectionManagerExe;D:\supermap-iobjectscpp\bin\bin_x64;D:\supermap-iobjectscpp\bin\bind_x64 + Path=D:\Program Files\Microsoft Visual Studio\2022\Professional\VC\Tools\MSVC\14.36.32532\bin\HostX64\x64;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\VC\VCPackages;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\TestWindow;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\TeamFoundation\Team Explorer;D:\Program Files\Microsoft Visual Studio\2022\Professional\MSBuild\Current\bin\Roslyn;D:\Program Files\Microsoft Visual Studio\2022\Professional\Team Tools\Performance Tools\x64;D:\Program Files\Microsoft Visual Studio\2022\Professional\Team Tools\Performance Tools;C:\Program Files (x86)\Windows Kits\10\bin\10.0.22000.0\\x64;C:\Program Files (x86)\Windows Kits\10\bin\\x64;D:\Program Files\Microsoft Visual Studio\2022\Professional\\MSBuild\Current\Bin\amd64;C:\Windows\Microsoft.NET\Framework64\v4.0.30319;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\Tools\;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\bin;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\libnvvp;D:\Program Files (x86)\VMware\VMware Workstation\bin\;C:\Program Files (x86)\Common Files\Oracle\Java\javapath;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Windows\System32\OpenSSH\;D:\Java\jdk1.8.0_361\bin;D:\Java\jdk1.8.0_361\jre\bin;D:\Program Files\nodejs\node_global;D:\Qt\5.15.2\msvc2019_64\bin;D:\Qt\5.15.2\mingw81_64\bin;D:\Qt\Tools\mingw810_64\bin;D:\Anaconda3;D:\Anaconda3\Scripts;D:\Anaconda3\Library\bin;C:\Program Files (x86)\NVIDIA Corporation\PhysX\Common;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\include;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\lib;D:\Program Files\CMake\bin;D:\Program Files\nodejs\;D:\Program Files\Git\cmd;C:\Users\wuche\AppData\Local\Microsoft\WindowsApps;D:\Program Files\Microsoft VS Code\bin;C:\Users\wuche\AppData\Roaming\npm;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\CMake\CMake\bin;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\CommonExtensions\Microsoft\CMake\Ninja;D:\Program Files\Microsoft Visual Studio\2022\Professional\Common7\IDE\VC\Linux\bin\ConnectionManagerExe;D:\supermap-iobjectscpp\bin\bin_x64;D:\supermap-iobjectscpp\bin\bind_x64;D:\supermap-iobjectscpp\MyProject\MapDisplay\bin\bin_x64;D:\supermap-iobjectscpp\MyProject\MapDisplay\bin\bind_x64 __VSCMD_PREINIT_PATH=C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\bin;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\libnvvp;D:\Program Files (x86)\VMware\VMware Workstation\bin\;C:\Program Files (x86)\Common Files\Oracle\Java\javapath;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem;C:\Windows\System32\WindowsPowerShell\v1.0\;C:\Windows\System32\OpenSSH\;D:\Java\jdk1.8.0_361\bin;D:\Java\jdk1.8.0_361\jre\bin;D:\Program Files\nodejs\node_global;D:\Qt\5.15.2\msvc2019_64\bin;D:\Qt\5.15.2\mingw81_64\bin;D:\Qt\Tools\mingw810_64\bin;D:\Anaconda3;D:\Anaconda3\Scripts;D:\Anaconda3\Library\bin;C:\Program Files (x86)\NVIDIA Corporation\PhysX\Common;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\include;C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.0\lib;D:\Program Files\CMake\bin;D:\Program Files\nodejs\;D:\Program Files\Git\cmd;C:\Users\wuche\AppData\Local\Microsoft\WindowsApps;D:\Program Files\Microsoft VS Code\bin;C:\Users\wuche\AppData\Roaming\npm;D:\supermap-iobjectscpp\bin\bin_x64 Release diff --git a/bin/bin_x64/erkir.dll b/bin/bin_x64/erkir.dll new file mode 100644 index 00000000..b08dd8ca Binary files /dev/null and b/bin/bin_x64/erkir.dll differ diff --git a/bin/bind_x64/erkird.dll b/bin/bind_x64/erkird.dll new file mode 100644 index 00000000..78e0cb82 Binary files /dev/null and b/bin/bind_x64/erkird.dll differ diff --git a/geocomputation.cpp b/geocomputation.cpp index 01dfd0ee..62372931 100644 --- a/geocomputation.cpp +++ b/geocomputation.cpp @@ -1,18 +1,47 @@ #include "geocomputation.h" #include "Stream/ugdefs.h" +#include "erkir/sphericalpoint.h" + GeoComputation::GeoComputation() { } +/*@brief 计算两点间球面距离 + * @param lon1(-180至180),lat1(-90到90) 经纬度点1(度) + * @param lon2,lat3 经纬度点2(度) + * @return 距离（m） + */ +double GeoComputation::getSphericalDistance(double lon1, double lat1, double lon2, double lat2) +{ + erkir::spherical::Point p1{ lat1, lon1 }; + erkir::spherical::Point p2{ lat2, lon2 }; + double distance = p1.distanceTo(p2); + return distance; +} + +/*@brief 根据起点坐标、方位角、距离，计算另一点坐标。 +* startPoint:起始点地理坐标点（lat(-90到90)，lon(-180,180)） +* bearing:方位角(0-360)（度） +* dist:两点之间距离(m) +*/ +QPointF GeoComputation::getDestinationPoint(double lon1, double lat1, double bearing, double dist) +{ + erkir::spherical::Point p{ lat1, lon1 }; + auto dest = p.destinationPoint(dist, bearing); // 51.5135°N, 000.0983°W + QPointF point(dest.longitude().degrees(),dest.latitude().degrees()); + return point; +} + +/*@brief 根据起点坐标、方位角、距离，计算另一点坐标。 +* 使用Vincenty's公式求解,使用WGS-84椭球 +* startPoint:起始点地理坐标点（lat(-90到90)，lon(-180,180)） +* bearing:方位角（度） +* dist:两点之间距离(km) +*/ QPointF GeoComputation::computeOffsetGeoPosition(double lon1, double lat1, double bearing, double dist) { - /*使用Vincenty's公式求解 - * startPoint:起始点地理坐标点（lat(-90到90)，lon(-180,180)） - * bearing:方位角（度） - * dist:两点之间距离(km) - */ QPointF endPoint; //角度转化为弧度 // qreal lon1 = (fmod(startPoint.x()+540,360)-180.0)*PI/180; diff --git a/geocomputation.h b/geocomputation.h index a0015fc2..e1386d81 100644 --- a/geocomputation.h +++ b/geocomputation.h @@ -11,9 +11,16 @@ class GeoComputation public: GeoComputation(); + //计算两点间球面距离 + double getSphericalDistance(double lon1, double lat1,double lon2, double lat2); + + //根据起始点、距离、方向角计算目标点坐标 + QPointF getDestinationPoint(double lon1, double lat1, double bearing, double dist); + QPointF computeOffsetGeoPosition(double lon1, double lat1, double bearing, double dist); UGPoint2D DPtoMP(UGMap* pMap, const QPoint &point); //显示坐标转地图坐标 + }; #endif // GEOCOMPUTATION_H diff --git a/geospatialanalysis.cpp b/geospatialanalysis.cpp new file mode 100644 index 00000000..205223ec --- /dev/null +++ b/geospatialanalysis.cpp @@ -0,0 +1,62 @@ + +#include "geospatialanalysis.h" +#include + +GeoSpatialAnalysis::GeoSpatialAnalysis() +{ + +} + +/*@brief 两点间的可视性 + * @param pntView 观察点 + * @param pntObject 目标点 + */ +UG3DAnalyst::SingleResult *GeoSpatialAnalysis::DoublePointVisibility(QMapControl* pMapControl,UGPoint2D pntView, UGPoint2D pntObject) +{ + //通过下面的方法首先判断点所在的位置是否存在高程数据； + UGDatasetRasterPtr viewPtDataset = pMapControl->getDemDataSet(pntView); + UGDatasetRasterPtr objectPtDataset = pMapControl->getDemDataSet(pntObject); + + if(viewPtDataset == NULL || objectPtDataset == NULL) + { + QMessageBox::about(NULL,"提示","无高程数据 "); + return NULL;//-1 + } + + if(viewPtDataset != objectPtDataset) + { + QMessageBox::about(NULL,"提示","观测点与被观测点不属于同一个高程数据集 "); + return NULL;//-1 + } + + UGDatasetRasterPtr datasetRaster = viewPtDataset; + + //3D分析分析实例对象 + UG3DAnalyst visAnalsyt; + UGPoint3D pntView3D,pntObject3D; + + //高程值需根据x,y查询数据集得到 + //得到当前测量点的高程值 + // 1.首先将坐标点转换为栅格数据集对应的行和列 + UGPoint pntviewImg; + double altitude; + datasetRaster->XYToImg(pntView,pntviewImg); + altitude = datasetRaster->GetValue(pntviewImg.x,pntviewImg.y); + pntView3D.x = pntView.x; + pntView3D.y = pntView.y; + pntView3D.z = altitude; + + UGPoint pntObjectImg; + //坐标点转换为栅格数据集对应的行和列 + datasetRaster->XYToImg(pntObject,pntObjectImg); + altitude = datasetRaster->GetValue(pntObjectImg.x,pntObjectImg.y); + pntObject3D.x = pntObject.x; + pntObject3D.y = pntObject.y; + pntObject3D.z = altitude; + + //进行通视分析查询； + UG3DAnalyst::SingleResult* result = visAnalsyt.IsVisible(datasetRaster,pntView3D,pntObject3D); + + return result; + +} diff --git a/geospatialanalysis.h b/geospatialanalysis.h new file mode 100644 index 00000000..65a83bc1 --- /dev/null +++ b/geospatialanalysis.h @@ -0,0 +1,32 @@ +#ifndef GEOSPATIALANALYSIS_H +#define GEOSPATIALANALYSIS_H +/* + * @cbwu + * @brief 空间分析类，实现各种空间分析算法 + */ + +#include "qmapcontrol.h" +#include "Engine/UGDataset.h" +#include "GridAnalyst/UG3DAnalyst.h" + +using namespace UGC; + +//MSVC编译器界面显示乱码问题 +#if _MSC_VER >= 1600 +#pragma execution_character_set("utf-8") +#endif + +class GeoSpatialAnalysis +{ +public: + GeoSpatialAnalysis(); + + //两点间的可视性 + UG3DAnalyst::SingleResult* DoublePointVisibility(QMapControl* pMapControl,UGPoint2D pntView, UGPoint2D pntObject); + +private: + + +}; + +#endif // GEOSPATIALANALYSIS_H diff --git a/include/erkir/README.md b/include/erkir/README.md new file mode 100644 index 00000000..2c22a23d --- /dev/null +++ b/include/erkir/README.md @@ -0,0 +1,148 @@ +# Երկիր (Erkir) - a C++ library for geodesic and trigonometric calculations + +[Erkir (armenian: Երկիր, means Earth)](https://github.com/vahancho/erkir) - is inspired +by and based on the great work of [Chris Veness](https://github.com/chrisveness), +the owner of the [Geodesy functions](https://github.com/chrisveness/geodesy) +project - provides a set of comprehensive API for geodesic and trigonometric calculations. +I would call it a C++ port of JavaScript functions provided by the mentioned Chris Veness' project, +however I designed the library to be more object oriented. Thus the code is organized a +little bit differently, but the implementation itself is preserved. + +[![Latest release](https://img.shields.io/github/v/release/vahancho/erkir?include_prereleases)](https://github.com/vahancho/erkir/releases) +[![Test (CMake)](https://github.com/vahancho/erkir/actions/workflows/cmake.yml/badge.svg)](https://github.com/vahancho/erkir/actions/workflows/cmake.yml) +[![codecov](https://codecov.io/gh/vahancho/erkir/branch/master/graph/badge.svg)](https://codecov.io/gh/vahancho/erkir) + +### Prerequisites + +There are no special requirements and dependencies except *C++11* compliant compiler. +The class is tested with *gcc 4.8.4* and *MSVC 15.x* (Visual Studio 2017). +The library is written with pure STL without any third party dependencies. +For more details see the CI badges (*GitHub Actions*) above. + +### Installation + +No installation required. Just incorporate header files from the *include/* and +source files from *src/* directories in your project and compile them. All library +classes are in *erkir* namespace. + +#### Integration with `CMake` projects + +However, if you use `CMake` and want to integrate the library into your project +you might want to install it first by invoking a `CMake` command from the build directory: + +``` +cmake --install . --prefix= --config=Release +``` + +Once the library is installed you can use it from in your project by adjusting its +`CMake` script. For example: + +``` +[..] +find_package(erkir REQUIRED) + +add_executable(example main.cpp) +target_link_libraries(example erkir) +[..] +``` + +### The API + +The code is virtually split into three domains (namespaces) that represent spherical +and ellipsoidal geodetic coordinates and cartesian (x/y/z) for geocentric ones: +`erkir::spherical`, `erkir::ellipsoidal` and `erkir::cartesian` correspondingly. +Spherical Earth model based calculations are accurate enough for most cases, however +in order to gain more precise measurements use `erkir::ellipsoidal` classes. + +`erkir::spherical::Point` class implements geodetic point on the basis of a spherical +earth (ignoring ellipsoidal effects). It uses formulae to calculate distances between +two points (using haversine formula), initial bearing from a point, final bearing to a point, etc. + +`erkir::ellipsoidal::Point` class represents geodetic point based on ellipsoidal +earth model. It includes ellipsoid parameters and datums for different coordinate +systems, and methods for converting between them and to Cartesian coordinates. + +`erkir::Vector3d` implements 3-d vector manipulation routines. With this class you +can perform basic operations with the vectors, such as calculate dot (scalar) product +of two vectors, multiply vectors, add and subtract them. + +`erkir::cartesian::Point` implements ECEF (earth-centered earth-fixed) geocentric +cartesian (x/y/z) coordinates. + +### Usage Examples: + +```cpp +#include "sphericalpoint.h" +#include "ellipsoidalpoint.h" + +int main(int argc, char **argv) +{ + // Calculate great-circle distance between two points. + erkir::spherical::Point p1{ 52.205, 0.119 }; + erkir::spherical::Point p2{ 48.857, 2.351 }; + auto d = p1.distanceTo(p2); // 404.3 km + + // Get destination point by given distance (shortest) and bearing from start point. + erkir::spherical::Point p3{ 51.4778, -0.0015 }; + auto dest = p3.destinationPoint(7794.0, 300.7); // 51.5135°N, 000.0983°W + + // Convert a point from one coordinates system to another. + erkir::ellipsoidal::Point pWGS84(51.4778, -0.0016, ellipsoidal::Datum::Type::WGS84); + auto pOSGB = pWGS84.toDatum(ellipsoidal::Datum::Type::OSGB36); // 51.4778°N, 000.0000°E + + // Convert to Cartesian coordinates. + auto cartesian = pWGS84.toCartesianPoint(); + + // Convert Cartesian point to a geodetic one. + auto geoPoint = cartesian->toGeoPoint(); + + return 0; +} +``` + +For more usage examples please refer to the unit tests at `/test/test.cpp` file. + +### Building and Testing + +There are unit tests. You can find them in the *test/* directory. +To run them you have to build and run the test application. For doing that you must invoke the following +commands from the terminal, assuming that compiler and environment are already configured: + +#### Linux (gcc) + +``` +mkdir build && cd build +cmake .. -DCMAKE_BUILD_TYPE=Release -DENABLE_TESTING=True +cmake --build . +ctest +``` + +#### Windows (MSVC Toolchain) + +``` +mkdir build && cd build +cmake .. -DENABLE_TESTING=True -A x64 +cmake --build . --config=Release +ctest -C Release +``` + +For x86 builds use `-A Win32` option instead. + +### Performance Tests + +I measured performance (on Intel Core i5 series processor) for some spherical geodesy +functions (`Point` class). I used similar approach as Chris Veness did in his tests, +i.e. called functions for 5000 random points or pairs of points. And here are my results: + +| Function | Avg. time/calculation (nanoseconds)| +| -------------------- |:----------------------------------:| +| Distance (haversine) | 162 | +| Initial bearing | 190 | +| Destination point | 227 | + +*of course timings are machine dependent* + +## See Also + +* [Movable Type Scripts Latitude/Longitude Calculations Reference](http://www.movable-type.co.uk/scripts/latlong.html) + diff --git a/include/erkir/cartesianpoint.h b/include/erkir/cartesianpoint.h new file mode 100644 index 00000000..bb1c5567 --- /dev/null +++ b/include/erkir/cartesianpoint.h @@ -0,0 +1,86 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2020 Vahan Aghajanyan * +* * +* Geodesy tools for conversions between (historical) datums * +* (c) Chris Veness 2005-2019 * +* www.movable-type.co.uk/scripts/latlong-convert-coords.html * +* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal-datum * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef CARTESIAN_POINT_H +#define CARTESIAN_POINT_H + +#include "datum.h" +#include "vector3d.h" + +#include + +namespace erkir +{ + +namespace ellipsoidal +{ + class Point; +} + +namespace cartesian +{ + +/// ECEF (earth-centered earth-fixed) geocentric Cartesian coordinates. +class ERKIR_EXPORT Point : public Vector3d +{ +public: + /// Creates Cartesian coordinate representing ECEF (earth-centric earth-fixed) point. + /*! + \param x X coordinate in metres (=> 0°N,0°E). + \param y Y coordinate in metres (=> 0°N,90°E). + \param z Z coordinate in metres (=> 90°N). + + \example auto coord = Cartesian(3980581.210, -111.159, 4966824.522); + */ + Point(double x, double y, double z, const ellipsoidal::Datum &datum = {ellipsoidal::Datum::Type::WGS84}); + + /// Converts 'this' (geocentric) cartesian (x/y/z) coordinate to (geodetic) latitude/longitude + /// point( based on the same datum, or WGS84 if unset ). + /*! + \returns {LatLon} Latitude/longitude point defined by cartesian coordinates. + + \example + auto c = cartesian::Point{4027893.924, 307041.993, 4919474.294}; + auto p = c.toGeoPoint(); // 50.7978°N, 004.3592°E + */ + std::unique_ptr toGeoPoint() const; + + /// Converts this point to the \p targetDatum. + Point &toDatum(ellipsoidal::Datum::Type targetDatum); + +private: + ellipsoidal::Datum m_datum; + }; + +} // cartesian + +} // erkir + +#endif // CARTESIAN_POINT_H + diff --git a/include/erkir/coordinate.h b/include/erkir/coordinate.h new file mode 100644 index 00000000..126f5857 --- /dev/null +++ b/include/erkir/coordinate.h @@ -0,0 +1,89 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2018 Vahan Aghajanyan * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef COORDINATE_H +#define COORDINATE_H + +#include "export.h" + +namespace erkir +{ + +//! Implements the geographical coordinate abstraction. +class ERKIR_EXPORT Coordinate +{ +public: + //! Constructs a coordinate by the given decimal degrees. + Coordinate(double degrees); + + //! Returns this coordinate's value in decimal degrees. + double degrees() const; + + //! Returns this coordinate's value in radians. + double radians() const; + + /// A helper function to convert radians to degrees. + static double toDegrees(double radians); + + /// A helper function to convert degrees to radians. + static double toRadians(double degrees); + + static double pi(); + + /// Constrain degrees to range 0..360.0 (e.g. for bearings); -1 => 359, 361 => 1. + static double wrap360(double degrees); + +private: + double m_degrees; +}; + +//! Implements the latitude - geographic coordinate that specifies the north–south +//! position of a point on the Earth's surface. +class ERKIR_EXPORT Latitude : public Coordinate +{ +public: + //! Constructs a latitude object. + /*! + \param degree Decimal degrees in the range from 0° to (+/–)90° + \throws std::out_of_range + */ + Latitude(double degree); +}; + +//! Implements the longitude - the measurement east or west of the prime meridian. +class ERKIR_EXPORT Longitude : public Coordinate +{ +public: + //! Constructs a longitude object. + /*! + \param degree Decimal degrees in the range from 0° to (+/–)180° + \throws std::out_of_range + */ + Longitude(double degree); +}; + +} + +#endif // COORDINATE_H + diff --git a/include/erkir/datum.h b/include/erkir/datum.h new file mode 100644 index 00000000..61e12ab3 --- /dev/null +++ b/include/erkir/datum.h @@ -0,0 +1,120 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2020 Vahan Aghajanyan * +* * +* Geodesy tools for conversions between (historical) datums * +* (c) Chris Veness 2005-2019 * +* www.movable-type.co.uk/scripts/latlong-convert-coords.html * +* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal-datum * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef DATUM_H +#define DATUM_H + +#include "export.h" + +namespace erkir +{ + +namespace cartesian +{ + class Point; +} + +namespace ellipsoidal +{ + +/// Implements a geodetic datum. +/*! + Note that precision of various datums will vary, and WGS-84 (original) is not defined to be + accurate to better than ±1 metre. No transformation should be assumed to be accurate to better + than a meter; for many datums somewhat less. +*/ +class ERKIR_EXPORT Datum +{ +public: + enum class Type + { + ED50, /*!< The older European Datum */ + Irl1975, + NAD27, /*!< The older North American Datum, of which NAD83 was basically a readjustment */ + NAD83, /*!< The North American Datum which is very similar to WGS 84 */ + NTF, /*!< Nouvelle Triangulation Francaise */ + OSGB36, /*!< Of the Ordnance Survey of Great Britain */ + Potsdam, /*!< The local datum of Germany with underlying Bessel ellipsoid */ + TokyoJapan, + WGS72, /*!< 72 of the World Geodetic System */ + WGS84 /*!< 84 of the World Geodetic System */ + }; + + //! Constructs a datum with the given \p type. WGS84 is the default datum. + Datum(Type type = Type::WGS84); + + //! Implements a reference ellipsoid. + struct Ellipsoid + { + enum class Type + { + WGS84, + Airy1830, + AiryModified, + Bessel1841, + Clarke1866, + Clarke1880IGN, + GRS80, + Intl1924, // aka Hayford + WGS72 + }; + + Ellipsoid(double a, double b, double f); + + /// Major axis (a). + double m_a{0.0}; + + /// Minor axis (b). + double m_b{0.0}; + + /// Flattening (f). + double m_f{0.0}; + }; + + /// Returns the reference ellipsoid for this datum. + const Ellipsoid &ellipsoid() const; + + /// Returns the type of this datum. + Type type() const; + + /// Converts the given cartesian \p point to the \p targetDatum. + void toDatum(cartesian::Point &point, Type targetDatum) const; + + /// Compares two datums. + bool operator==(const Datum &other) const; + +private: + Type m_type{Type::WGS84}; +}; + +} // ellipsoidal + +} // erkir + +#endif // DATUM_H diff --git a/include/erkir/ellipsoidalpoint.h b/include/erkir/ellipsoidalpoint.h new file mode 100644 index 00000000..5120ffed --- /dev/null +++ b/include/erkir/ellipsoidalpoint.h @@ -0,0 +1,198 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2018-2020 Vahan Aghajanyan * +* * +* Geodesy tools for conversions between (historical) datums * +* (c) Chris Veness 2005-2019 * +* www.movable-type.co.uk/scripts/latlong-convert-coords.html * +* www.movable-type.co.uk/scripts/geodesy-library.html#latlon-ellipsoidal-datum * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef ELLIPSOIDAL_POINT_H +#define ELLIPSOIDAL_POINT_H + +#include + +#include "point.h" +#include "datum.h" + + +namespace erkir +{ + +namespace cartesian +{ + class Point; +} + +namespace ellipsoidal +{ + +//! Implements geodetic point based on ellipsoidal earth model. +/*! + Includes ellipsoid parameters and datums for different coordinate systems, and methods for + converting between them and to Cartesian coordinates. +*/ +class ERKIR_EXPORT Point : public erkir::Point +{ +public: + //! Constructs a point with the given \p latitude, \p longitude \p height above ellipsoid in metres and \p datum. + Point(const Latitude &latitude, const Longitude &longitude, double height = 0.0, + const Datum &datum = {Datum::Type::WGS84}); + + /// Return the datum. + Datum datum() const; + + /// Returns height above the ellipsoid. + double height() const; + + /// Converts 'this' point's coordinate system to new one. + /*! + \param toDatum Datum this coordinate is to be converted to. + \returns Reference to this point converted to new datum. + + \example + Point pWGS84(51.47788, -0.00147, Datum::Type::WGS84); + auto pOSGB = pWGS84.toDatum(Datum::Type::OSGB36); // 51.4773°N, 000.0001°E + */ + Point &toDatum(Datum::Type toDatum); + + /// Converts 'this' point from (geodetic) coordinates to (geocentric) Cartesian (x/y/z) coordinates. + /*! + \returns Cartesian point equivalent to lat/lon point, with x, y, z in metres from earth centre. + */ + std::unique_ptr toCartesianPoint(); + + /*! + Returns the distance between 'this' point and destination point along a geodesic on the + surface of the ellipsoid, using Vincenty inverse solution. + + \param point Latitude/longitude of destination point. + \returns Distance in metres between points or NaN if failed to converge. + \example + auto p1 = Point(50.06632, -5.71475); + auto p2 = Point(58.64402, -3.07009); + auto d = p1.distanceTo(p2); // 969,954.166 m + */ + double distanceTo(const Point &point) const; + + /*! + Returns the destination point having travelled the given distance along a geodesic given by + initial bearing from 'this' point, using Vincenty direct solution. + + \param distance Distance travelled along the geodesic in metres. + \param initialBearing Initial bearing in degrees from north. + \returns Destination point. + \example + auto p1 = Point(-37.95103, 144.42487); + auto p2 = p1.destinationPoint(54972.271, 306.86816); // 37.6528°S, 143.9265°E + */ + Point destinationPoint(double distance, double initialBearing) const; + + /*! + Returns the initial bearing (forward azimuth) to travel along a geodesic from ‘this’ point to + the given point, using Vincenty inverse solution. + + \param point Latitude/longitude of destination point. + \returns Initial bearing in degrees from north (0°..360°) or NaN if failed to converge. + \example + auto p1 = Point(50.06632, -5.71475); + auto p2 = Point(58.64402, -3.07009); + auto b1 = p1.initialBearingTo(p2); // 9.1419° + */ + double initialBearingTo(const Point &point) const; + + /*! + Returns the final bearing (reverse azimuth) having travelled along a geodesic from 'this' + point to the given point, using Vincenty inverse solution. + + \param point Latitude/longitude of destination point. + \returns Final bearing in degrees from north (0°..360°) or NaN if failed to converge. + + \example + auto p1 = Point(50.06632, -5.71475); + auto p2 = Point(58.64402, -3.07009); + auto b2 = p1.finalBearingTo(p2); // 11.2972° + */ + double finalBearingTo(const Point &point) const; + + /*! + Returns the final bearing (reverse azimuth) having travelled along a geodesic given by initial + bearing for a given distance from 'this' point, using Vincenty direct solution. + + \param distance Distance travelled along the geodesic in metres. + \param initialBearing Initial bearing in degrees from north. + \returns Final bearing in degrees from north (0°..360°). + + \example + auto p1 = Point(-37.95103, 144.42487); + auto b2 = p1.finalBearingOn(54972.271, 306.86816); // 307.1736° + */ + double finalBearingOn(double distance, double initialBearing) const; + +private: + enum class DirectField + { + Point, + FinalBearing + }; + + //!Vincenty direct calculation. + /*! + Ellipsoid parameters are taken from datum of 'this' point. Height is ignored. + + \param distance Distance along bearing in metres. + \param initialBearing Initial bearing in degrees from north. + \returns Object including point (destination point), finalBearing. + \throws Formula failed to converge. + */ + std::tuple direct(double distance, double initialBearing) const; + + enum class InverseField + { + Distance, + InitialBearing, + FinalBearing + }; + + //! Vincenty inverse calculation. + /*! + Ellipsoid parameters are taken from datum of 'this' point. Height is ignored. + + \param point Latitude/longitude of destination point. + \returns Object including distance, initialBearing, finalBearing. + \throws Invalid point. + \throws Points must be on surface of ellipsoid. + \throws Formula failed to converge. + */ + std::tuple inverse(const Point &point) const; + + double m_height{ 0.0 }; + Datum m_datum; +}; + +} // ellipsoidal + +} // erkir + +#endif // ELLIPSOIDAL_POINT_H + diff --git a/include/erkir/export.h b/include/erkir/export.h new file mode 100644 index 00000000..edcb0b3d --- /dev/null +++ b/include/erkir/export.h @@ -0,0 +1,14 @@ +#ifndef __EXPORT_H_ +#define __EXPORT_H_ + +#ifdef _WIN32 + #ifdef MAKEDLL + # define ERKIR_EXPORT __declspec(dllexport) + #else + # define ERKIR_EXPORT __declspec(dllimport) + #endif +#else + #define ERKIR_EXPORT +#endif + +#endif // __EXPORT_H_ diff --git a/include/erkir/point.h b/include/erkir/point.h new file mode 100644 index 00000000..d23bf8a2 --- /dev/null +++ b/include/erkir/point.h @@ -0,0 +1,70 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2018-2019 Vahan Aghajanyan * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef POINT_H +#define POINT_H + +#include "coordinate.h" +#include "export.h" + +#include + +namespace erkir +{ + +/// Base class for all types of geodetic points. +class ERKIR_EXPORT Point +{ +public: + //! Constructs an invalid point object. + Point(); + + //! Constructs a point with the given \p latitude and \p longitude. + Point(const Latitude &latitude, const Longitude &longitude); + + //! Returns the latitude of this point. + const Latitude &latitude() const; + + //! Returns the longitude of this point. + const Longitude &longitude() const; + + //! Returns true if this point is a valid one and false otherwise. + bool isValid() const; + + /// Returns true if two points are equal and false otherwise. + bool operator==(const Point &other) const; + + /// Returns true if two points are different and false otherwise. + bool operator!=(const Point &other) const; + +private: + Latitude m_latitude; + Longitude m_longitude; + bool m_isValid; +}; + +} // erkir + +#endif // POINT_H + diff --git a/include/erkir/sphericalpoint.h b/include/erkir/sphericalpoint.h new file mode 100644 index 00000000..4a411611 --- /dev/null +++ b/include/erkir/sphericalpoint.h @@ -0,0 +1,285 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2018-2019 Vahan Aghajanyan * +* * +* Latitude/longitude spherical geodesy tools (c) Chris Veness 2002-2018 * +* www.movable-type.co.uk/scripts/latlong.html * +* www.movable-type.co.uk/scripts/geodesy/docs/module-latlon-spherical.html * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef SPHERICAL_POINT_H +#define SPHERICAL_POINT_H + +#include "point.h" +#include + +namespace erkir +{ + +namespace spherical +{ + +//! Implements the geographical point. +/*! + All formulae in this class are for calculations on the basis of a spherical earth + (ignoring ellipsoidal effects). +*/ +class ERKIR_EXPORT Point : public erkir::Point +{ +public: + //! Constructs an invalid point object. + Point(); + + //! Constructs a point with the given \p latitude and \p longitude. + Point(const Latitude &latitude, const Longitude &longitude); + + /// Returns the distance from this point to destination point(using haversine formula). + /*! + This function uses calculations on the basis of a spherical earth + (ignoring ellipsoidal effects) – which is accurate enough for most purposes. + It uses the 'haversine' formula to calculate the great-circle distance between + two points – that is, the shortest distance over the earth's surface – giving + an 'as-the-crow-flies' distance between the points. + + \param point Latitude/longitude of destination point. + \param radius (Mean)radius of earth(defaults to radius in 6371e3 metres). + \returns Distance between this point and destination point, in same units as radius. + + \example + Point p1{ 52.205, 0.119 }; + Point p2{ 48.857, 2.351 }; + auto d = p1.distanceTo(p2); // 404.3 km + */ + double distanceTo(const Point &point, double radius = 6371e3) const; + + /// Returns the(initial) bearing from 'this' point to destination point. + /*! + \param point Latitude/longitude of destination point. + \returns Initial bearing in degrees from north. + + \example + Point p1{ 52.205, 0.119 }; + Point p2{ 48.857, 2.351 }; + auto b1 = p1.bearingTo(p2); // 156.2° + */ + double bearingTo(const Point &point) const; + + /// Returns final bearing arriving at destination point from 'this' point. + /*! + Returns final bearing arriving at destination point from 'this' point; the final bearing + will differ from the initial bearing by autoying degrees according to distance and latitude. + + \param point - Latitude/longitude of destination point. + \returns Final bearing in degrees from north. + + \example + Point p1{52.205, 0.119}; + Point p2{48.857, 2.351}; + auto b2 = p1.finalBearingTo(p2); // 157.9° + */ + double finalBearingTo(const Point &point) const; + + /// Returns the midpoint between 'this' point and the supplied point. + /*! + \param point - Latitude/longitude of destination point. + \returns Midpoint between this point and the supplied point. + + \example + Point p1{52.205, 0.119}; + Point p2{48.857, 2.351}; + auto pMid = p1.midpointTo(p2); // 50.5363°N, 001.2746°E + */ + Point midpointTo(const Point &point) const; + + /// Returns the point at given fraction between 'this' point and specified point. + /*! + \param point Latitude/longitude of destination point. + \param fraction Fraction between the two points (0 = this point, 1.0 = specified point). + \returns Intermediate point between this point and destination point. + + \example + Point p1{52.205, 0.119}; + Point p2{48.857, 2.351}; + auto pMid = p1.intermediatePointTo(p2, 0.25); // 51.3721°N, 000.7073°E + */ + Point intermediatePointTo(const Point &point, double fraction) const; + + /// Returns the destination point from 'this' point. + /*! + Returns the destination point from 'this' point having travelled the given distance on the + given initial bearing (bearing normally autoies around path followed). + + \param distance Distance travelled, in same units as earth radius (default: metres). + \param bearing Initial bearing in degrees from north. + \param radius (Mean) radius of earth (defaults to radius in 6371e3 metres). + \returns Destination point. + + \example + Point p1{51.4778, -0.0015}; + Point p2 = p1.destinationPoint(7794, 300.7); // 51.5135°N, 000.0983°W + */ + Point destinationPoint(double distance, double bearing, double radius = 6371e3) const; + + /// Returns the point of intersection of two paths defined by point and bearing. + /*! + @param p1 First point. + @param brng1 Initial bearing from first point in degrees. + @param p2 Second point. + @param brng2 Initial bearing from second point in degrees. + @returns Destination point (an invalid point if no unique intersection defined). + + @example + Point p1{51.8853, 0.2545} + auto brng1 = 108.547; + Point p2{49.0034, 2.5735} + auto brng2 = 32.435; + auto pInt = intersection(p1, brng1, p2, brng2); // 50.9078°N, 004.5084°E + */ + static Point intersection(const Point &p1, double brng1, + const Point &p2, double brng2); + + /// Returns (signed) distance from 'this' point to great circle defined by start - point and end - point. + /*! + \param pathStart Start point of great circle path. + \param pathEnd End point of great circle path. + \param radius (Mean) radius of earth (defaults to radius in 6371e3 metres). + \returns Distance to great circle (negative if to left, positive if to right of path). + + \example + Point pCurrent{53.2611, -0.7972}; + Point p1{53.3206, -1.7297}; + Point p2{53.1887, 0.1334}; + auto d = pCurrent.crossTrackDistanceTo(p1, p2); // -307.5 m + */ + double crossTrackDistanceTo(const Point &pathStart, const Point &pathEnd, + double radius = 6371e3) const; + + /// Returns how far 'this' point is along a path from start-point. + /*! + Returns how far 'this' point is along a path from start-point, heading towards end-point. + That is, if a perpendicular is drawn from 'this' point to the (great circle) path, the along-track + distance is the distance from the start point to where the perpendicular crosses the path. + + \param pathStart Start point of great circle path. + \param pathEnd End point of great circle path. + \param radius (Mean) radius of earth (defaults to radius in 6371e3 metres). + \returns Distance along great circle to point nearest 'this' point. + + \example + Point pCurrent{53.2611, -0.7972}; + Point p1{53.3206, -1.7297}; + Point p2{53.1887, 0.1334}; + auto d = pCurrent.alongTrackDistanceTo(p1, p2); // 62.331 km + */ + double alongTrackDistanceTo(const Point &pathStart, const Point &pathEnd, + double radius = 6371e3) const; + + /// Returns maximum latitude reached when travelling on a great circle. + /*! + Returns maximum latitude reached when travelling on a great circle on given bearing from this + point ('Clairaut's formula'). Negate the result for the minimum latitude (in the Southern + hemisphere). + + The maximum latitude is independent of longitude; it will be the same for all points on a given + latitude. + + \param bearing Initial bearing. + \param latitude Starting latitude. + */ + double maxLatitude(double bearing) const; + + /// Returns the distance travelling from 'this' point to destination point along a rhumb line. + /*! + A 'rhumb line' (or loxodrome) is a path of constant bearing, which crosses + all meridians at the same angle. Sailors used to (and sometimes still) navigate + along rhumb lines since it is easier to follow a constant compass bearing than + to be continually adjusting the bearing, as is needed to follow a great circle. + Rhumb lines are straight lines on a Mercator Projec­tion map (also helpful for + naviga­tion). Rhumb lines are generally longer than great-circle (orthodrome) routes. + + \param point Latitude/longitude of destination point. + \param radius (Mean) radius of earth (defaults to radius in 6371e3 metres). + \returns Distance between this point and destination point (same units as radius). + + \example + Point p1{51.127, 1.338}; + Point p2{50.964, 1.853}; + auto d = p1.rhumbDistanceTo(p2); // 40.31 km + */ + double rhumbDistanceTo(const Point &point, double radius = 6371e3) const; + + /// Returns the bearing from 'this' point to destination point along a rhumb line. + /*! + \param point Latitude/longitude of destination point. + \returns Bearing in degrees from north. + + \example + Point p1{51.127, 1.338}; + Point p2{50.964, 1.853}; + auto d = p1.rhumbBearingTo(p2); // 116.7 + */ + double rhumbBearingTo(const Point &point) const; + + /// Returns the destination point having travelled along a rhumb line from 'this' point the given distance on the given bearing. + /*! + \param distance Distance travelled, in same units as earth radius (default: metres). + \param bearing Bearing in degrees from north. + \param radius (Mean)radius of earth(defaults to radius in 6371e3 metres). + \returns Destination point. + + \example + Point p1{51.127, 1.338}; + auto p2 = p1.rhumbDestinationPoint(40300, 116.7); // 50.9642°N, 001.8530°E + */ + Point rhumbDestinationPoint(double distance, double bearing, double radius = 6371e3) const; + + /// Returns the loxodromic midpoint (along a rhumb line) between 'this' point and second point. + /*! + \param point Latitude/longitude of second point. + \returns Midpoint between this point and second point. + + \example + Point p1{51.127, 1.338}; + Point p2{50.964, 1.853}; + auto pMid = p1.rhumbMidpointTo(p2); // 51.0455°N, 001.5957°E + */ + Point rhumbMidpointTo(const Point &point) const; + + /// Calculates the area of a spherical polygon where the sides of the polygon are great circle arcs joining the vertices. + /*! + \param polygon Array of points defining vertices of the polygon + \param radius (Mean)radius of earth(defaults to radius in 6371e3 metres). + \returns The area of the polygon, in the same units as radius. + + \example + std::vector polygon = {{0,0}, {1,0}, {0,1}}; + auto area = Point::areaOf(polygon); // 6.18e9 mІ + */ + static double areaOf(const std::vector &polygon, double radius = 6371e3); +}; + +} // spherical + +} // erkir + +#endif // SPHERICAL_POINT_H + diff --git a/include/erkir/vector3d.h b/include/erkir/vector3d.h new file mode 100644 index 00000000..59a58eec --- /dev/null +++ b/include/erkir/vector3d.h @@ -0,0 +1,156 @@ +/********************************************************************************** +* MIT License * +* * +* Copyright (c) 2018-2019 Vahan Aghajanyan * +* * +* Vector handling functions (c) Chris Veness 2011-2016 * +* www.movable-type.co.uk/scripts/geodesy/docs/module-vector3d.html * +* * +* Permission is hereby granted, free of charge, to any person obtaining a copy * +* of this software and associated documentation files (the "Software"), to deal * +* in the Software without restriction, including without limitation the rights * +* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * +* copies of the Software, and to permit persons to whom the Software is * +* furnished to do so, subject to the following conditions: * +* * +* The above copyright notice and this permission notice shall be included in all * +* copies or substantial portions of the Software. * +* * +* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * +* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * +* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * +* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * +* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * +* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * +* SOFTWARE. * +***********************************************************************************/ + +#ifndef VECTOR3D_H +#define VECTOR3D_H + +#include "export.h" + +namespace erkir +{ + +/// Implements 3-d vector manipulation routines. +/*! + In a geodesy context, these vectors may be used to represent: + - n-vector representing a normal to point on Earth's surface + - earth-centered, earth fixed vector (= Gade's 'p-vector') + - great circle normal to vector (on spherical earth model) + - motion vector on Earth's surface + - etc + + Functions return vectors as return results, so that operations can be chained. + \example auto v = v1.cross(v2).dot(v3) // equivalent to v1 × v2 . v3 +*/ +class ERKIR_EXPORT Vector3d +{ +public: + /// Creates an invalid 3-d vector. + Vector3d(); + + /// Creates a 3-d vector. + /*! + The vector may be normalised, or use x/y/z values for eg height relative to the sphere or + ellipsoid, distance from earth centre, etc. + + \param x X component of vector. + \param y Y component of vector. + \param z Z component of vector. + */ + Vector3d(double x, double y, double z); + + double x() const; + double y() const; + double z() const; + bool isValid() const; + + /// Dot (scalar) product of two vectors. + double dot(const Vector3d &v) const; + + /// Multiplies vector by the supplied vector using cross (vector) product. + Vector3d cross(const Vector3d &v) const; + + /// Length (magnitude or norm) of 'this' vector + /*! + \returns Magnitude of this vector. + */ + double length() const; + + /// Normalizes a vector to its unit vector + /*! + If the vector is already unit or is zero magnitude, this is a no-op. + \returns Normalised version of this vector. + */ + Vector3d unit() const; + + /// Calculates the angle between 'this' vector and supplied vector. + /*! + \param v Supplied vector + \param n Plane normal: if supplied, angle is -PI..+PI, signed +ve if this->v is + clockwise looking along n, -ve in opposite direction (if not supplied, angle is always 0..PI). + \returns Angle (in radians) between this vector and supplied vector. + */ + double angleTo(const Vector3d &v, const Vector3d &n = Vector3d()) const; + + /// Rotates 'this' point around an axis by a specified angle. + /*! + \param axis The axis being rotated around. + \param theta The angle of rotation (in radians). + \returns The rotated point. + */ + Vector3d rotateAround(const Vector3d &axis, double theta) const; + +private: + double m_x{ 0.0 }; + double m_y{ 0.0 }; + double m_z{ 0.0 }; + + bool m_isValid; +}; + +/// Vector addition +inline Vector3d operator + (const Vector3d &v1, const Vector3d &v2) +{ + return Vector3d(v1.x() + v2.x(), v1.y() + v2.y(), v1.z() + v2.z()); +} + +/// Vector subtraction +inline Vector3d operator - (const Vector3d &v1, const Vector3d & v2) +{ + return Vector3d(v1.x() - v2.x(), v1.y() - v2.y(), v1.z() - v2.z()); +} + +/// Unary negation of a vector +/*! + Negates a vector to point in the opposite direction. +*/ +inline Vector3d operator - (const Vector3d &v1) +{ + return Vector3d(-v1.x(), -v1.y(), -v1.z()); +} + +/// Multiplication of vector and scalar +inline Vector3d operator * (const Vector3d &v1, double s) +{ + return Vector3d(v1.x() * s, v1.y() * s, v1.z() * s); +} + +/// Division of vector and scalar +inline Vector3d operator / (const Vector3d &v1, double s) +{ + return Vector3d(v1.x() / s, v1.y() / s, v1.z() / s); +} + +/// Multiplication of scalar and vector +inline Vector3d operator * (double s, const Vector3d &v1) +{ + return v1 * s; +} + +} // erkir + +#endif // VECTOR3D_H + diff --git a/lib/lib_x64/erkir.lib b/lib/lib_x64/erkir.lib new file mode 100644 index 00000000..1ef55954 Binary files /dev/null and b/lib/lib_x64/erkir.lib differ diff --git a/lib/libd_x64/erkird.lib b/lib/libd_x64/erkird.lib new file mode 100644 index 00000000..80e3e797 Binary files /dev/null and b/lib/libd_x64/erkird.lib differ diff --git a/mainwindow.cpp b/mainwindow.cpp index 8716814d..581104b1 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -3,6 +3,7 @@ #include "FileParser/UGModelConfigParams.h" #include "computeoffsetpositiondialog.h" #include "displayroutedialog.h" +#include "geospatialanalysis.h" #include "importscenedatadialog.h" #include "maplocationdialog.h" #include "qregularexpression.h" @@ -104,6 +105,8 @@ #include "Symbol/UGMarkerSymbolLib.h" +#include "erkir/sphericalpoint.h" + //工作空间及数据源默认存储路径及别名 const QString wsName = "MapWorkspace"; //工作空间别名 QString wsFullName; //工作空间存储路径 @@ -380,11 +383,29 @@ MainWindow::MainWindow(QWidget *parent) bool s = m_pWorkspace->GetResources()->GetMarkerSymbolLib()->IsIDExisted(318); qDebug()<<"***********s2:"<addDEMDataset(demName); + + UGPoint2D ptView(119.703784,32.228956); + UGPoint2D ptObject(119.707676,32.227302); + GeoSpatialAnalysis geoSAnalysis; + UG3DAnalyst::SingleResult* result = geoSAnalysis.DoublePointVisibility(qMapControl,ptView,ptObject); + if(!result==NULL) + { + qDebug()<<"*******************bVisible:"<bVisible; + qDebug()<<"*******************x:"<pnt3D.x; + qDebug()<<"*******************y:"<pnt3D.y; + qDebug()<<"*******************z:"<pnt3D.z; + } + + erkir::spherical::Point p1{ 52.205, 0.119 }; + erkir::spherical::Point p2{ 48.857, 2.351 }; + auto d = p1.distanceTo(p2); // 404.3 km -// int id = ms->GetSymbolAt(0)->GetID(); -// QString name = Translator::UGStr2QStr(ms->GetSymbolAt(0)->GetName()); -// qDebug()<<"****************:"; + erkir::spherical::Point p3{ 51.4778, -0.0015 }; + auto dest = p3.destinationPoint(7794.0, 300.7); // 51.5135°N, 000.0983°W + qDebug()<<"*******************z:"; /* qDebug()<<"****************:" <GetDataSource(0)->GetDataset(line3DSetAlias).get(); @@ -547,8 +568,9 @@ void MainWindow::initMapControl() { //实例化一个地图窗口对象，即QMapControl对象 qMapControl = new QMapControl; - qMapControl->GetUGMapWnd()->SetAfterGeometryAddedFunc(AfterGeometryAddedCallback,(UGlong)qMapControl); + qMapControl->GetMap()->SetWorkspace(m_pWorkspace); + qMapControl->GetUGMapWnd()->SetAfterGeometryAddedFunc(AfterGeometryAddedCallback,(UGlong)qMapControl); qMapControl->GetUGMapWnd()->SetAfterPointInputFunc(AfterPointInputCallback,(UGlong)qMapControl); qMapControl->GetUGMapWnd()->SetGeometrySelectedFunc(GeometrySelectedCallback,(UGlong)qMapControl); // qMapControl->GetUGMapWnd()->SetSingleGeometrySelectedFunc(SingleGeometrySelectedCallback,(UGlong)qMapControl); @@ -2468,8 +2490,8 @@ void MainWindow::Edit() // qDebug()<<"*********************1:"<GetUGLayer3Ds()->GetInnerCount(); } - } + #pragma endregion } /****************地图量算操作**********************/ diff --git a/mapdatamaneger.cpp b/mapdatamaneger.cpp index ad7925b6..7b4b08a7 100644 --- a/mapdatamaneger.cpp +++ b/mapdatamaneger.cpp @@ -57,6 +57,8 @@ UGDatasetVectorPtr MapDataManeger::createVectorSet(UGDataSource *ds, UGDataset:: return ds->CreateDatasetVector(vectorInfo); } + + //UGLayer3D *MapDataManeger::createTemporaryLayer3D(UGDataSource *ds, UGLayer3Ds *pLayers, UGint nUGLayerType, UGString datasetName) //{