Remove Osgb36 (UKgrid) support. Such conversions should be done with the proper tools.
This commit is contained in:
Christian Schmitt
parent
c1d0a92ea8
commit
abfea1cd2d
15 changed files with 8 additions and 1079 deletions
|
|
@ -1,8 +1,5 @@
|
||||||
|
|
||||||
|
|
||||||
include_directories(${PROJECT_SOURCE_DIR}/src/Lib)
|
include_directories(${PROJECT_SOURCE_DIR}/src/Lib)
|
||||||
include_directories(${PROJECT_SOURCE_DIR}/src/BuildTiles)
|
include_directories(${PROJECT_SOURCE_DIR}/src/BuildTiles)
|
||||||
|
|
||||||
add_subdirectory(Osgb36)
|
|
||||||
add_subdirectory(Parallel)
|
add_subdirectory(Parallel)
|
||||||
add_subdirectory(Main)
|
add_subdirectory(Main)
|
||||||
|
|
|
||||||
|
|
@ -27,7 +27,6 @@ set_target_properties(tg-construct PROPERTIES
|
||||||
"DEFAULT_USGS_MAPFILE=\"${PKGDATADIR}/usgsmap.txt\";DEFAULT_PRIORITIES_FILE=\"${PKGDATADIR}/default_priorities.txt\"" )
|
"DEFAULT_USGS_MAPFILE=\"${PKGDATADIR}/usgsmap.txt\";DEFAULT_PRIORITIES_FILE=\"${PKGDATADIR}/default_priorities.txt\"" )
|
||||||
|
|
||||||
target_link_libraries(tg-construct
|
target_link_libraries(tg-construct
|
||||||
Osgb36
|
|
||||||
Polygon Geometry
|
Polygon Geometry
|
||||||
Array landcover poly2tri
|
Array landcover poly2tri
|
||||||
${GDAL_LIBRARY}
|
${GDAL_LIBRARY}
|
||||||
|
|
|
||||||
|
|
@ -143,7 +143,6 @@ static void usage( const string name ) {
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " --nudge=<float>");
|
SG_LOG(SG_GENERAL, SG_ALERT, " --nudge=<float>");
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " --priorities=<filename>");
|
SG_LOG(SG_GENERAL, SG_ALERT, " --priorities=<filename>");
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " --usgs-map=<filename>");
|
SG_LOG(SG_GENERAL, SG_ALERT, " --usgs-map=<filename>");
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " --useUKgrid");
|
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " --ignore-landmass");
|
SG_LOG(SG_GENERAL, SG_ALERT, " --ignore-landmass");
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " ] <load directory...>");
|
SG_LOG(SG_GENERAL, SG_ALERT, " ] <load directory...>");
|
||||||
exit(-1);
|
exit(-1);
|
||||||
|
|
@ -166,10 +165,6 @@ int main(int argc, char **argv) {
|
||||||
vector<string> debug_shape_defs;
|
vector<string> debug_shape_defs;
|
||||||
vector<string> debug_area_defs;
|
vector<string> debug_area_defs;
|
||||||
|
|
||||||
// flag indicating whether UK grid should be used for in-UK
|
|
||||||
// texture coordinate generation
|
|
||||||
bool useUKgrid = false;
|
|
||||||
|
|
||||||
bool ignoreLandmass = false;
|
bool ignoreLandmass = false;
|
||||||
|
|
||||||
sglog().setLogLevels( SG_ALL, SG_INFO );
|
sglog().setLogLevels( SG_ALL, SG_INFO );
|
||||||
|
|
@ -208,8 +203,6 @@ int main(int argc, char **argv) {
|
||||||
priorities_file = arg.substr(13);
|
priorities_file = arg.substr(13);
|
||||||
} else if (arg.find("--usgs-map=") == 0) {
|
} else if (arg.find("--usgs-map=") == 0) {
|
||||||
usgs_map_file = arg.substr(11);
|
usgs_map_file = arg.substr(11);
|
||||||
} else if (arg.find("--useUKgrid") == 0) {
|
|
||||||
useUKgrid = true;
|
|
||||||
} else if (arg.find("--ignore-landmass") == 0) {
|
} else if (arg.find("--ignore-landmass") == 0) {
|
||||||
ignoreLandmass = true;
|
ignoreLandmass = true;
|
||||||
} else if (arg.find("--debug-dir=") == 0) {
|
} else if (arg.find("--debug-dir=") == 0) {
|
||||||
|
|
@ -267,7 +260,7 @@ int main(int argc, char **argv) {
|
||||||
all_stages = new TGConstruct();
|
all_stages = new TGConstruct();
|
||||||
all_stages->set_cover( cover );
|
all_stages->set_cover( cover );
|
||||||
all_stages->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
all_stages->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
||||||
all_stages->set_options( useUKgrid, ignoreLandmass, nudge );
|
all_stages->set_options( ignoreLandmass, nudge );
|
||||||
all_stages->set_bucket( b );
|
all_stages->set_bucket( b );
|
||||||
all_stages->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
all_stages->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
||||||
|
|
||||||
|
|
@ -293,7 +286,7 @@ int main(int argc, char **argv) {
|
||||||
all_stages = new TGConstruct();
|
all_stages = new TGConstruct();
|
||||||
all_stages->set_cover( cover );
|
all_stages->set_cover( cover );
|
||||||
all_stages->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
all_stages->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
||||||
all_stages->set_options( useUKgrid, ignoreLandmass, nudge );
|
all_stages->set_options( ignoreLandmass, nudge );
|
||||||
all_stages->set_bucket( b_min );
|
all_stages->set_bucket( b_min );
|
||||||
all_stages->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
all_stages->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
||||||
|
|
||||||
|
|
@ -325,7 +318,7 @@ int main(int argc, char **argv) {
|
||||||
stage1 = new TGConstruct();
|
stage1 = new TGConstruct();
|
||||||
stage1->set_cover( cover );
|
stage1->set_cover( cover );
|
||||||
stage1->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
stage1->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
||||||
stage1->set_options( useUKgrid, ignoreLandmass, nudge );
|
stage1->set_options( ignoreLandmass, nudge );
|
||||||
stage1->set_bucket( b_cur );
|
stage1->set_bucket( b_cur );
|
||||||
stage1->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
stage1->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
||||||
|
|
||||||
|
|
@ -354,7 +347,7 @@ int main(int argc, char **argv) {
|
||||||
stage2 = new TGConstruct();
|
stage2 = new TGConstruct();
|
||||||
stage2->set_cover( cover );
|
stage2->set_cover( cover );
|
||||||
stage2->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
stage2->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
||||||
stage2->set_options( useUKgrid, ignoreLandmass, nudge );
|
stage2->set_options( ignoreLandmass, nudge );
|
||||||
stage2->set_bucket( b_cur );
|
stage2->set_bucket( b_cur );
|
||||||
stage2->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
stage2->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
||||||
|
|
||||||
|
|
@ -384,7 +377,7 @@ int main(int argc, char **argv) {
|
||||||
stage3 = new TGConstruct();
|
stage3 = new TGConstruct();
|
||||||
stage3->set_cover( cover );
|
stage3->set_cover( cover );
|
||||||
stage3->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
stage3->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
||||||
stage3->set_options( useUKgrid, ignoreLandmass, nudge );
|
stage3->set_options( ignoreLandmass, nudge );
|
||||||
stage3->set_bucket( b_cur );
|
stage3->set_bucket( b_cur );
|
||||||
stage3->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
stage3->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
||||||
|
|
||||||
|
|
@ -408,7 +401,7 @@ int main(int argc, char **argv) {
|
||||||
all_stages = new TGConstruct();
|
all_stages = new TGConstruct();
|
||||||
all_stages->set_cover( cover );
|
all_stages->set_cover( cover );
|
||||||
all_stages->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
all_stages->set_paths( work_dir, share_dir, output_dir, load_dirs );
|
||||||
all_stages->set_options( useUKgrid, ignoreLandmass, nudge );
|
all_stages->set_options( ignoreLandmass, nudge );
|
||||||
all_stages->set_bucket( b );
|
all_stages->set_bucket( b );
|
||||||
all_stages->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
all_stages->set_debug( debug_dir, debug_area_defs, debug_shape_defs );
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -50,7 +50,6 @@
|
||||||
|
|
||||||
// Constructor
|
// Constructor
|
||||||
TGConstruct::TGConstruct():
|
TGConstruct::TGConstruct():
|
||||||
useUKGrid(false),
|
|
||||||
ignoreLandmass(false),
|
ignoreLandmass(false),
|
||||||
debug_all(false),
|
debug_all(false),
|
||||||
ds_id((void*)-1)
|
ds_id((void*)-1)
|
||||||
|
|
@ -77,8 +76,7 @@ void TGConstruct::set_paths( const std::string work, const std::string share, co
|
||||||
load_dirs = load;
|
load_dirs = load;
|
||||||
}
|
}
|
||||||
|
|
||||||
void TGConstruct::set_options( bool uk_grid, bool ignore_lm, double n ) {
|
void TGConstruct::set_options( bool ignore_lm, double n ) {
|
||||||
useUKGrid = uk_grid;
|
|
||||||
ignoreLandmass = ignore_lm;
|
ignoreLandmass = ignore_lm;
|
||||||
nudge = n;
|
nudge = n;
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -87,10 +87,6 @@ private:
|
||||||
|
|
||||||
const static double gSnap = 0.00000001; // approx 1 mm
|
const static double gSnap = 0.00000001; // approx 1 mm
|
||||||
|
|
||||||
// flag indicating whether to align texture coords within the UK
|
|
||||||
// with the UK grid
|
|
||||||
bool useUKGrid;
|
|
||||||
|
|
||||||
// flag indicating whether to ignore the landmass
|
// flag indicating whether to ignore the landmass
|
||||||
bool ignoreLandmass;
|
bool ignoreLandmass;
|
||||||
|
|
||||||
|
|
@ -136,7 +132,6 @@ private:
|
||||||
int LoadLandclassPolys( void );
|
int LoadLandclassPolys( void );
|
||||||
// Load Data Helpers
|
// Load Data Helpers
|
||||||
bool load_poly(const std::string& path);
|
bool load_poly(const std::string& path);
|
||||||
bool load_osgb36_poly(const std::string& path);
|
|
||||||
void add_poly(int area, const TGPolygon &poly, std::string material);
|
void add_poly(int area, const TGPolygon &poly, std::string material);
|
||||||
|
|
||||||
// Clip Data
|
// Clip Data
|
||||||
|
|
@ -239,7 +234,7 @@ public:
|
||||||
inline void set_load_dirs( const std::vector<std::string> ld ) { load_dirs = ld; }
|
inline void set_load_dirs( const std::vector<std::string> ld ) { load_dirs = ld; }
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
void set_options( bool uk_grid, bool ignore_lm, double n );
|
void set_options( bool ignore_lm, double n );
|
||||||
|
|
||||||
#if 0
|
#if 0
|
||||||
// UK grid flag
|
// UK grid flag
|
||||||
|
|
|
||||||
|
|
@ -30,7 +30,6 @@
|
||||||
#include <simgear/debug/logstream.hxx>
|
#include <simgear/debug/logstream.hxx>
|
||||||
|
|
||||||
#include <Geometry/poly_support.hxx>
|
#include <Geometry/poly_support.hxx>
|
||||||
#include <Osgb36/osgb36.hxx>
|
|
||||||
|
|
||||||
#include "tgconstruct.hxx"
|
#include "tgconstruct.hxx"
|
||||||
|
|
||||||
|
|
@ -272,102 +271,6 @@ bool TGConstruct::load_poly(const string& path) {
|
||||||
return true;
|
return true;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
// Load a polygon definition file containing osgb36 Eastings and Northings
|
|
||||||
// and convert them to WGS84 Latitude and Longitude
|
|
||||||
bool TGConstruct::load_osgb36_poly(const string& path) {
|
|
||||||
string poly_name;
|
|
||||||
AreaType poly_type;
|
|
||||||
int contours, count, i, j;
|
|
||||||
int hole_flag;
|
|
||||||
double startx, starty, x, y, lastx, lasty;
|
|
||||||
|
|
||||||
SG_LOG( SG_CLIPPER, SG_INFO, "Loading " << path << " ..." );
|
|
||||||
|
|
||||||
sg_gzifstream in( path );
|
|
||||||
|
|
||||||
if ( !in ) {
|
|
||||||
SG_LOG( SG_CLIPPER, SG_ALERT, "Cannot open file: " << path );
|
|
||||||
exit(-1);
|
|
||||||
}
|
|
||||||
|
|
||||||
TGPolygon poly;
|
|
||||||
|
|
||||||
Point3D p;
|
|
||||||
Point3D OSRef;
|
|
||||||
in >> skipcomment;
|
|
||||||
while ( !in.eof() ) {
|
|
||||||
in >> poly_name;
|
|
||||||
SG_LOG( SG_CLIPPER, SG_INFO, "poly name = " << poly_name);
|
|
||||||
poly_type = get_area_type( poly_name );
|
|
||||||
SG_LOG( SG_CLIPPER, SG_INFO, "poly type (int) = " << (int)poly_type);
|
|
||||||
in >> contours;
|
|
||||||
SG_LOG( SG_CLIPPER, SG_INFO, "num contours = " << contours);
|
|
||||||
|
|
||||||
poly.erase();
|
|
||||||
|
|
||||||
for ( i = 0; i < contours; ++i ) {
|
|
||||||
in >> count;
|
|
||||||
|
|
||||||
if ( count < 3 ) {
|
|
||||||
SG_LOG( SG_CLIPPER, SG_ALERT, "Polygon with less than 3 data points." );
|
|
||||||
exit(-1);
|
|
||||||
}
|
|
||||||
|
|
||||||
in >> hole_flag;
|
|
||||||
|
|
||||||
in >> startx;
|
|
||||||
in >> starty;
|
|
||||||
OSRef = Point3D(startx, starty, -9999.0);
|
|
||||||
|
|
||||||
//Convert from OSGB36 Eastings/Northings to WGS84 Lat/Lon
|
|
||||||
//Note that startx and starty themselves must not be altered since we compare them with unaltered lastx and lasty later
|
|
||||||
p = OSGB36ToWGS84(OSRef);
|
|
||||||
|
|
||||||
poly.add_node( i, p );
|
|
||||||
nodes.unique_add( p );
|
|
||||||
|
|
||||||
SG_LOG( SG_CLIPPER, SG_BULK, "0 = " << startx << ", " << starty );
|
|
||||||
|
|
||||||
for ( j = 1; j < count - 1; ++j ) {
|
|
||||||
in >> x;
|
|
||||||
in >> y;
|
|
||||||
OSRef = Point3D( x, y, -9999.0 );
|
|
||||||
p = OSGB36ToWGS84(OSRef);
|
|
||||||
|
|
||||||
poly.add_node( i, p );
|
|
||||||
nodes.unique_add( p );
|
|
||||||
SG_LOG( SG_CLIPPER, SG_BULK, j << " = " << x << ", " << y );
|
|
||||||
}
|
|
||||||
|
|
||||||
in >> lastx;
|
|
||||||
in >> lasty;
|
|
||||||
|
|
||||||
if ( (fabs(startx - lastx) < SG_EPSILON) &&
|
|
||||||
(fabs(starty - lasty) < SG_EPSILON) ) {
|
|
||||||
// last point same as first, discard
|
|
||||||
} else {
|
|
||||||
OSRef = Point3D( lastx, lasty, -9999.0 );
|
|
||||||
p = OSGB36ToWGS84(OSRef);
|
|
||||||
|
|
||||||
poly.add_node( i, p );
|
|
||||||
nodes.unique_add( p );
|
|
||||||
SG_LOG( SG_CLIPPER, SG_BULK, count - 1 << " = " << lastx << ", " << lasty );
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// TODO : Make like OGR
|
|
||||||
int area = (int)poly_type;
|
|
||||||
string material = get_area_name( area );
|
|
||||||
add_poly(area, poly, material);
|
|
||||||
// END TODO
|
|
||||||
|
|
||||||
in >> skipcomment;
|
|
||||||
}
|
|
||||||
|
|
||||||
return true;
|
|
||||||
}
|
|
||||||
|
|
||||||
// load all 2d polygons from the specified load disk directories and
|
// load all 2d polygons from the specified load disk directories and
|
||||||
// clip against each other to resolve any overlaps
|
// clip against each other to resolve any overlaps
|
||||||
int TGConstruct::LoadLandclassPolys( void ) {
|
int TGConstruct::LoadLandclassPolys( void ) {
|
||||||
|
|
@ -403,10 +306,6 @@ int TGConstruct::LoadLandclassPolys( void ) {
|
||||||
(lext == "fit") || (lext == "fit.gz") || (lext == "ind"))
|
(lext == "fit") || (lext == "fit.gz") || (lext == "ind"))
|
||||||
{
|
{
|
||||||
// skipped!
|
// skipped!
|
||||||
} else if (lext == "osgb36") {
|
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " Loading osgb36 poly definition file " << p.file());
|
|
||||||
load_osgb36_poly( p.str() );
|
|
||||||
++count;
|
|
||||||
} else {
|
} else {
|
||||||
load_poly( p.str() );
|
load_poly( p.str() );
|
||||||
SG_LOG(SG_GENERAL, SG_ALERT, " Loaded " << p.file());
|
SG_LOG(SG_GENERAL, SG_ALERT, " Loaded " << p.file());
|
||||||
|
|
|
||||||
1
src/BuildTiles/Osgb36/.gitignore
vendored
1
src/BuildTiles/Osgb36/.gitignore
vendored
|
|
@ -1 +0,0 @@
|
||||||
testosgb36
|
|
||||||
|
|
@ -1,8 +0,0 @@
|
||||||
|
|
||||||
|
|
||||||
add_library(Osgb36 STATIC
|
|
||||||
osgb36.cxx osgb36.hxx osgbtc.cxx osgbtc.hxx uk.cxx uk.hxx
|
|
||||||
)
|
|
||||||
|
|
||||||
add_executable(testosgb36 testosgb36.cxx)
|
|
||||||
target_link_libraries(testosgb36 Osgb36)
|
|
||||||
|
|
@ -1,583 +0,0 @@
|
||||||
// osgb36.cxx -- Routines to convert UK OS coordinates to and
|
|
||||||
// from world WGS84 coordinates.
|
|
||||||
//
|
|
||||||
// Written by David Luff, started December 2000.
|
|
||||||
//
|
|
||||||
// Copyright (C) 2000 David C. Luff - david.luff@nottingham.ac.uk
|
|
||||||
//
|
|
||||||
// This program is free software; you can redistribute it and/or
|
|
||||||
// modify it under the terms of the GNU General Public License as
|
|
||||||
// published by the Free Software Foundation; either version 2 of the
|
|
||||||
// License, or (at your option) any later version.
|
|
||||||
//
|
|
||||||
// This program is distributed in the hope that it will be useful, but
|
|
||||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
||||||
// General Public License for more details.
|
|
||||||
//
|
|
||||||
// You should have received a copy of the GNU General Public License
|
|
||||||
// along with this program; if not, write to the Free Software
|
|
||||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
||||||
|
|
||||||
#include <simgear/compiler.h>
|
|
||||||
|
|
||||||
#include <iostream>
|
|
||||||
#include <math.h>
|
|
||||||
|
|
||||||
#include "osgb36.hxx"
|
|
||||||
|
|
||||||
using std::cout;
|
|
||||||
|
|
||||||
double DEG_TO_RAD = 2.0 * 3.14159265358979323846264338327950288419716939967511 / 360.0;
|
|
||||||
double RAD_TO_DEG = 360.0 / (2.0 * 3.14159265358979323846264338327950288419716939967511);
|
|
||||||
|
|
||||||
/********************************************************************
|
|
||||||
*
|
|
||||||
* FORWARD DECLARATIONS
|
|
||||||
*
|
|
||||||
********************************************************************/
|
|
||||||
//*******************************************************************
|
|
||||||
//
|
|
||||||
// NOTE All functions which take or return Latitude and Longitude
|
|
||||||
// take and return degrees. Conversion to and from radians
|
|
||||||
// is performed internaly. XYZ coordinates and ellipsoid
|
|
||||||
// heights are returned and required in meters.
|
|
||||||
//
|
|
||||||
//*******************************************************************
|
|
||||||
|
|
||||||
//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference)
|
|
||||||
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon);
|
|
||||||
|
|
||||||
//Convert OSGB36 Eastings and Northings to OSGB36 Lat and Lon
|
|
||||||
Point3D ConvertEastingsNorthingsToLatLon(Point3D GridRef);
|
|
||||||
|
|
||||||
//Convert OSGB36 coordinates from polar to cartesian (x,y,z) form
|
|
||||||
Point3D ConvertAiry1830PolarToCartesian(Point3D LatLon);
|
|
||||||
|
|
||||||
//Convert WGS84 coordinates from polar to cartesian (x,y,z) form
|
|
||||||
Point3D ConvertGRS80PolarToCartesian(Point3D LatLon);
|
|
||||||
|
|
||||||
//Convert OSGB36 coordinates from cartesian (x,y,z) to polar form
|
|
||||||
Point3D ConvertAiry1830CartesianToPolar(Point3D XYZCoord);
|
|
||||||
|
|
||||||
//Convert WGS84 coordinates from cartesian (x,y,z) to polar form
|
|
||||||
Point3D ConvertGRS80CartesianToPolar(Point3D XYZCoord);
|
|
||||||
|
|
||||||
//Transform a point in WGS84 cartesian coordinates to OSGB36 cartesian coordinates
|
|
||||||
//Uses the Helmert Transformation with coefficients from www.gps.gov.org
|
|
||||||
//Only accurate to around 5m since OSGB36 is an inhomogenous TRF
|
|
||||||
Point3D ConvertWGS84ToOSGB36(Point3D WGS84);
|
|
||||||
|
|
||||||
//Transform a point in OSGB36 cartesian coordinates to WGS84 cartesian coordinates
|
|
||||||
//Uses the Helmert Transformation with coefficients from www.gps.gov.org
|
|
||||||
//Only accurate to around 5m since OSGB36 is an inhomogenous TRF
|
|
||||||
Point3D ConvertOSGB36ToWGS84(Point3D OSGB36);
|
|
||||||
|
|
||||||
/****************************************************************************/
|
|
||||||
|
|
||||||
|
|
||||||
// Start of implementation proper.
|
|
||||||
|
|
||||||
// Convert WGS84 Lat/Lon to OSGB36 Eastings and Northings
|
|
||||||
Point3D WGS84ToOSGB36(Point3D p) {
|
|
||||||
p = ConvertGRS80PolarToCartesian(p);
|
|
||||||
p = ConvertWGS84ToOSGB36(p);
|
|
||||||
p = ConvertAiry1830CartesianToPolar(p);
|
|
||||||
p = ConvertLatLonToEastingsNorthings(p);
|
|
||||||
return(p);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Convert OSGB36 Eastings and Northings to WGS84 Lat/Lon
|
|
||||||
Point3D OSGB36ToWGS84(Point3D p) {
|
|
||||||
p = ConvertEastingsNorthingsToLatLon(p);
|
|
||||||
p = ConvertAiry1830PolarToCartesian(p);
|
|
||||||
p = ConvertOSGB36ToWGS84(p);
|
|
||||||
p = ConvertGRS80CartesianToPolar(p);
|
|
||||||
return(p);
|
|
||||||
}
|
|
||||||
|
|
||||||
//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference)
|
|
||||||
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon)
|
|
||||||
{
|
|
||||||
// ellipsoid constants
|
|
||||||
double a; //semi-major axis (m)
|
|
||||||
double b; //semi-minor axis (m)
|
|
||||||
double e_squared; //ellipsoid squared eccentricity constant
|
|
||||||
|
|
||||||
// Airy 1830
|
|
||||||
a = 6377563.396;
|
|
||||||
b = 6356256.910;
|
|
||||||
|
|
||||||
// GRS80 (aka WGS84)
|
|
||||||
// a = 6378137.000;
|
|
||||||
// b = 6356752.3141;
|
|
||||||
|
|
||||||
e_squared = ((a*a)-(b*b))/(a*a);
|
|
||||||
|
|
||||||
// Projection constants
|
|
||||||
double N_zero; //northing of true origin (m)
|
|
||||||
double E_zero; //easting of true origin (m)
|
|
||||||
double F_zero; //scale factor on central meridian
|
|
||||||
double lat_zero; //latitude of true origin (rad)
|
|
||||||
double lon_zero; //longitude of true origin and central meridian (rad)
|
|
||||||
|
|
||||||
// OSGB36 (UK National Grid)
|
|
||||||
N_zero = -100000;
|
|
||||||
E_zero = 400000;
|
|
||||||
F_zero = 0.9996012717;
|
|
||||||
lat_zero = 49 * DEG_TO_RAD;
|
|
||||||
lon_zero = -2 * DEG_TO_RAD;
|
|
||||||
|
|
||||||
double n;
|
|
||||||
double v;
|
|
||||||
double rho;
|
|
||||||
double neta_squared;
|
|
||||||
double M;
|
|
||||||
double I;
|
|
||||||
double II;
|
|
||||||
double III;
|
|
||||||
double IIIA;
|
|
||||||
double IV;
|
|
||||||
double V;
|
|
||||||
double VI;
|
|
||||||
|
|
||||||
double lon = LatLon.lon() * DEG_TO_RAD;
|
|
||||||
double lat = LatLon.lat() * DEG_TO_RAD;
|
|
||||||
|
|
||||||
n = (a-b)/(a+b);
|
|
||||||
v = a*F_zero*pow( (1-e_squared*sin(lat)*sin(lat)) , -0.5 );
|
|
||||||
rho = a*F_zero*(1-e_squared)*pow((1-e_squared*sin(lat)*sin(lat)) , -1.5 );
|
|
||||||
neta_squared = v/rho - 1.0;
|
|
||||||
|
|
||||||
// terms in the M equation to make it more manageable (eqn A11)
|
|
||||||
double term1;
|
|
||||||
double term2;
|
|
||||||
double term3;
|
|
||||||
double term4;
|
|
||||||
|
|
||||||
term1 = (1.0 + n + (5.0/4.0)*n*n + (5.0/4.0)*n*n*n)*(lat-lat_zero);
|
|
||||||
term2 = ((3.0*n) + (3.0*n*n) + ((21.0/8.0)*n*n*n))*(sin(lat-lat_zero))*(cos(lat+lat_zero));
|
|
||||||
term3 = (((15.0/8.0)*n*n) + ((15.0/8.0)*n*n*n))*(sin(2.0*(lat-lat_zero)))*(cos(2.0*(lat+lat_zero)));
|
|
||||||
term4 = ((35.0/24.0)*n*n*n)*(sin(3.0*(lat-lat_zero)))*(cos(3.0*(lat+lat_zero)));
|
|
||||||
|
|
||||||
M = b * F_zero * (term1 - term2 + term3 - term4); //Eqn A11
|
|
||||||
|
|
||||||
I = M + N_zero;
|
|
||||||
II = (v/2.0)*sin(lat)*cos(lat);
|
|
||||||
III = (v/24.0)*sin(lat)*pow(cos(lat),3)*(5.0 - pow(tan(lat),2) + 9.0*neta_squared);
|
|
||||||
IIIA = (v/720.0)*sin(lat)*pow(cos(lat),5)*(61.0 - 58.0*pow(tan(lat),2) + pow(tan(lat),4));
|
|
||||||
IV = v*cos(lat);
|
|
||||||
V = (v/6.0)*pow(cos(lat),3)*(v/rho - pow(tan(lat),2));
|
|
||||||
VI = (v/120.0)*pow(cos(lat),5)*(5.0 - 18.0*pow(tan(lat),2) + pow(tan(lat),4) + 14.0*neta_squared - 58.0*pow(tan(lat),2)*neta_squared);
|
|
||||||
|
|
||||||
double N; //calculated Northing
|
|
||||||
double E; //calculated Easting
|
|
||||||
|
|
||||||
N = I + II*(lon-lon_zero)*(lon-lon_zero) + III*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero) + IIIA*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero); //eqn A12
|
|
||||||
E = E_zero + IV*(lon-lon_zero) + V*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero) + VI*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero)*(lon-lon_zero); //eqn A13
|
|
||||||
|
|
||||||
Point3D OSref;
|
|
||||||
|
|
||||||
OSref.setx(E);
|
|
||||||
OSref.sety(N);
|
|
||||||
OSref.setz(0.0); // Not used at present but could be used to return ODN height in future
|
|
||||||
|
|
||||||
return(OSref);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
Point3D ConvertEastingsNorthingsToLatLon(Point3D GridRef)
|
|
||||||
{
|
|
||||||
// ellipsoid constants
|
|
||||||
double a; //semi-major axis (m)
|
|
||||||
double b; //semi-minor axis (m)
|
|
||||||
double e_squared; //ellipsoid squared eccentricity constant
|
|
||||||
|
|
||||||
// Airy 1830
|
|
||||||
a = 6377563.396;
|
|
||||||
b = 6356256.910;
|
|
||||||
|
|
||||||
// GRS80 (aka WGS84)
|
|
||||||
// a = 6378137.000;
|
|
||||||
// b = 6356752.3141;
|
|
||||||
|
|
||||||
e_squared = ((a*a)-(b*b))/(a*a);
|
|
||||||
|
|
||||||
// Projection constants
|
|
||||||
double N_zero; //northing of true origin (m)
|
|
||||||
double E_zero; //easting of true origin (m)
|
|
||||||
double F_zero; //scale factor on central meridian
|
|
||||||
double lat_zero; //latitude of true origin (rad)
|
|
||||||
double lon_zero; //longitude of true origin and central meridian (rad)
|
|
||||||
|
|
||||||
// OSGB36 (UK National Grid)
|
|
||||||
N_zero = -100000;
|
|
||||||
E_zero = 400000;
|
|
||||||
F_zero = 0.9996012717;
|
|
||||||
lat_zero = 49 * DEG_TO_RAD;
|
|
||||||
lon_zero = -2 * DEG_TO_RAD;
|
|
||||||
|
|
||||||
double n;
|
|
||||||
double M;
|
|
||||||
|
|
||||||
n = (a-b)/(a+b);
|
|
||||||
|
|
||||||
double E = GridRef.x();
|
|
||||||
double N = GridRef.y();
|
|
||||||
|
|
||||||
// terms in the M equation to make it more manageable (eqn A11)
|
|
||||||
double term1;
|
|
||||||
double term2;
|
|
||||||
double term3;
|
|
||||||
double term4;
|
|
||||||
|
|
||||||
double lat_prime;
|
|
||||||
|
|
||||||
//calculate initial lat_prime
|
|
||||||
lat_prime = ((N-N_zero)/(a*F_zero)) + lat_zero;
|
|
||||||
|
|
||||||
int i;
|
|
||||||
|
|
||||||
for(i=0;i<100;i++)
|
|
||||||
{
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
term1 = (1.0 + n + ((5.0/4.0)*n*n) + ((5.0/4.0)*n*n*n))*(lat_prime-lat_zero);
|
|
||||||
term2 = ((3.0*n) + (3.0*n*n) + ((21.0/8.0)*n*n*n))*(sin(lat_prime-lat_zero))*(cos(lat_prime+lat_zero));
|
|
||||||
term3 = (((15.0/8.0)*n*n) + ((15.0/8.0)*n*n*n))*(sin(2.0*(lat_prime-lat_zero)))*(cos(2.0*(lat_prime+lat_zero)));
|
|
||||||
term4 = ((35.0/24.0)*n*n*n)*(sin(3.0*(lat_prime-lat_zero)))*(cos(3.0*(lat_prime+lat_zero)));
|
|
||||||
|
|
||||||
M = b * F_zero * (term1 - term2 + term3 - term4); //Eqn A11
|
|
||||||
|
|
||||||
if((N - N_zero - M) < 0.001)
|
|
||||||
break; //N - N_zero - M is less than 1mm
|
|
||||||
else //recompute lat_prime and keep iterating until it is
|
|
||||||
lat_prime += ((N - N_zero - M)/(a*F_zero));
|
|
||||||
|
|
||||||
if(i == 99)
|
|
||||||
cout << "WARNING: Iteration failed to converge in ConvertEastingsNorthingsToLatLon\n";
|
|
||||||
}
|
|
||||||
|
|
||||||
double v;
|
|
||||||
double rho;
|
|
||||||
double neta_squared;
|
|
||||||
double VII;
|
|
||||||
double VIII;
|
|
||||||
double IX;
|
|
||||||
double X;
|
|
||||||
double XI;
|
|
||||||
double XII;
|
|
||||||
double XIIA;
|
|
||||||
|
|
||||||
v = a*F_zero*pow( (1-e_squared*sin(lat_prime)*sin(lat_prime)) , -0.5 );
|
|
||||||
rho = a*F_zero*(1-e_squared)*pow((1-e_squared*sin(lat_prime)*sin(lat_prime)) , -1.5 );
|
|
||||||
neta_squared = v/rho - 1.0;
|
|
||||||
|
|
||||||
VII = tan(lat_prime) / (2.0*rho*v);
|
|
||||||
VIII = (tan(lat_prime)/(24.0*rho*v*v*v)) * (5.0 + 3.0*pow(tan(lat_prime),2.0) + neta_squared - 9.0*pow(tan(lat_prime),2.0)*neta_squared);
|
|
||||||
IX = (tan(lat_prime)/(720.0*rho*v*v*v*v*v)) * (61.0 + 90.0*pow(tan(lat_prime),2.0) + 45.0*pow(tan(lat_prime),4.0));
|
|
||||||
X = (1/cos(lat_prime)) / v;
|
|
||||||
XI = ((1/cos(lat_prime)) / (6.0 * pow(v,3))) * (v/rho + 2*pow(tan(lat_prime),2));
|
|
||||||
XII = ((1/cos(lat_prime)) / (120.0 * pow(v,5))) * (5.0 + 28.0*pow(tan(lat_prime),2) + 24.0*pow(tan(lat_prime),4));
|
|
||||||
XIIA = ((1/cos(lat_prime)) / (5040.0 * pow(v,7))) * (61.0 + 662.0*pow(tan(lat_prime),2) + 1320.0*pow(tan(lat_prime),4) + 720.0*pow(tan(lat_prime),6));
|
|
||||||
|
|
||||||
double lat;
|
|
||||||
double lon;
|
|
||||||
|
|
||||||
lat = lat_prime - VII*pow((E-E_zero),2) + VIII*pow((E-E_zero),4) + IX*pow((E-E_zero),6);
|
|
||||||
lon = lon_zero + X*(E-E_zero) - XI*pow((E-E_zero),3) + XII*pow((E-E_zero),5) - XIIA*pow((E-E_zero),7);
|
|
||||||
|
|
||||||
Point3D LatLon;
|
|
||||||
LatLon.setlon(lon * RAD_TO_DEG);
|
|
||||||
LatLon.setlat(lat * RAD_TO_DEG);
|
|
||||||
LatLon.setelev(0.0);
|
|
||||||
|
|
||||||
return(LatLon);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
Point3D ConvertAiry1830PolarToCartesian(Point3D LatLon)
|
|
||||||
{
|
|
||||||
|
|
||||||
double lon = LatLon.lon() * DEG_TO_RAD;
|
|
||||||
double lat = LatLon.lat() * DEG_TO_RAD;
|
|
||||||
double H = LatLon.elev() * DEG_TO_RAD;
|
|
||||||
|
|
||||||
// ellipsoid constants
|
|
||||||
double a; //semi-major axis (m)
|
|
||||||
double b; //semi-minor axis (m)
|
|
||||||
double e_squared; //ellipsoid squared eccentricity constant
|
|
||||||
|
|
||||||
// Airy 1830
|
|
||||||
a = 6377563.396;
|
|
||||||
b = 6356256.910;
|
|
||||||
|
|
||||||
e_squared = ((a*a)-(b*b))/(a*a);
|
|
||||||
|
|
||||||
double v;
|
|
||||||
|
|
||||||
v = a*pow( (1-e_squared*sin(lat)*sin(lat)) , -0.5 );
|
|
||||||
|
|
||||||
double x;
|
|
||||||
double y;
|
|
||||||
double z;
|
|
||||||
Point3D Cartesian;
|
|
||||||
|
|
||||||
x = (v+H)*cos(lat)*cos(lon);
|
|
||||||
y = (v+H)*cos(lat)*sin(lon);
|
|
||||||
z = ((1-e_squared)*v + H)*sin(lat);
|
|
||||||
|
|
||||||
Cartesian.setx(x);
|
|
||||||
Cartesian.sety(y);
|
|
||||||
Cartesian.setz(z);
|
|
||||||
|
|
||||||
return(Cartesian);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
Point3D ConvertGRS80PolarToCartesian(Point3D LatLon)
|
|
||||||
{
|
|
||||||
|
|
||||||
double lon = LatLon.lon() * DEG_TO_RAD;
|
|
||||||
double lat = LatLon.lat() * DEG_TO_RAD;
|
|
||||||
double H = LatLon.elev() * DEG_TO_RAD;
|
|
||||||
|
|
||||||
// ellipsoid constants
|
|
||||||
double a; //semi-major axis (m)
|
|
||||||
double b; //semi-minor axis (m)
|
|
||||||
double e_squared; //ellipsoid squared eccentricity constant
|
|
||||||
|
|
||||||
|
|
||||||
// GRS80 (aka WGS84)
|
|
||||||
a = 6378137.000;
|
|
||||||
b = 6356752.3141;
|
|
||||||
|
|
||||||
e_squared = ((a*a)-(b*b))/(a*a);
|
|
||||||
|
|
||||||
double v;
|
|
||||||
|
|
||||||
v = a*pow( (1-e_squared*sin(lat)*sin(lat)) , -0.5 );
|
|
||||||
|
|
||||||
double x;
|
|
||||||
double y;
|
|
||||||
double z;
|
|
||||||
Point3D Cartesian;
|
|
||||||
|
|
||||||
x = (v+H)*cos(lat)*cos(lon);
|
|
||||||
y = (v+H)*cos(lat)*sin(lon);
|
|
||||||
z = ((1-e_squared)*v + H)*sin(lat);
|
|
||||||
|
|
||||||
Cartesian.setx(x);
|
|
||||||
Cartesian.sety(y);
|
|
||||||
Cartesian.setz(z);
|
|
||||||
|
|
||||||
return(Cartesian);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
Point3D ConvertAiry1830CartesianToPolar(Point3D XYZCoord)
|
|
||||||
{
|
|
||||||
double x = XYZCoord.x();
|
|
||||||
double y = XYZCoord.y();
|
|
||||||
double z = XYZCoord.z();
|
|
||||||
double p;
|
|
||||||
double lat1;
|
|
||||||
double lat2;
|
|
||||||
double lat;
|
|
||||||
double lon;
|
|
||||||
double H;
|
|
||||||
Point3D LatLon;
|
|
||||||
|
|
||||||
// ellipsoid constants
|
|
||||||
double a; //semi-major axis (m)
|
|
||||||
double b; //semi-minor axis (m)
|
|
||||||
double e_squared; //ellipsoid squared eccentricity constant
|
|
||||||
|
|
||||||
// Airy 1830
|
|
||||||
a = 6377563.396;
|
|
||||||
b = 6356256.910;
|
|
||||||
|
|
||||||
e_squared = ((a*a)-(b*b))/(a*a);
|
|
||||||
|
|
||||||
double v;
|
|
||||||
|
|
||||||
lon = atan(y/x);
|
|
||||||
|
|
||||||
p = sqrt(x*x + y*y);
|
|
||||||
|
|
||||||
//lat is obtained iteratively
|
|
||||||
//obtain initial value of lat
|
|
||||||
lat1 = atan(z/(p*(1-e_squared)));
|
|
||||||
|
|
||||||
// cout << "Initial value of lat = " << lat1 << '\n';
|
|
||||||
|
|
||||||
int i;
|
|
||||||
|
|
||||||
for(i=0;i<100;i++)
|
|
||||||
{
|
|
||||||
v = a*pow( (1-e_squared*sin(lat1)*sin(lat1)) , -0.5 );
|
|
||||||
|
|
||||||
lat2 = atan((z + e_squared*v*sin(lat1))/p);
|
|
||||||
|
|
||||||
// cout << "lat2 = " << lat2 << '\n';
|
|
||||||
|
|
||||||
if(fabs(lat2-lat1) < 0.00000001)
|
|
||||||
{
|
|
||||||
lat = lat2;
|
|
||||||
// cout << i << " iterations required in ConvertAiry1830CartesianToPolar\n";
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
lat1 = lat2;
|
|
||||||
}
|
|
||||||
|
|
||||||
if(i == 99)
|
|
||||||
cout << "WARNING: Iteration failed to converge in ConvertAiry1830CartesianToPolar\n";
|
|
||||||
}
|
|
||||||
|
|
||||||
H = p/cos(lat) - v;
|
|
||||||
|
|
||||||
LatLon.setlon(lon * RAD_TO_DEG);
|
|
||||||
LatLon.setlat(lat * RAD_TO_DEG);
|
|
||||||
LatLon.setelev(H);
|
|
||||||
|
|
||||||
return(LatLon);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
Point3D ConvertGRS80CartesianToPolar(Point3D XYZCoord)
|
|
||||||
{
|
|
||||||
double x = XYZCoord.x();
|
|
||||||
double y = XYZCoord.y();
|
|
||||||
double z = XYZCoord.z();
|
|
||||||
double p;
|
|
||||||
double lat1;
|
|
||||||
double lat2;
|
|
||||||
double lat;
|
|
||||||
double lon;
|
|
||||||
double H;
|
|
||||||
Point3D LatLon;
|
|
||||||
|
|
||||||
// ellipsoid constants
|
|
||||||
double a; //semi-major axis (m)
|
|
||||||
double b; //semi-minor axis (m)
|
|
||||||
double e_squared; //ellipsoid squared eccentricity constant
|
|
||||||
|
|
||||||
// GRS80 (aka WGS84)
|
|
||||||
a = 6378137.000;
|
|
||||||
b = 6356752.3141;
|
|
||||||
|
|
||||||
e_squared = ((a*a)-(b*b))/(a*a);
|
|
||||||
|
|
||||||
double v;
|
|
||||||
|
|
||||||
lon = atan(y/x);
|
|
||||||
|
|
||||||
p = sqrt(x*x + y*y);
|
|
||||||
|
|
||||||
//lat is obtained iteratively
|
|
||||||
//obtain initial value of lat
|
|
||||||
lat1 = atan(z/(p*(1-e_squared)));
|
|
||||||
|
|
||||||
int i;
|
|
||||||
|
|
||||||
for(i=0;i<100;i++)
|
|
||||||
{
|
|
||||||
v = a*pow( (1-e_squared*sin(lat1)*sin(lat1)) , -0.5 );
|
|
||||||
|
|
||||||
lat2 = atan((z + e_squared*v*sin(lat1))/p);
|
|
||||||
|
|
||||||
if(fabs(lat2-lat1) < 0.00000001)
|
|
||||||
{
|
|
||||||
lat = lat2;
|
|
||||||
// cout << i << " iterations required in ConvertGRS80CartesianToPolar\n";
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
lat1 = lat2;
|
|
||||||
}
|
|
||||||
|
|
||||||
if(i == 99)
|
|
||||||
cout << "WARNING: Iteration failed to converge in ConvertGRS80CartesianToPolar\n";
|
|
||||||
}
|
|
||||||
|
|
||||||
H = p/cos(lat) - v;
|
|
||||||
|
|
||||||
LatLon.setlon(lon * RAD_TO_DEG);
|
|
||||||
LatLon.setlat(lat * RAD_TO_DEG);
|
|
||||||
LatLon.setelev(H);
|
|
||||||
|
|
||||||
return(LatLon);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
//Transform a point in WGS84 cartesian coordinates to OSGB36 cartesian coordinates
|
|
||||||
//Uses the Helmert Transformation with coefficients from www.gps.gov.org
|
|
||||||
//Only accurate to around 5m since OSGB36 is an imhomogenous TRF
|
|
||||||
Point3D ConvertWGS84ToOSGB36(Point3D WGS84)
|
|
||||||
{
|
|
||||||
Point3D OSGB36;
|
|
||||||
|
|
||||||
double tx = -446.448; // (m)
|
|
||||||
double ty = 125.157;
|
|
||||||
double tz = -542.060;
|
|
||||||
double rx = -0.1502; // (sec)
|
|
||||||
double ry = -0.2470;
|
|
||||||
double rz = -0.8421;
|
|
||||||
double s = 20.4894; // scale factor in ppm
|
|
||||||
|
|
||||||
// Convert rotations from seconds of arc to radians
|
|
||||||
rx = rx / 3600.0 * DEG_TO_RAD;
|
|
||||||
ry = ry / 3600.0 * DEG_TO_RAD;
|
|
||||||
rz = rz / 3600.0 * DEG_TO_RAD;
|
|
||||||
|
|
||||||
s /= 1000000;
|
|
||||||
|
|
||||||
// cout << "(1+s) = " << (1+s) << '\n';
|
|
||||||
|
|
||||||
OSGB36.setx(tx + (1+s)*WGS84.x() - rz*WGS84.y() + ry*WGS84.z());
|
|
||||||
OSGB36.sety(ty + rz*WGS84.x() + (1+s)*WGS84.y() - rx*WGS84.z());
|
|
||||||
OSGB36.setz(tz - ry*WGS84.x() + rx*WGS84.y() + (1+s)*WGS84.z());
|
|
||||||
|
|
||||||
return(OSGB36);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
//Reverse transformation achieved by reversing the signs of coefficients for the above.
|
|
||||||
//Only valid for small angles, and has same accuracy limitations
|
|
||||||
Point3D ConvertOSGB36ToWGS84(Point3D OSGB36)
|
|
||||||
{
|
|
||||||
Point3D WGS84;
|
|
||||||
|
|
||||||
double tx = 446.448; // (m)
|
|
||||||
double ty = -125.157;
|
|
||||||
double tz = 542.060;
|
|
||||||
double rx = 0.1502; // (sec)
|
|
||||||
double ry = 0.2470;
|
|
||||||
double rz = 0.8421;
|
|
||||||
double s = -20.4894; // scale factor in ppm
|
|
||||||
|
|
||||||
// Convert rotations from seconds to radians
|
|
||||||
rx = rx / 3600.0 * DEG_TO_RAD;
|
|
||||||
ry = ry / 3600.0 * DEG_TO_RAD;
|
|
||||||
rz = rz / 3600.0 * DEG_TO_RAD;
|
|
||||||
|
|
||||||
s /= 1000000;
|
|
||||||
|
|
||||||
WGS84.setx(tx + (1+s)*OSGB36.x() - rz*OSGB36.y() + ry*OSGB36.z());
|
|
||||||
WGS84.sety(ty + rz*OSGB36.x() + (1+s)*OSGB36.y() - rx*OSGB36.z());
|
|
||||||
WGS84.setz(tz - ry*OSGB36.x() + rx*OSGB36.y() + (1+s)*OSGB36.z());
|
|
||||||
|
|
||||||
return(WGS84);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -1,36 +0,0 @@
|
||||||
// osgb36.hxx -- Routines to convert UK OS coordinates to and
|
|
||||||
// from world WGS84 coordinates.
|
|
||||||
//
|
|
||||||
// Written by David Luff, started December 2000.
|
|
||||||
//
|
|
||||||
// Copyright (C) 2000 David C. Luff - david.luff@nottingham.ac.uk
|
|
||||||
//
|
|
||||||
// This program is free software; you can redistribute it and/or
|
|
||||||
// modify it under the terms of the GNU General Public License as
|
|
||||||
// published by the Free Software Foundation; either version 2 of the
|
|
||||||
// License, or (at your option) any later version.
|
|
||||||
//
|
|
||||||
// This program is distributed in the hope that it will be useful, but
|
|
||||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
||||||
// General Public License for more details.
|
|
||||||
//
|
|
||||||
// You should have received a copy of the GNU General Public License
|
|
||||||
// along with this program; if not, write to the Free Software
|
|
||||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
||||||
|
|
||||||
#include <Geometry/point3d.hxx>
|
|
||||||
|
|
||||||
//*******************************************************************
|
|
||||||
//
|
|
||||||
// NOTE All functions which take or return Latitude and Longitude
|
|
||||||
// take and return degrees. XYZ coordinates and ellipsoid
|
|
||||||
// heights are returned and required in meters.
|
|
||||||
//
|
|
||||||
//*******************************************************************
|
|
||||||
|
|
||||||
// Convert WGS84 Lat/Lon (degrees) to OSGB36 Eastings and Northings (meters)
|
|
||||||
Point3D WGS84ToOSGB36(Point3D p);
|
|
||||||
|
|
||||||
// Convert OSGB36 Eastings and Northings (meters) to WGS84 Lat/Lon (degrees)
|
|
||||||
Point3D OSGB36ToWGS84(Point3D p);
|
|
||||||
|
|
@ -1,168 +0,0 @@
|
||||||
#ifdef HAVE_CONFIG_H
|
|
||||||
# include <config.h>
|
|
||||||
#endif
|
|
||||||
|
|
||||||
#include "osgb36.hxx"
|
|
||||||
#include "uk.hxx"
|
|
||||||
|
|
||||||
#include "osgbtc.hxx"
|
|
||||||
|
|
||||||
|
|
||||||
// traverse the specified fan/strip/list of vertices and attempt to
|
|
||||||
// calculate "none stretching" texture coordinates
|
|
||||||
point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
|
|
||||||
const int_list& fan, double scale )
|
|
||||||
{
|
|
||||||
// find min/max of fan
|
|
||||||
Point3D t, tmin, tmax;
|
|
||||||
Point3D OSGridRef;
|
|
||||||
bool first = true;
|
|
||||||
|
|
||||||
int i;
|
|
||||||
|
|
||||||
for ( i = 0; i < (int)fan.size(); ++i ) {
|
|
||||||
|
|
||||||
OSGridRef = WGS84ToOSGB36(geod_nodes[ fan[i] ]);
|
|
||||||
|
|
||||||
t = UK_basic_tex_coord( OSGridRef );
|
|
||||||
// cout << "basic_tex_coord = " << t << endl;
|
|
||||||
|
|
||||||
if ( first ) {
|
|
||||||
tmin = tmax = t;
|
|
||||||
first = false;
|
|
||||||
} else {
|
|
||||||
if ( t.x() < tmin.x() ) {
|
|
||||||
tmin.setx( t.x() );
|
|
||||||
}
|
|
||||||
if ( t.y() < tmin.y() ) {
|
|
||||||
tmin.sety( t.y() );
|
|
||||||
}
|
|
||||||
if ( t.x() > tmax.x() ) {
|
|
||||||
tmax.setx( t.x() );
|
|
||||||
}
|
|
||||||
if ( t.y() > tmax.y() ) {
|
|
||||||
tmax.sety( t.y() );
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// cout << "tmin.x = " << tmin.x() << " tmax.x = " << tmax.x() << " tmin.y = " << tmin.y() << " tmax.y = " << tmax.y() << '\n';
|
|
||||||
}
|
|
||||||
|
|
||||||
double dx = fabs( tmax.x() - tmin.x() );
|
|
||||||
double dy = fabs( tmax.y() - tmin.y() );
|
|
||||||
// cout << "dx = " << dx << " dy = " << dy << endl;
|
|
||||||
|
|
||||||
bool do_shift = false;
|
|
||||||
// Point3D mod_shift;
|
|
||||||
if ( (dx > HALF_MAX_TEX_COORD) || (dy > HALF_MAX_TEX_COORD) ) {
|
|
||||||
// structure is too big, we'll just have to shift it so that
|
|
||||||
// tmin = (0,0). This messes up subsequent texture scaling,
|
|
||||||
// but is the best we can do.
|
|
||||||
// cout << "SHIFTING" << endl;
|
|
||||||
do_shift = true;
|
|
||||||
if ( tmin.x() < 0 ) {
|
|
||||||
tmin.setx( (double)( (int)tmin.x() - 1 ) );
|
|
||||||
} else {
|
|
||||||
tmin.setx( (int)tmin.x() );
|
|
||||||
}
|
|
||||||
if ( tmin.y() < 0 ) {
|
|
||||||
tmin.sety( (double)( (int)tmin.y() - 1 ) );
|
|
||||||
} else {
|
|
||||||
tmin.sety( (int)tmin.y() );
|
|
||||||
}
|
|
||||||
// cout << "found tmin = " << tmin << endl;
|
|
||||||
} else {
|
|
||||||
if ( tmin.x() < 0 ) {
|
|
||||||
tmin.setx( ( (int)(tmin.x() / HALF_MAX_TEX_COORD) - 1 )
|
|
||||||
* HALF_MAX_TEX_COORD );
|
|
||||||
} else {
|
|
||||||
tmin.setx( ( (int)(tmin.x() / HALF_MAX_TEX_COORD) )
|
|
||||||
* HALF_MAX_TEX_COORD );
|
|
||||||
}
|
|
||||||
if ( tmin.y() < 0 ) {
|
|
||||||
tmin.sety( ( (int)(tmin.y() / HALF_MAX_TEX_COORD) - 1 )
|
|
||||||
* HALF_MAX_TEX_COORD );
|
|
||||||
} else {
|
|
||||||
tmin.sety( ( (int)(tmin.y() / HALF_MAX_TEX_COORD) )
|
|
||||||
* HALF_MAX_TEX_COORD );
|
|
||||||
}
|
|
||||||
#if 0
|
|
||||||
// structure is small enough ... we can mod it so we can
|
|
||||||
// properly scale the texture coordinates later.
|
|
||||||
// cout << "MODDING" << endl;
|
|
||||||
double x1 = fmod(tmin.x(), MAX_TEX_COORD);
|
|
||||||
while ( x1 < 0 ) { x1 += MAX_TEX_COORD; }
|
|
||||||
|
|
||||||
double y1 = fmod(tmin.y(), MAX_TEX_COORD);
|
|
||||||
while ( y1 < 0 ) { y1 += MAX_TEX_COORD; }
|
|
||||||
|
|
||||||
double x2 = fmod(tmax.x(), MAX_TEX_COORD);
|
|
||||||
while ( x2 < 0 ) { x2 += MAX_TEX_COORD; }
|
|
||||||
|
|
||||||
double y2 = fmod(tmax.y(), MAX_TEX_COORD);
|
|
||||||
while ( y2 < 0 ) { y2 += MAX_TEX_COORD; }
|
|
||||||
|
|
||||||
// At this point we know that the object is < 16 wide in
|
|
||||||
// texture coordinate space. If the modulo of the tmin is >
|
|
||||||
// the mod of the tmax at this point, then we know that the
|
|
||||||
// starting tex coordinate for the tmax > 16 so we can shift
|
|
||||||
// everything down by 16 and get it within the 0-32 range.
|
|
||||||
|
|
||||||
if ( x1 > x2 ) {
|
|
||||||
mod_shift.setx( HALF_MAX_TEX_COORD );
|
|
||||||
} else {
|
|
||||||
mod_shift.setx( 0.0 );
|
|
||||||
}
|
|
||||||
|
|
||||||
if ( y1 > y2 ) {
|
|
||||||
mod_shift.sety( HALF_MAX_TEX_COORD );
|
|
||||||
} else {
|
|
||||||
mod_shift.sety( 0.0 );
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
// cout << "mod_shift = " << mod_shift << endl;
|
|
||||||
}
|
|
||||||
|
|
||||||
// generate tex_list
|
|
||||||
Point3D adjusted_t;
|
|
||||||
point_list tex;
|
|
||||||
tex.clear();
|
|
||||||
for ( i = 0; i < (int)fan.size(); ++i ) {
|
|
||||||
|
|
||||||
OSGridRef = WGS84ToOSGB36(geod_nodes[ fan[i] ]);
|
|
||||||
|
|
||||||
t = UK_basic_tex_coord(OSGridRef);
|
|
||||||
// cout << "second t = " << t << endl;
|
|
||||||
|
|
||||||
// if ( do_shift ) {
|
|
||||||
adjusted_t = t - tmin;
|
|
||||||
// cout << "adjusted_t = " << adjusted_t << '\n';
|
|
||||||
#if 0
|
|
||||||
} else {
|
|
||||||
adjusted_t.setx( fmod(t.x() + mod_shift.x(), MAX_TEX_COORD) );
|
|
||||||
while ( adjusted_t.x() < 0 ) {
|
|
||||||
adjusted_t.setx( adjusted_t.x() + MAX_TEX_COORD );
|
|
||||||
}
|
|
||||||
adjusted_t.sety( fmod(t.y() + mod_shift.y(), MAX_TEX_COORD) );
|
|
||||||
while ( adjusted_t.y() < 0 ) {
|
|
||||||
adjusted_t.sety( adjusted_t.y() + MAX_TEX_COORD );
|
|
||||||
}
|
|
||||||
// cout << "adjusted_t " << adjusted_t << endl;
|
|
||||||
}
|
|
||||||
#endif
|
|
||||||
if ( adjusted_t.x() < SG_EPSILON ) {
|
|
||||||
adjusted_t.setx( 0.0 );
|
|
||||||
}
|
|
||||||
if ( adjusted_t.y() < SG_EPSILON ) {
|
|
||||||
adjusted_t.sety( 0.0 );
|
|
||||||
}
|
|
||||||
adjusted_t.setz( 0.0 );
|
|
||||||
//cout << "adjusted_t after check for SG_EPSILON = " << adjusted_t << endl;
|
|
||||||
|
|
||||||
tex.push_back( adjusted_t );
|
|
||||||
}
|
|
||||||
// char dcl_pause2;
|
|
||||||
// cin >> dcl_pause2;
|
|
||||||
|
|
||||||
return tex;
|
|
||||||
}
|
|
||||||
|
|
@ -1,31 +0,0 @@
|
||||||
#ifndef OSBTEXCOORD_H
|
|
||||||
#define OSBTEXCOORD_H
|
|
||||||
|
|
||||||
|
|
||||||
#include <simgear/bucket/newbucket.hxx>
|
|
||||||
#include <Geometry/point3d.hxx>
|
|
||||||
#include <simgear/math/sg_types.hxx>
|
|
||||||
|
|
||||||
|
|
||||||
#define FG_STANDARD_TEXTURE_DIMENSION 1000.0 // meters
|
|
||||||
#define MAX_TEX_COORD 8.0
|
|
||||||
#define HALF_MAX_TEX_COORD ( MAX_TEX_COORD / 2.0 )
|
|
||||||
|
|
||||||
|
|
||||||
// return the basic unshifted/unmoded texture coordinate for a Easting/Northing (UK) grid reference
|
|
||||||
inline Point3D UK_basic_tex_coord( const Point3D& p )
|
|
||||||
{
|
|
||||||
// cout << "Point in dcl_basic_tex_coord is " << p << '\n';
|
|
||||||
return Point3D( p.x() / FG_STANDARD_TEXTURE_DIMENSION ,
|
|
||||||
p.y() / FG_STANDARD_TEXTURE_DIMENSION ,
|
|
||||||
0.0 );
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
// traverse the specified fan/strip/list of vertices and attempt to
|
|
||||||
// calculate "none stretching" texture coordinates
|
|
||||||
point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
|
|
||||||
const int_list& fan, double scale );
|
|
||||||
|
|
||||||
#endif // OSBTEXCOORD_H
|
|
||||||
|
|
@ -1,39 +0,0 @@
|
||||||
#include <simgear/compiler.h>
|
|
||||||
#include <Geometry/point3d.hxx>
|
|
||||||
|
|
||||||
#include <iostream>
|
|
||||||
#include <stdlib.h>
|
|
||||||
|
|
||||||
#include "osgb36.hxx"
|
|
||||||
|
|
||||||
using std::cout;
|
|
||||||
|
|
||||||
static void Usage() {
|
|
||||||
cout << "Usage is testosgb36 <lon> <lat>\n";
|
|
||||||
exit(-1);
|
|
||||||
}
|
|
||||||
|
|
||||||
//Test harness for WGS84 to OSGB36 conversion.
|
|
||||||
int main(int argc, char** argv)
|
|
||||||
{
|
|
||||||
if(argc != 3) Usage();
|
|
||||||
|
|
||||||
cout << "Test OSGB36\n";
|
|
||||||
|
|
||||||
Point3D p1;
|
|
||||||
p1.setlon(atof(argv[1]));
|
|
||||||
p1.setlat(atof(argv[2]));
|
|
||||||
p1.setelev(0.0);
|
|
||||||
|
|
||||||
cout << "WGS84 position is " << p1 << '\n';
|
|
||||||
|
|
||||||
Point3D p2 = WGS84ToOSGB36(p1);
|
|
||||||
|
|
||||||
cout << "OS coords are " << p2 << '\n';
|
|
||||||
|
|
||||||
p1 = OSGB36ToWGS84(p2);
|
|
||||||
|
|
||||||
cout << "Conversion back to WGS84 gives " << p1 << '\n';
|
|
||||||
|
|
||||||
return(0);
|
|
||||||
}
|
|
||||||
|
|
@ -1,61 +0,0 @@
|
||||||
// uk.cxx -- Routines to determine whether a point is within the UK.
|
|
||||||
//
|
|
||||||
// Written by David Luff, started Janurary 2001.
|
|
||||||
//
|
|
||||||
// Copyright (C) 2001 David C. Luff - david.luff@nottingham.ac.uk
|
|
||||||
//
|
|
||||||
// This program is free software; you can redistribute it and/or
|
|
||||||
// modify it under the terms of the GNU General Public License as
|
|
||||||
// published by the Free Software Foundation; either version 2 of the
|
|
||||||
// License, or (at your option) any later version.
|
|
||||||
//
|
|
||||||
// This program is distributed in the hope that it will be useful, but
|
|
||||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
||||||
// General Public License for more details.
|
|
||||||
//
|
|
||||||
// You should have received a copy of the GNU General Public License
|
|
||||||
// along with this program; if not, write to the Free Software
|
|
||||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
||||||
|
|
||||||
#include "uk.hxx"
|
|
||||||
|
|
||||||
//Returns true if a point is within the mainland UK, excluding
|
|
||||||
//Northern Ireland. Requires lat and lon to be passed in degrees.
|
|
||||||
bool isInUK(Point3D p)
|
|
||||||
{
|
|
||||||
bool inUK = false;
|
|
||||||
double lat = p.lat();
|
|
||||||
double lon = p.lon();
|
|
||||||
//The entire UK (excepting a small portion of Northern Ireland which we arn't doing anyway)
|
|
||||||
//falls within the box of -8 -> 2 deg lon, and 49 -> 60 deg lat
|
|
||||||
//We'll check for within this box first, and if so then we'll refine
|
|
||||||
//further to avoid bits of France and Ireland.
|
|
||||||
if((lat >= 49.0) && (lat <= 60.0) && (lon >= -8.0) && (lon <= 2.0))
|
|
||||||
{
|
|
||||||
//we might be within the UK
|
|
||||||
//set inUK = true and then test for the exceptions
|
|
||||||
inUK = true;
|
|
||||||
|
|
||||||
//check for France (Normandy/Calais)
|
|
||||||
if((lon >= -2.0) && (lat <=50.0))
|
|
||||||
//Normandy
|
|
||||||
inUK = false;
|
|
||||||
if((lon > 1.0) && (lat <= 51.0))
|
|
||||||
//Calais area
|
|
||||||
inUK = false;
|
|
||||||
|
|
||||||
//Check for Ireland
|
|
||||||
if((lon <= -6.0) && (lat > 51.0) && (lat <= 55.5))
|
|
||||||
//Ireland
|
|
||||||
inUK = false;
|
|
||||||
//Check for last bit of NI
|
|
||||||
if((lat > 54.0) && (lat <= 55.1) && (lon <= -5.3333))
|
|
||||||
inUK = false;
|
|
||||||
|
|
||||||
}
|
|
||||||
else
|
|
||||||
inUK = false;
|
|
||||||
|
|
||||||
return(inUK);
|
|
||||||
}
|
|
||||||
|
|
@ -1,25 +0,0 @@
|
||||||
// uk.hxx -- Routines to determine whether a point is within the UK.
|
|
||||||
//
|
|
||||||
// Written by David Luff, started Janurary 2001.
|
|
||||||
//
|
|
||||||
// Copyright (C) 2001 David C. Luff - david.luff@nottingham.ac.uk
|
|
||||||
//
|
|
||||||
// This program is free software; you can redistribute it and/or
|
|
||||||
// modify it under the terms of the GNU General Public License as
|
|
||||||
// published by the Free Software Foundation; either version 2 of the
|
|
||||||
// License, or (at your option) any later version.
|
|
||||||
//
|
|
||||||
// This program is distributed in the hope that it will be useful, but
|
|
||||||
// WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
||||||
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
||||||
// General Public License for more details.
|
|
||||||
//
|
|
||||||
// You should have received a copy of the GNU General Public License
|
|
||||||
// along with this program; if not, write to the Free Software
|
|
||||||
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
||||||
|
|
||||||
#include <Geometry/point3d.hxx>
|
|
||||||
|
|
||||||
//Returns true if a point is within the mainland UK, excluding
|
|
||||||
//Northern Ireland. Requires lat and lon to be passed in degrees.
|
|
||||||
bool isInUK(Point3D p);
|
|
||||||
Loading…
Reference in a new issue