David Luff:
The patches deal with three separate issues, all rolled up into one tarball: Currently, arrayfit always appends .arr.gz onto the name passed on the command line, meaning that only tile names can be passed. The patch strips off .arr or .arr.gz if present prior to it's appending, meaning that tile names or filenames can be passed on the command line. The interface to the OSGB36 conversion functions is cleaned up a lot. I can't believe I originally wrote it in such an ugly manner! A lot of console output (> 5000 lines per tile) is removed from the final construction process, meaning that the output left can actually be read.
This commit is contained in:
parent
a71d62c859
commit
791b1bb4fc
10 changed files with 152 additions and 202 deletions
|
@ -209,9 +209,6 @@ bool TGClipper::load_osgb36_polys(const string& path) {
|
||||||
|
|
||||||
Point3D p;
|
Point3D p;
|
||||||
Point3D OSRef;
|
Point3D OSRef;
|
||||||
Point3D OSLatLon;
|
|
||||||
Point3D OSCartesian;
|
|
||||||
Point3D WGS84Cartesian;
|
|
||||||
in >> skipcomment;
|
in >> skipcomment;
|
||||||
while ( !in.eof() ) {
|
while ( !in.eof() ) {
|
||||||
in >> poly_name;
|
in >> poly_name;
|
||||||
|
@ -237,14 +234,11 @@ bool TGClipper::load_osgb36_polys(const string& path) {
|
||||||
in >> startx;
|
in >> startx;
|
||||||
in >> starty;
|
in >> starty;
|
||||||
OSRef = Point3D(startx, starty, -9999.0);
|
OSRef = Point3D(startx, starty, -9999.0);
|
||||||
|
|
||||||
//Convert from OSGB36 Eastings/Northings to WGS84 Lat/Lon
|
//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
|
//Note that startx and starty themselves must not be altered since we compare them with unaltered lastx and lasty later
|
||||||
OSLatLon = ConvertEastingsNorthingsToLatLon(OSRef);
|
p = OSGB36ToWGS84(OSRef);
|
||||||
OSCartesian = ConvertAiry1830PolarToCartesian(OSLatLon);
|
|
||||||
WGS84Cartesian = ConvertOSGB36ToWGS84(OSCartesian);
|
|
||||||
p = ConvertGRS80CartesianToPolar(WGS84Cartesian);
|
|
||||||
|
|
||||||
poly.add_node( i, p );
|
poly.add_node( i, p );
|
||||||
SG_LOG( SG_CLIPPER, SG_BULK, "0 = "
|
SG_LOG( SG_CLIPPER, SG_BULK, "0 = "
|
||||||
<< startx << ", " << starty );
|
<< startx << ", " << starty );
|
||||||
|
@ -253,12 +247,8 @@ bool TGClipper::load_osgb36_polys(const string& path) {
|
||||||
in >> x;
|
in >> x;
|
||||||
in >> y;
|
in >> y;
|
||||||
OSRef = Point3D( x, y, -9999.0 );
|
OSRef = Point3D( x, y, -9999.0 );
|
||||||
|
p = OSGB36ToWGS84(OSRef);
|
||||||
OSLatLon = ConvertEastingsNorthingsToLatLon(OSRef);
|
|
||||||
OSCartesian = ConvertAiry1830PolarToCartesian(OSLatLon);
|
|
||||||
WGS84Cartesian = ConvertOSGB36ToWGS84(OSCartesian);
|
|
||||||
p = ConvertGRS80CartesianToPolar(WGS84Cartesian);
|
|
||||||
|
|
||||||
poly.add_node( i, p );
|
poly.add_node( i, p );
|
||||||
SG_LOG( SG_CLIPPER, SG_BULK, j << " = " << x << ", " << y );
|
SG_LOG( SG_CLIPPER, SG_BULK, j << " = " << x << ", " << y );
|
||||||
}
|
}
|
||||||
|
@ -271,12 +261,8 @@ bool TGClipper::load_osgb36_polys(const string& path) {
|
||||||
// last point same as first, discard
|
// last point same as first, discard
|
||||||
} else {
|
} else {
|
||||||
OSRef = Point3D( lastx, lasty, -9999.0 );
|
OSRef = Point3D( lastx, lasty, -9999.0 );
|
||||||
|
p = OSGB36ToWGS84(OSRef);
|
||||||
OSLatLon = ConvertEastingsNorthingsToLatLon(OSRef);
|
|
||||||
OSCartesian = ConvertAiry1830PolarToCartesian(OSLatLon);
|
|
||||||
WGS84Cartesian = ConvertOSGB36ToWGS84(OSCartesian);
|
|
||||||
p = ConvertGRS80CartesianToPolar(WGS84Cartesian);
|
|
||||||
|
|
||||||
poly.add_node( i, p );
|
poly.add_node( i, p );
|
||||||
SG_LOG( SG_CLIPPER, SG_BULK, count - 1 << " = "
|
SG_LOG( SG_CLIPPER, SG_BULK, count - 1 << " = "
|
||||||
<< lastx << ", " << lasty );
|
<< lastx << ", " << lasty );
|
||||||
|
|
|
@ -574,17 +574,17 @@ static void fix_point_heights( TGConstruct& c, const TGArray& array )
|
||||||
point_list raw_nodes = c.get_tri_nodes().get_node_list();
|
point_list raw_nodes = c.get_tri_nodes().get_node_list();
|
||||||
|
|
||||||
for ( i = 0; i < (int)raw_nodes.size(); ++i ) {
|
for ( i = 0; i < (int)raw_nodes.size(); ++i ) {
|
||||||
cout << " fixing = " << raw_nodes[i] << " ";
|
//cout << " fixing = " << raw_nodes[i] << " ";
|
||||||
int index = c.get_fixed_elevations().find( raw_nodes[i] );
|
int index = c.get_fixed_elevations().find( raw_nodes[i] );
|
||||||
if ( index >= 0 ) {
|
if ( index >= 0 ) {
|
||||||
// found an elevation we want to preserve
|
// found an elevation we want to preserve
|
||||||
z = c.get_fixed_elevations().get_node(index).z();
|
z = c.get_fixed_elevations().get_node(index).z();
|
||||||
cout << " forced = " << z << endl;
|
//cout << " forced = " << z << endl;
|
||||||
} else {
|
} else {
|
||||||
// interpolate point from DEM data.
|
// interpolate point from DEM data.
|
||||||
z = array.altitude_from_grid( raw_nodes[i].x() * 3600.0,
|
z = array.altitude_from_grid( raw_nodes[i].x() * 3600.0,
|
||||||
raw_nodes[i].y() * 3600.0 );
|
raw_nodes[i].y() * 3600.0 );
|
||||||
cout << " interpolated = " << z << endl;
|
//cout << " interpolated = " << z << endl;
|
||||||
}
|
}
|
||||||
raw_nodes[i].setz( z );
|
raw_nodes[i].setz( z );
|
||||||
}
|
}
|
||||||
|
@ -747,14 +747,14 @@ static void fix_land_cover_assignments( TGConstruct& c ) {
|
||||||
// a different coverage for each vertex, just pick
|
// a different coverage for each vertex, just pick
|
||||||
// from the middle/average
|
// from the middle/average
|
||||||
Point3D average = ( p1 + p2 + p3 ) / 3.0;
|
Point3D average = ( p1 + p2 + p3 ) / 3.0;
|
||||||
cout << " average triangle center = " << average;
|
//cout << " average triangle center = " << average;
|
||||||
new_area = get_area_type( c.get_cover(),
|
new_area = get_area_type( c.get_cover(),
|
||||||
average.x() - quarter_cover_size,
|
average.x() - quarter_cover_size,
|
||||||
average.y() - quarter_cover_size,
|
average.y() - quarter_cover_size,
|
||||||
cover_size, cover_size );
|
cover_size, cover_size );
|
||||||
}
|
}
|
||||||
|
|
||||||
cout << " new attrib = " << get_area_name( new_area ) << endl;
|
//cout << " new attrib = " << get_area_name( new_area ) << endl;
|
||||||
c.set_tri_attribute( i, new_area );
|
c.set_tri_attribute( i, new_area );
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -912,7 +912,7 @@ static point_list gen_face_normals( TGConstruct& c ) {
|
||||||
triele_list tri_elements = c.get_tri_elements();
|
triele_list tri_elements = c.get_tri_elements();
|
||||||
for ( int i = 0; i < (int)tri_elements.size(); i++ ) {
|
for ( int i = 0; i < (int)tri_elements.size(); i++ ) {
|
||||||
Point3D p = calc_normal(c, i );
|
Point3D p = calc_normal(c, i );
|
||||||
cout << p << endl;
|
//cout << p << endl;
|
||||||
face_normals.push_back( p );
|
face_normals.push_back( p );
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -950,7 +950,7 @@ static point_list gen_point_normals( TGConstruct& c ) {
|
||||||
}
|
}
|
||||||
|
|
||||||
average /= total_area;
|
average /= total_area;
|
||||||
cout << "average = " << average << endl;
|
//cout << "average = " << average << endl;
|
||||||
|
|
||||||
point_normals.push_back( average );
|
point_normals.push_back( average );
|
||||||
}
|
}
|
||||||
|
|
|
@ -423,6 +423,7 @@ void TGMatch::split_tile( TGConstruct& c ) {
|
||||||
if ( !se_flag ) { cout << " se corner = " << se_node << endl; }
|
if ( !se_flag ) { cout << " se corner = " << se_node << endl; }
|
||||||
if ( !ne_flag ) { cout << " ne corner = " << ne_node << endl; }
|
if ( !ne_flag ) { cout << " ne corner = " << ne_node << endl; }
|
||||||
if ( !nw_flag ) { cout << " nw corner = " << nw_node << endl; }
|
if ( !nw_flag ) { cout << " nw corner = " << nw_node << endl; }
|
||||||
|
/*
|
||||||
if ( !north_flag ) {
|
if ( !north_flag ) {
|
||||||
cout << " north nodes = " << north_nodes.size() << endl;
|
cout << " north nodes = " << north_nodes.size() << endl;
|
||||||
for ( i = 0; i < (int)north_nodes.size(); ++i ) {
|
for ( i = 0; i < (int)north_nodes.size(); ++i ) {
|
||||||
|
@ -451,6 +452,7 @@ void TGMatch::split_tile( TGConstruct& c ) {
|
||||||
for ( i = 0; i < (int)body_nodes.size(); ++i ) {
|
for ( i = 0; i < (int)body_nodes.size(); ++i ) {
|
||||||
cout << " " << body_nodes[i] << endl;
|
cout << " " << body_nodes[i] << endl;
|
||||||
}
|
}
|
||||||
|
*/
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
|
@ -4,7 +4,7 @@ libOsgb36_a_SOURCES = osgb36.cxx osgb36.hxx osgbtc.cxx osgbtc.hxx uk.cxx uk.hxx
|
||||||
|
|
||||||
noinst_PROGRAMS = testosgb36
|
noinst_PROGRAMS = testosgb36
|
||||||
|
|
||||||
testosgb36_SOURCES = testosgb36.cxx
|
testosgb36_SOURCES = testosgb36.cxx osgb36.cxx
|
||||||
|
|
||||||
INCLUDES = \
|
INCLUDES = \
|
||||||
-I$(top_srcdir) \
|
-I$(top_srcdir) \
|
||||||
|
|
|
@ -19,11 +19,83 @@
|
||||||
// along with this program; if not, write to the Free Software
|
// along with this program; if not, write to the Free Software
|
||||||
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
||||||
|
|
||||||
|
#include <simgear/compiler.h>
|
||||||
|
|
||||||
|
#include STL_IOSTREAM
|
||||||
|
#include <math.h>
|
||||||
|
|
||||||
#include "osgb36.hxx"
|
#include "osgb36.hxx"
|
||||||
|
|
||||||
|
SG_USING_STD(cout);
|
||||||
|
|
||||||
double DEG_TO_RAD = 2.0 * 3.14159265358979323846264338327950288419716939967511 / 360.0;
|
double DEG_TO_RAD = 2.0 * 3.14159265358979323846264338327950288419716939967511 / 360.0;
|
||||||
double RAD_TO_DEG = 360.0 / (2.0 * 3.14159265358979323846264338327950288419716939967511);
|
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)
|
//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference)
|
||||||
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon)
|
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon)
|
||||||
{
|
{
|
||||||
|
|
|
@ -19,50 +19,18 @@
|
||||||
// along with this program; if not, write to the Free Software
|
// along with this program; if not, write to the Free Software
|
||||||
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
|
||||||
|
|
||||||
#include <simgear/compiler.h>
|
|
||||||
|
|
||||||
#include STL_IOSTREAM
|
|
||||||
#include <math.h>
|
|
||||||
#include <simgear/math/point3d.hxx>
|
#include <simgear/math/point3d.hxx>
|
||||||
|
|
||||||
SG_USING_STD(cout);
|
|
||||||
SG_USING_STD(endl);
|
|
||||||
|
|
||||||
//*******************************************************************
|
//*******************************************************************
|
||||||
//
|
//
|
||||||
// NOTE All functions which take or return Latitude and Longitude
|
// NOTE All functions which take or return Latitude and Longitude
|
||||||
// take and return degrees. Conversion to and from radians
|
// take and return degrees. XYZ coordinates and ellipsoid
|
||||||
// is performed internaly. XYZ coordinates and ellipsoid
|
|
||||||
// heights are returned and required in meters.
|
// heights are returned and required in meters.
|
||||||
//
|
//
|
||||||
//*******************************************************************
|
//*******************************************************************
|
||||||
|
|
||||||
|
// Convert WGS84 Lat/Lon (degrees) to OSGB36 Eastings and Northings (meters)
|
||||||
|
Point3D WGS84ToOSGB36(Point3D p);
|
||||||
|
|
||||||
//Convert OSGB36 Lat and Lon to OSGB36 Eastings and Northings (Grid Reference)
|
// Convert OSGB36 Eastings and Northings (meters) to WGS84 Lat/Lon (degrees)
|
||||||
Point3D ConvertLatLonToEastingsNorthings(Point3D LatLon);
|
Point3D OSGB36ToWGS84(Point3D p);
|
||||||
|
|
||||||
//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);
|
|
||||||
|
|
||||||
|
|
|
@ -9,67 +9,20 @@
|
||||||
point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
|
point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
|
||||||
const int_list& fan, double scale )
|
const int_list& fan, double scale )
|
||||||
{
|
{
|
||||||
// cout << "calculating texture coordinates for a specific fan of size = "
|
|
||||||
// << fan.size() << endl;
|
|
||||||
|
|
||||||
// calculate perimeter based on center of this degree (not center
|
|
||||||
// of bucket)
|
|
||||||
/* double clat = (int)b.get_center_lat();
|
|
||||||
if ( clat > 0 ) {
|
|
||||||
clat = (int)clat + 0.5;
|
|
||||||
} else {
|
|
||||||
clat = (int)clat - 0.5;
|
|
||||||
}
|
|
||||||
*/
|
|
||||||
// double clat_rad = clat * DEG_TO_RAD;
|
|
||||||
// double cos_lat = cos( clat_rad );
|
|
||||||
// double local_radius = cos_lat * EQUATORIAL_RADIUS_M;
|
|
||||||
// double local_perimeter = 2.0 * local_radius * FG_PI;
|
|
||||||
// double degree_width = local_perimeter / 360.0;
|
|
||||||
|
|
||||||
// cout << "clat = " << clat << endl;
|
|
||||||
// cout << "clat (radians) = " << clat_rad << endl;
|
|
||||||
// cout << "cos(lat) = " << cos_lat << endl;
|
|
||||||
// cout << "local_radius = " << local_radius << endl;
|
|
||||||
// cout << "local_perimeter = " << local_perimeter << endl;
|
|
||||||
// cout << "degree_width = " << degree_width << endl;
|
|
||||||
|
|
||||||
// double perimeter = 2.0 * EQUATORIAL_RADIUS_M * FG_PI;
|
|
||||||
// double degree_height = perimeter / 360.0;
|
|
||||||
// cout << "degree_height = " << degree_height << endl;
|
|
||||||
|
|
||||||
// find min/max of fan
|
// find min/max of fan
|
||||||
Point3D tmin, tmax, WGS84LatLon, t;
|
Point3D t, tmin, tmax;
|
||||||
Point3D WGS84xyz, OSGB36xyz, OSGB36LatLon;
|
|
||||||
Point3D OSGridRef;
|
Point3D OSGridRef;
|
||||||
bool first = true;
|
bool first = true;
|
||||||
|
|
||||||
int i;
|
int i;
|
||||||
|
|
||||||
for ( i = 0; i < (int)fan.size(); ++i ) {
|
for ( i = 0; i < (int)fan.size(); ++i ) {
|
||||||
WGS84LatLon = geod_nodes[ fan[i] ];
|
|
||||||
// cout << "point WGS84LatLon = " << WGS84LatLon << endl;
|
OSGridRef = WGS84ToOSGB36(geod_nodes[ fan[i] ]);
|
||||||
WGS84xyz = ConvertGRS80PolarToCartesian(WGS84LatLon);
|
|
||||||
// cout << "WGS84XYZ = " << WGS84XYZ.x << ", " << WGS84XYZ.y << ", " << WGS84XYZ.z << '\n';
|
|
||||||
OSGB36xyz = ConvertWGS84ToOSGB36(WGS84xyz);
|
|
||||||
// cout << "OSXYZ = " << OSXYZ.x << ", " << OSXYZ.y << ", " << OSXYZ.z << '\n';
|
|
||||||
OSGB36LatLon = ConvertAiry1830CartesianToPolar(OSGB36xyz);
|
|
||||||
// cout << "OSGB36LatLon = " << OSGB36LatLon.x() << ", " << OSGB36LatLon.y() << ", " << OSGB36LatLon.z() << '\n';
|
|
||||||
OSGridRef = ConvertLatLonToEastingsNorthings(OSGB36LatLon);
|
|
||||||
// cout << "OS Eastings and Northings = " << OSGridRef << '\n';
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
t = UK_basic_tex_coord( OSGridRef );
|
t = UK_basic_tex_coord( OSGridRef );
|
||||||
// cout << "basic_tex_coord = " << t << endl;
|
// cout << "basic_tex_coord = " << t << endl;
|
||||||
|
|
||||||
// cout << "p = " << p << '\n';
|
|
||||||
// cout << "t = " << t << '\n';
|
|
||||||
// cout << "Degree width = " << degree_width << '\n';
|
|
||||||
// cout << "Degree height = " << degree_height << '\n';
|
|
||||||
// char dcl_pause;
|
|
||||||
// cin >> dcl_pause;
|
|
||||||
|
|
||||||
if ( first ) {
|
if ( first ) {
|
||||||
tmin = tmax = t;
|
tmin = tmax = t;
|
||||||
first = false;
|
first = false;
|
||||||
|
@ -171,17 +124,8 @@ point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
|
||||||
point_list tex;
|
point_list tex;
|
||||||
tex.clear();
|
tex.clear();
|
||||||
for ( i = 0; i < (int)fan.size(); ++i ) {
|
for ( i = 0; i < (int)fan.size(); ++i ) {
|
||||||
WGS84LatLon = geod_nodes[ fan[i] ];
|
|
||||||
// cout << "About to do second run through tex coords\n";
|
OSGridRef = WGS84ToOSGB36(geod_nodes[ fan[i] ]);
|
||||||
// cout << "point WGS84LatLon = " << WGS84LatLon << endl;
|
|
||||||
WGS84xyz = ConvertGRS80PolarToCartesian(WGS84LatLon);
|
|
||||||
// cout << "WGS84XYZ = " << WGS84XYZ.x << ", " << WGS84XYZ.y << ", " << WGS84XYZ.z << '\n';
|
|
||||||
OSGB36xyz = ConvertWGS84ToOSGB36(WGS84xyz);
|
|
||||||
// cout << "OSXYZ = " << OSXYZ.x << ", " << OSXYZ.y << ", " << OSXYZ.z << '\n';
|
|
||||||
OSGB36LatLon = ConvertAiry1830CartesianToPolar(OSGB36xyz);
|
|
||||||
// cout << "OSGB36LatLon = " << OSGB36LatLon.x() << ", " << OSGB36LatLon.y() << ", " << OSGB36LatLon.z() << '\n';
|
|
||||||
OSGridRef = ConvertLatLonToEastingsNorthings(OSGB36LatLon);
|
|
||||||
// cout << "OS Eastings and Northings = " << OSGridRef << '\n';
|
|
||||||
|
|
||||||
t = UK_basic_tex_coord(OSGridRef);
|
t = UK_basic_tex_coord(OSGridRef);
|
||||||
// cout << "second t = " << t << endl;
|
// cout << "second t = " << t << endl;
|
||||||
|
@ -209,7 +153,7 @@ point_list UK_calc_tex_coords( const SGBucket& b, const point_list& geod_nodes,
|
||||||
adjusted_t.sety( 0.0 );
|
adjusted_t.sety( 0.0 );
|
||||||
}
|
}
|
||||||
adjusted_t.setz( 0.0 );
|
adjusted_t.setz( 0.0 );
|
||||||
cout << "adjusted_t after check for SG_EPSILON = " << adjusted_t << endl;
|
//cout << "adjusted_t after check for SG_EPSILON = " << adjusted_t << endl;
|
||||||
|
|
||||||
tex.push_back( adjusted_t );
|
tex.push_back( adjusted_t );
|
||||||
}
|
}
|
||||||
|
|
|
@ -1,74 +1,38 @@
|
||||||
|
#include <simgear/compiler.h>
|
||||||
#include <simgear/math/point3d.hxx>
|
#include <simgear/math/point3d.hxx>
|
||||||
|
|
||||||
|
#include STL_IOSTREAM
|
||||||
|
|
||||||
#include "osgb36.hxx"
|
#include "osgb36.hxx"
|
||||||
|
|
||||||
//Test harness for various WGS84 and OSGB36 conversion functions
|
SG_USING_STD(cout);
|
||||||
int main()
|
|
||||||
|
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)
|
||||||
{
|
{
|
||||||
cout << "Test OSGB36\n";
|
if(argc != 3) Usage();
|
||||||
/* Point3D initial_position;
|
|
||||||
Point3D final_position;
|
cout << "Test OSGB36\n";
|
||||||
|
|
||||||
initial_position.setlat(52.0 + ((39.0 + 27.2531/60.0)/60.0));
|
Point3D p1;
|
||||||
initial_position.setlon(1.0 + ((43.0 + 4.5177/60.0)/60.0));
|
p1.setlon(atof(argv[1]));
|
||||||
initial_position.setelev(0.0);
|
p1.setlat(atof(argv[2]));
|
||||||
|
p1.setelev(0.0);
|
||||||
final_position = ConvertLatLonToEastingsNorthings(initial_position);
|
|
||||||
|
cout << "WGS84 position is " << p1 << '\n';
|
||||||
cout << "Eastings = " << final_position.x() << '\n';
|
|
||||||
cout << "Northings = " << final_position.y() << '\n';
|
Point3D p2 = WGS84ToOSGB36(p1);
|
||||||
|
|
||||||
Point3D check_position;
|
cout << "OS coords are " << p2 << '\n';
|
||||||
|
|
||||||
check_position = ConvertEastingsNorthingsToLatLon(final_position);
|
p1 = OSGB36ToWGS84(p2);
|
||||||
|
|
||||||
cout << "Lat = " << check_position.lat() << '\n';
|
cout << "Conversion back to WGS84 gives " << p1 << '\n';
|
||||||
cout << "Lon = " << check_position.lon() << '\n';
|
|
||||||
cout << "Elipsiodal height = " << check_position.elev() << '\n';
|
|
||||||
|
|
||||||
check_position.setelev(24.7);
|
|
||||||
|
|
||||||
Point3D cartesian_position = ConvertAiry1830PolarToCartesian(check_position);
|
|
||||||
|
|
||||||
cout << "X = " << cartesian_position.x() << '\n';
|
|
||||||
cout << "Y = " << cartesian_position.y() << '\n';
|
|
||||||
cout << "Z = " << cartesian_position.z() << '\n';
|
|
||||||
|
|
||||||
cout << "\nTesting ConvertAiry1830CartesianToPolar.....\n";
|
|
||||||
|
|
||||||
Point3D XYZPosition;
|
|
||||||
|
|
||||||
XYZPosition.setx(3874938.849);
|
|
||||||
XYZPosition.sety(116218.624);
|
|
||||||
XYZPosition.setz(5047168.208);
|
|
||||||
|
|
||||||
Point3D LatLonPosition = ConvertAiry1830CartesianToPolar(XYZPosition);
|
|
||||||
|
|
||||||
cout << "Lat = " << LatLonPosition.lat() << '\n';
|
|
||||||
cout << "Lon = " << LatLonPosition.lon() << '\n';
|
|
||||||
cout << "Elipsiodal height = " << LatLonPosition.elev() << '\n';
|
|
||||||
|
|
||||||
cout << "\nNow attempting to convert WGS84 lat & lon to OS grid reference .....\n";
|
|
||||||
|
|
||||||
Point3D WGS84LatLon;
|
|
||||||
Point3D WGS84XYZ;
|
|
||||||
Point3D OSXYZ;
|
|
||||||
Point3D OSLatLon;
|
|
||||||
Point3D OSGridRef;
|
|
||||||
|
|
||||||
WGS84LatLon.setlat(52.0);
|
|
||||||
WGS84LatLon.setlon(-2.0);
|
|
||||||
WGS84LatLon.setelev(0.0);
|
|
||||||
|
|
||||||
WGS84XYZ = ConvertGRS80PolarToCartesian(WGS84LatLon);
|
|
||||||
cout << "WGS84XYZ = " << WGS84XYZ.x() << ", " << WGS84XYZ.y() << ", " << WGS84XYZ.z() << '\n';
|
|
||||||
OSXYZ = ConvertWGS84ToOSGB36(WGS84XYZ);
|
|
||||||
cout << "OSXYZ = " << OSXYZ.x() << ", " << OSXYZ.y() << ", " << OSXYZ.z() << '\n';
|
|
||||||
OSLatLon = ConvertAiry1830CartesianToPolar(OSXYZ);
|
|
||||||
cout << "OSLatLon = " << OSLatLon.lat() << ", " << OSLatLon.lon() << ", " << OSLatLon.elev() << '\n';
|
|
||||||
OSGridRef = ConvertLatLonToEastingsNorthings(OSLatLon);
|
|
||||||
|
|
||||||
cout << "OS Grid Reference = " << OSGridRef.x() << ", " << OSGridRef.y() << "\n\n";
|
|
||||||
*/
|
|
||||||
return(0);
|
return(0);
|
||||||
}
|
}
|
||||||
|
|
|
@ -133,9 +133,9 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes,
|
||||||
y_err = fabs(current->y() - (m * current->x() + b));
|
y_err = fabs(current->y() - (m * current->x() + b));
|
||||||
|
|
||||||
if ( y_err < FG_PROXIMITY_EPSILON ) {
|
if ( y_err < FG_PROXIMITY_EPSILON ) {
|
||||||
cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl;
|
//cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl;
|
||||||
cout << p0 << " < " << *current << " < "
|
//cout << p0 << " < " << *current << " < "
|
||||||
<< p1 << endl;
|
// << p1 << endl;
|
||||||
found_extra = true;
|
found_extra = true;
|
||||||
if ( y_err < y_err_min ) {
|
if ( y_err < y_err_min ) {
|
||||||
extra_index = counter;
|
extra_index = counter;
|
||||||
|
@ -182,9 +182,9 @@ void TGTriSegments::unique_divide_and_add( const point_list& nodes,
|
||||||
// }
|
// }
|
||||||
|
|
||||||
if ( x_err < FG_PROXIMITY_EPSILON ) {
|
if ( x_err < FG_PROXIMITY_EPSILON ) {
|
||||||
cout << "FOUND EXTRA SEGMENT NODE (X)" << endl;
|
//cout << "FOUND EXTRA SEGMENT NODE (X)" << endl;
|
||||||
cout << p0 << " < " << *current << " < "
|
//cout << p0 << " < " << *current << " < "
|
||||||
<< p1 << endl;
|
// << p1 << endl;
|
||||||
found_extra = true;
|
found_extra = true;
|
||||||
if ( x_err < x_err_min ) {
|
if ( x_err < x_err_min ) {
|
||||||
extra_index = counter;
|
extra_index = counter;
|
||||||
|
|
|
@ -146,6 +146,20 @@ int main( int argc, char **argv ) {
|
||||||
if ( ! infile.str().length() ) {
|
if ( ! infile.str().length() ) {
|
||||||
usage( argv[0] );
|
usage( argv[0] );
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// Strip off .arr or .arr.gz if part of the input file name
|
||||||
|
string s = infile.str();
|
||||||
|
if(s.length() >= 7) {
|
||||||
|
if(s.substr(s.length() - 7) == ".arr.gz") {
|
||||||
|
s = s.substr(0, s.length() - 7);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if(s.length() >= 4) {
|
||||||
|
if(s.substr(s.length() - 4) == ".arr") {
|
||||||
|
s = s.substr(0, s.length() - 4);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
infile.set(s);
|
||||||
|
|
||||||
SGPath outfile = infile;
|
SGPath outfile = infile;
|
||||||
outfile.concat( ".fit.gz" );
|
outfile.concat( ".fit.gz" );
|
||||||
|
|
Loading…
Reference in a new issue