diff --git a/src/BuildTiles/CMakeLists.txt b/src/BuildTiles/CMakeLists.txt index 0fa541e3..5787061e 100644 --- a/src/BuildTiles/CMakeLists.txt +++ b/src/BuildTiles/CMakeLists.txt @@ -1,8 +1,5 @@ - - include_directories(${PROJECT_SOURCE_DIR}/src/Lib) include_directories(${PROJECT_SOURCE_DIR}/src/BuildTiles) -add_subdirectory(Osgb36) add_subdirectory(Parallel) add_subdirectory(Main) diff --git a/src/BuildTiles/Main/CMakeLists.txt b/src/BuildTiles/Main/CMakeLists.txt index e5f44bc8..4210bb99 100644 --- a/src/BuildTiles/Main/CMakeLists.txt +++ b/src/BuildTiles/Main/CMakeLists.txt @@ -27,7 +27,6 @@ set_target_properties(tg-construct PROPERTIES "DEFAULT_USGS_MAPFILE=\"${PKGDATADIR}/usgsmap.txt\";DEFAULT_PRIORITIES_FILE=\"${PKGDATADIR}/default_priorities.txt\"" ) target_link_libraries(tg-construct - Osgb36 Polygon Geometry Array landcover poly2tri ${GDAL_LIBRARY} diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index 9cc0ec7e..65be2329 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -143,7 +143,6 @@ static void usage( const string name ) { SG_LOG(SG_GENERAL, SG_ALERT, " --nudge="); SG_LOG(SG_GENERAL, SG_ALERT, " --priorities="); SG_LOG(SG_GENERAL, SG_ALERT, " --usgs-map="); - SG_LOG(SG_GENERAL, SG_ALERT, " --useUKgrid"); SG_LOG(SG_GENERAL, SG_ALERT, " --ignore-landmass"); SG_LOG(SG_GENERAL, SG_ALERT, " ] "); exit(-1); @@ -166,10 +165,6 @@ int main(int argc, char **argv) { vector debug_shape_defs; vector debug_area_defs; - // flag indicating whether UK grid should be used for in-UK - // texture coordinate generation - bool useUKgrid = false; - bool ignoreLandmass = false; sglog().setLogLevels( SG_ALL, SG_INFO ); @@ -208,8 +203,6 @@ int main(int argc, char **argv) { priorities_file = arg.substr(13); } else if (arg.find("--usgs-map=") == 0) { usgs_map_file = arg.substr(11); - } else if (arg.find("--useUKgrid") == 0) { - useUKgrid = true; } else if (arg.find("--ignore-landmass") == 0) { ignoreLandmass = true; } else if (arg.find("--debug-dir=") == 0) { @@ -267,7 +260,7 @@ int main(int argc, char **argv) { all_stages = new TGConstruct(); all_stages->set_cover( cover ); 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_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->set_cover( cover ); 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_debug( debug_dir, debug_area_defs, debug_shape_defs ); @@ -325,7 +318,7 @@ int main(int argc, char **argv) { stage1 = new TGConstruct(); stage1->set_cover( cover ); 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_debug( debug_dir, debug_area_defs, debug_shape_defs ); @@ -354,7 +347,7 @@ int main(int argc, char **argv) { stage2 = new TGConstruct(); stage2->set_cover( cover ); 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_debug( debug_dir, debug_area_defs, debug_shape_defs ); @@ -384,7 +377,7 @@ int main(int argc, char **argv) { stage3 = new TGConstruct(); stage3->set_cover( cover ); 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_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->set_cover( cover ); 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_debug( debug_dir, debug_area_defs, debug_shape_defs ); diff --git a/src/BuildTiles/Main/tgconstruct.cxx b/src/BuildTiles/Main/tgconstruct.cxx index 557d14b9..cd01b92c 100644 --- a/src/BuildTiles/Main/tgconstruct.cxx +++ b/src/BuildTiles/Main/tgconstruct.cxx @@ -50,7 +50,6 @@ // Constructor TGConstruct::TGConstruct(): - useUKGrid(false), ignoreLandmass(false), debug_all(false), ds_id((void*)-1) @@ -77,8 +76,7 @@ void TGConstruct::set_paths( const std::string work, const std::string share, co load_dirs = load; } -void TGConstruct::set_options( bool uk_grid, bool ignore_lm, double n ) { - useUKGrid = uk_grid; +void TGConstruct::set_options( bool ignore_lm, double n ) { ignoreLandmass = ignore_lm; nudge = n; } diff --git a/src/BuildTiles/Main/tgconstruct.hxx b/src/BuildTiles/Main/tgconstruct.hxx index f29c4847..ff90dee5 100644 --- a/src/BuildTiles/Main/tgconstruct.hxx +++ b/src/BuildTiles/Main/tgconstruct.hxx @@ -87,10 +87,6 @@ private: 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 bool ignoreLandmass; @@ -136,7 +132,6 @@ private: int LoadLandclassPolys( void ); // Load Data Helpers 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); // Clip Data @@ -239,7 +234,7 @@ public: inline void set_load_dirs( const std::vector ld ) { load_dirs = ld; } #endif - void set_options( bool uk_grid, bool ignore_lm, double n ); + void set_options( bool ignore_lm, double n ); #if 0 // UK grid flag diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx index 22e60e46..fb990284 100644 --- a/src/BuildTiles/Main/tgconstruct_poly.cxx +++ b/src/BuildTiles/Main/tgconstruct_poly.cxx @@ -30,7 +30,6 @@ #include #include -#include #include "tgconstruct.hxx" @@ -272,102 +271,6 @@ bool TGConstruct::load_poly(const string& path) { 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 // clip against each other to resolve any overlaps int TGConstruct::LoadLandclassPolys( void ) { @@ -403,10 +306,6 @@ int TGConstruct::LoadLandclassPolys( void ) { (lext == "fit") || (lext == "fit.gz") || (lext == "ind")) { // 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 { load_poly( p.str() ); SG_LOG(SG_GENERAL, SG_ALERT, " Loaded " << p.file()); diff --git a/src/BuildTiles/Osgb36/.gitignore b/src/BuildTiles/Osgb36/.gitignore deleted file mode 100644 index 5e3afa89..00000000 --- a/src/BuildTiles/Osgb36/.gitignore +++ /dev/null @@ -1 +0,0 @@ -testosgb36 diff --git a/src/BuildTiles/Osgb36/CMakeLists.txt b/src/BuildTiles/Osgb36/CMakeLists.txt deleted file mode 100644 index a8a91bf8..00000000 --- a/src/BuildTiles/Osgb36/CMakeLists.txt +++ /dev/null @@ -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) diff --git a/src/BuildTiles/Osgb36/osgb36.cxx b/src/BuildTiles/Osgb36/osgb36.cxx deleted file mode 100644 index 9ddbb6c4..00000000 --- a/src/BuildTiles/Osgb36/osgb36.cxx +++ /dev/null @@ -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 - -#include -#include - -#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); -} - - - - - - - diff --git a/src/BuildTiles/Osgb36/osgb36.hxx b/src/BuildTiles/Osgb36/osgb36.hxx deleted file mode 100644 index fb176c55..00000000 --- a/src/BuildTiles/Osgb36/osgb36.hxx +++ /dev/null @@ -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 - -//******************************************************************* -// -// 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); diff --git a/src/BuildTiles/Osgb36/osgbtc.cxx b/src/BuildTiles/Osgb36/osgbtc.cxx deleted file mode 100644 index ab77b47e..00000000 --- a/src/BuildTiles/Osgb36/osgbtc.cxx +++ /dev/null @@ -1,168 +0,0 @@ -#ifdef HAVE_CONFIG_H -# include -#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; -} diff --git a/src/BuildTiles/Osgb36/osgbtc.hxx b/src/BuildTiles/Osgb36/osgbtc.hxx deleted file mode 100644 index 7438b7e6..00000000 --- a/src/BuildTiles/Osgb36/osgbtc.hxx +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef OSBTEXCOORD_H -#define OSBTEXCOORD_H - - -#include -#include -#include - - -#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 diff --git a/src/BuildTiles/Osgb36/testosgb36.cxx b/src/BuildTiles/Osgb36/testosgb36.cxx deleted file mode 100644 index baca6d5c..00000000 --- a/src/BuildTiles/Osgb36/testosgb36.cxx +++ /dev/null @@ -1,39 +0,0 @@ -#include -#include - -#include -#include - -#include "osgb36.hxx" - -using std::cout; - -static void Usage() { - cout << "Usage is testosgb36 \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); -} diff --git a/src/BuildTiles/Osgb36/uk.cxx b/src/BuildTiles/Osgb36/uk.cxx deleted file mode 100644 index 754f36ce..00000000 --- a/src/BuildTiles/Osgb36/uk.cxx +++ /dev/null @@ -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); -} diff --git a/src/BuildTiles/Osgb36/uk.hxx b/src/BuildTiles/Osgb36/uk.hxx deleted file mode 100644 index 6ee76585..00000000 --- a/src/BuildTiles/Osgb36/uk.hxx +++ /dev/null @@ -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 - -//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);