新增通视分析函数;新增第三方库地理距离计算库Erkir,利用该库封装了两个距离函数.

master
cbwu 2 years ago
parent cb3a310221
commit fa6b6d00a4

@ -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\
}
}

@ -1,6 +1,6 @@
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE QtCreatorProject>
<!-- Written by QtCreator 10.0.2, 2023-06-30T15:09:24. -->
<!-- Written by QtCreator 10.0.2, 2023-07-12T13:14:28. -->
<qtcreator>
<data>
<variable>EnvironmentId</variable>
@ -174,7 +174,7 @@
<valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.CustomParsers"/>
<value type="bool" key="ProjectExplorer.BuildConfiguration.ParseStandardOutput">false</value>
<valuelist type="QVariantList" key="ProjectExplorer.BuildConfiguration.UserEnvironmentChanges">
<value type="QString">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</value>
<value type="QString">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</value>
<value type="QString">__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</value>
</valuelist>
<value type="QString" key="ProjectExplorer.ProjectConfiguration.DisplayName">Release</value>

Binary file not shown.

Binary file not shown.

@ -1,18 +1,47 @@
#include "geocomputation.h"
#include "Stream/ugdefs.h"
#include "erkir/sphericalpoint.h"
GeoComputation::GeoComputation()
{
}
/*@brief 计算两点间球面距离
* @param lon1(-180180),lat1(-9090) 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(-9090)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(-9090)lon(-180,180)
* bearing:
* dist:(km)
*/
QPointF GeoComputation::computeOffsetGeoPosition(double lon1, double lat1, double bearing, double dist)
{
/*使用Vincenty's公式求解
* startPoint:lat(-9090)lon(-180,180)
* bearing:
* dist:(km)
*/
QPointF endPoint;
//角度转化为弧度
// qreal lon1 = (fmod(startPoint.x()+540,360)-180.0)*PI/180;

@ -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

@ -0,0 +1,62 @@
#include "geospatialanalysis.h"
#include <QMessageBox>
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;
}

@ -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

@ -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=<install_path> --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)

@ -0,0 +1,86 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2020 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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 <memory>
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<ellipsoidal::Point> 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

@ -0,0 +1,89 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2018 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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 northsouth
//! 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

@ -0,0 +1,120 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2020 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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

@ -0,0 +1,198 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2018-2020 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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 <memory>
#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<cartesian::Point> 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<Point, double> 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<double, double, double> inverse(const Point &point) const;
double m_height{ 0.0 };
Datum m_datum;
};
} // ellipsoidal
} // erkir
#endif // ELLIPSOIDAL_POINT_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_

@ -0,0 +1,70 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2018-2019 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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 <vector>
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

@ -0,0 +1,285 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2018-2019 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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 <vector>
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<Point> polygon = {{0,0}, {1,0}, {0,1}};
auto area = Point::areaOf(polygon); // 6.18e9 mІ
*/
static double areaOf(const std::vector<Point> &polygon, double radius = 6371e3);
};
} // spherical
} // erkir
#endif // SPHERICAL_POINT_H

@ -0,0 +1,156 @@
/**********************************************************************************
* MIT License *
* *
* Copyright (c) 2018-2019 Vahan Aghajanyan <vahancho@gmail.com> *
* *
* 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

Binary file not shown.

Binary file not shown.

@ -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:"<<s;
QString demName = "ASTGTMV003_N32E119_dem";
qMapControl->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:"<<result->bVisible;
qDebug()<<"*******************x:"<<result->pnt3D.x;
qDebug()<<"*******************y:"<<result->pnt3D.y;
qDebug()<<"*******************z:"<<result->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()<<"****************:" <<QDir("./").absolutePath();
UGDatasetVector* dv = (UGDatasetVector*) m_pWorkspace->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:"<<pSceneControl->GetUGLayer3Ds()->GetInnerCount();
}
}
#pragma endregion }
/****************地图量算操作**********************/

@ -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)
//{

Loading…
Cancel
Save