diff --git a/src/Airports/GenAirports/build.cxx b/src/Airports/GenAirports/build.cxx index cbbae7e8..84f66fc4 100644 --- a/src/Airports/GenAirports/build.cxx +++ b/src/Airports/GenAirports/build.cxx @@ -61,10 +61,6 @@ #include "texparams.hxx" -static const double tgAirportEpsilon = FG_EPSILON / 10.0; -// static const double tgAirportEpsilon = FG_EPSILON * 10.0; - - // calculate texture coordinates for runway section using the provided // texturing parameters. Returns a mirror polygon to the runway, // except each point is the texture coordinate of the corresponding @@ -155,197 +151,6 @@ static FGPolygon rwy_section_tex_coords( const FGPolygon& in_poly, } -// Find a point in the given node list that lies between start and -// end, return true if something found, false if nothing found. -bool find_intermediate_node( const Point3D& start, const Point3D& end, - const point_list& nodes, Point3D *result ) -{ - bool found_node = false; - double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min; - - Point3D p0 = start; - Point3D p1 = end; - - // cout << " add_intermediate_nodes()" << endl; - printf(" %.7f %.7f %.7f <=> %.7f %.7f %.7f\n", - p0.x(), p0.y(), p0.z(), p1.x(), p1.y(), p1.z() ); - - double xdist = fabs(p0.x() - p1.x()); - double ydist = fabs(p0.y() - p1.y()); - // cout << "xdist = " << xdist << " ydist = " << ydist << endl; - x_err_min = xdist + 1.0; - y_err_min = ydist + 1.0; - - if ( xdist > ydist ) { - // cout << "use y = mx + b" << endl; - - // sort these in a sensible order - Point3D p_min, p_max; - if ( p0.x() < p1.x() ) { - p_min = p0; - p_max = p1; - } else { - p_min = p1; - p_max = p0; - } - - m = (p_min.y() - p_max.y()) / (p_min.x() - p_max.x()); - b = p_max.y() - m * p_max.x(); - - // cout << "m = " << m << " b = " << b << endl; - - for ( int i = 0; i < (int)nodes.size(); ++i ) { - // cout << i << endl; - Point3D current = nodes[i]; - - if ( (current.x() > (p_min.x() + FG_EPSILON)) - && (current.x() < (p_max.x() - FG_EPSILON)) ) { - - // printf( "found a potential candidate %.7f %.7f %.7f\n", - // current.x(), current.y(), current.z() ); - - y_err = fabs(current.y() - (m * current.x() + b)); - // cout << "y_err = " << y_err << endl; - - if ( y_err < tgAirportEpsilon ) { - // cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl; - // cout << p_min << " < " << current << " < " - // << p_max << endl; - found_node = true; - if ( y_err < y_err_min ) { - *result = current; - y_err_min = y_err; - } - } - } - } - } else { - // cout << "use x = m1 * y + b1" << endl; - - // sort these in a sensible order - Point3D p_min, p_max; - if ( p0.y() < p1.y() ) { - p_min = p0; - p_max = p1; - } else { - p_min = p1; - p_max = p0; - } - - m1 = (p_min.x() - p_max.x()) / (p_min.y() - p_max.y()); - b1 = p_max.x() - m1 * p_max.y(); - - // cout << " m1 = " << m1 << " b1 = " << b1 << endl; - // printf( " m = %.8f b = %.8f\n", 1/m1, -b1/m1); - - // cout << " should = 0 = " - // << fabs(p_min.x() - (m1 * p_min.y() + b1)) << endl; - // cout << " should = 0 = " - // << fabs(p_max.x() - (m1 * p_max.y() + b1)) << endl; - - for ( int i = 0; i < (int)nodes.size(); ++i ) { - Point3D current = nodes[i]; - - if ( (current.y() > (p_min.y() + FG_EPSILON)) - && (current.y() < (p_max.y() - FG_EPSILON)) ) { - - // printf( "found a potential candidate %.7f %.7f %.7f\n", - // current.x(), current.y(), current.z() ); - - x_err = fabs(current.x() - (m1 * current.y() + b1)); - // cout << "x_err = " << x_err << endl; - - // if ( temp ) { - // cout << " (" << counter << ") x_err = " << x_err << endl; - // } - - if ( x_err < tgAirportEpsilon ) { - // cout << "FOUND EXTRA SEGMENT NODE (X)" << endl; - // cout << p_min << " < " << current << " < " - // << p_max << endl; - found_node = true; - if ( x_err < x_err_min ) { - *result = current; - x_err_min = x_err; - } - } - } - } - } - - return found_node; -} - - -// Attempt to reduce degeneracies where a subsequent point of a -// polygon lies *on* a previous line segment. These artifacts are -// occasionally introduced by the gpc polygon clipper. -static point_list reduce_contour_degeneracy( const point_list& contour ) { - point_list result = contour; - - Point3D p0, p1, bad_node; - bool done = false; - - while ( !done ) { - // traverse the contour until we find the first bad node or - // hit the end of the contour - cout << " ... not done ... " << endl; - bool bad = false; - - int i = 0; - while ( i < (int)result.size() - 1 && !bad ) { - p0 = result[i]; - p1 = result[i+1]; - - bad = find_intermediate_node( p0, p1, result, &bad_node ); - ++i; - } - if ( !bad ) { - // do the end/start connecting segment - p0 = result[result.size() - 1]; - p1 = result[0]; - - bad = find_intermediate_node( p0, p1, result, &bad_node ); - } - - if ( bad ) { - // remove bad node from contour - point_list tmp; tmp.clear(); - for ( int j = 0; j < (int)result.size(); ++j ) { - if ( result[j] == bad_node ) { - // skip - } else { - tmp.push_back( result[j] ); - } - } - result = tmp; - } else { - done = true; - } - } - - return result; -} - -// Search each segment of each contour for degenerate points (i.e. out -// of order points that lie coincident on other segments - -static FGPolygon reduce_degeneracy( const FGPolygon& poly ) { - FGPolygon result; - - for ( int i = 0; i < poly.contours(); ++i ) { - point_list contour = poly.get_contour(i); - contour = reduce_contour_degeneracy( contour ); - result.add_contour( contour, poly.get_hole_flag(i) ); - - // maintain original hole flag setting - // result.set_hole_flag( i, poly.get_hole_flag( i ) ); - } - - return result; -} - - // Divide segment if there are other existing points on it, return the // new polygon void add_intermediate_nodes( int contour, const Point3D& start, @@ -426,50 +231,6 @@ static FGPolygon add_nodes_to_poly( const FGPolygon& poly, } -// remove any duplicate nodes -static FGPolygon remove_dups( const FGPolygon &poly ) { - FGPolygon result; - result.erase(); - - for ( int i = 0; i < poly.contours(); ++i ) { - Point3D last = poly.get_pt( i, poly.contour_size(i) - 1 ); - bool all_same = true; - for ( int j = 0; j < poly.contour_size(i); ++j ) { - // cout << " " << i << " " << j << endl; - Point3D cur = poly.get_pt( i, j ); - if ( cur == last ) { - // skip - // cout << "skipping a duplicate point" << endl; - } else { - result.add_node( i, cur ); - all_same = false; - last = cur; - } - } - - // make sure the last point doesn't equal the previous or the first. - Point3D begin = poly.get_pt( i, 0 ); - Point3D end = poly.get_pt( i, poly.contour_size(i) - 1 ); - if ( begin == end ) { - // skip - cout << "begin == end!" << endl; - // exit(-1); - } - - if ( !all_same ) { - int flag = poly.get_hole_flag( i ); - result.set_hole_flag( i, flag ); - } else { - // too small an area ... add a token point to the contour - // to keep other things happy - result.add_node( i, begin ); - } - } - - return result; -} - - // Traverse a polygon and split edges until they are less than max_len // (specified in meters) static FGPolygon split_long_edges( const FGPolygon &poly, double max_len ) { @@ -543,24 +304,6 @@ static FGPolygon split_long_edges( const FGPolygon &poly, double max_len ) { } -// remove any degenerate contours -static FGPolygon remove_bad_contours( const FGPolygon &poly ) { - FGPolygon result; - result.erase(); - - for ( int i = 0; i < poly.contours(); ++i ) { - point_list contour = poly.get_contour( i ); - if ( contour.size() >= 3 ) { - // good - int flag = poly.get_hole_flag( i ); - result.add_contour( contour, flag ); - } - } - - return result; -} - - // fix node elevations point_list calc_elevations( const string& root, const point_list& geod_nodes ) { bool done = false; diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index a810d089..d196578b 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -171,13 +171,43 @@ static int actual_load_polys( const string& dir, } -// generate polygons from land-cover raster. +// Merge a polygon with an existing one if possible, append a new +// one otherwise; this function is used by actual_load_landcover, below, +// to reduce the number of separate polygons. +static void inline add_to_polys (vector &list, + const FGPolygon &poly) +{ + int len = list.size(); + FGPolygon temp; + + // Look for an adjacent polygon + if (len > 0) { + for (int i = 0; i < len; i++) { + temp = polygon_union(list[i], poly); + if (temp.contours() == 1) { // no new contours: they are adjacent + cout << "Merging with existing land-cover polygon" << endl; + list[i] = temp; + return; + } + } + } + + // No adjacent polygons available; append + // this one to the list. + cout << "Adding new land-cover polygon" << endl; + list.push_back(poly); +} + + +// Generate polygons from land-cover raster. Horizontally- or vertically- +// adjacent polygons will be merged automatically. static int actual_load_landcover ( LandCover &cover, FGConstruct & c, FGClipper &clipper ) { int count = 0; double lon, lat; - FGPolygon poly; + vector poly_list[FG_MAX_AREA_TYPES]; + FGPolygon poly; // working polygon // Get the top corner of the tile lon = @@ -187,31 +217,53 @@ static int actual_load_landcover ( LandCover &cover, FGConstruct & c, cout << "DPM: tile at " << lon << ',' << lat << endl; - // FIXME: this may still be wrong + + // Figure out how many units wide and + // high this tile is; each unit is + // 30 arc seconds. int x_span = int(120 * bucket_span(lat)); // arcsecs of longitude int y_span = int(120 * FG_BUCKET_SPAN); // arcsecs of latitude for (int x = 0; x < x_span; x++) { for (int y = 0; y < y_span; y++) { + // Figure out the boundaries of the + // 30 arcsec square. double x1 = lon + (x * (1.0/120.0)); double y1 = lat + (y * (1.0/120.0)); double x2 = x1 + (1.0/120.0); double y2 = y1 + (1.0/120.0); + // Look up the land cover for the square int cover_value = cover.getValue(x1 + (1.0/240.0), y1 + (1.0/240.0)); cout << " position: " << x1 << ',' << y1 << ',' << cover.getDescUSGS(cover_value) << endl; AreaType area = translateUSGSCover(cover_value); if (area != DefaultArea) { + // Create a square polygon and merge + // it into the list. poly.erase(); poly.add_node(0, Point3D(x1, y1, 0.0)); poly.add_node(0, Point3D(x1, y2, 0.0)); poly.add_node(0, Point3D(x2, y2, 0.0)); poly.add_node(0, Point3D(x2, y1, 0.0)); - clipper.add_poly(area, poly); - count++; + add_to_polys(poly_list[area], poly); } } } + // Now that we're finished looking up + // land cover, we have a list of lists + // of polygons, one (possibly-empty) + // list for each area type. Add the + // remaining polygons to the clipper. + for (int i = 0; i < FG_MAX_AREA_TYPES; i++) { + int len = poly_list[i].size(); + for (int j = 0; j < len; j++) { + clipper.add_poly(i, poly_list[i][j]); + count++; + } + } + + // Return the number of polygons + // actually read. return count; } diff --git a/src/BuildTiles/Triangulate/triangle.cxx b/src/BuildTiles/Triangulate/triangle.cxx index 51d30d58..56a9d4ad 100644 --- a/src/BuildTiles/Triangulate/triangle.cxx +++ b/src/BuildTiles/Triangulate/triangle.cxx @@ -123,6 +123,13 @@ FGTriangle::build( const point_list& corner_list, // also inside any interior contours // new way + + // try to make sure our polygons aren't goofy + gpc_poly = remove_dups( gpc_poly ); + gpc_poly = reduce_degeneracy( gpc_poly ); + gpc_poly = remove_dups( gpc_poly ); + gpc_poly = remove_bad_contours( gpc_poly ); + calc_points_inside( gpc_poly ); #if 0 diff --git a/src/Lib/Geometry/poly_support.cxx b/src/Lib/Geometry/poly_support.cxx index 3b2dfb32..069e1cb9 100644 --- a/src/Lib/Geometry/poly_support.cxx +++ b/src/Lib/Geometry/poly_support.cxx @@ -22,6 +22,8 @@ // $Id$ +#include + #include #include #include @@ -888,3 +890,261 @@ void calc_points_inside( FGPolygon& p ) { // contour/hole calc_point_inside( ct, p ); } + + +// remove duplicate nodes in a polygon should they exist. Returns the +// fixed polygon +FGPolygon remove_dups( const FGPolygon &poly ) { + FGPolygon result; + result.erase(); + + for ( int i = 0; i < poly.contours(); ++i ) { + Point3D last = poly.get_pt( i, poly.contour_size(i) - 1 ); + bool all_same = true; + for ( int j = 0; j < poly.contour_size(i); ++j ) { + // cout << " " << i << " " << j << endl; + Point3D cur = poly.get_pt( i, j ); + if ( cur == last ) { + // skip + // cout << "skipping a duplicate point" << endl; + } else { + result.add_node( i, cur ); + all_same = false; + last = cur; + } + } + + // make sure the last point doesn't equal the previous or the first. + Point3D begin = poly.get_pt( i, 0 ); + Point3D end = poly.get_pt( i, poly.contour_size(i) - 1 ); + if ( begin == end ) { + // skip + cout << "begin == end!" << endl; + // exit(-1); + } + + if ( !all_same ) { + int flag = poly.get_hole_flag( i ); + result.set_hole_flag( i, flag ); + } else { + // too small an area ... add a token point to the contour + // to keep other things happy + result.add_node( i, begin ); + } + } + + return result; +} + + +static const double tgAirportEpsilon = FG_EPSILON / 10.0; +// static const double tgAirportEpsilon = FG_EPSILON * 10.0; + + +// Find a point in the given node list that lies between start and +// end, return true if something found, false if nothing found. +bool find_intermediate_node( const Point3D& start, const Point3D& end, + const point_list& nodes, Point3D *result ) +{ + bool found_node = false; + double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min; + + Point3D p0 = start; + Point3D p1 = end; + + // cout << " add_intermediate_nodes()" << endl; + printf(" %.7f %.7f %.7f <=> %.7f %.7f %.7f\n", + p0.x(), p0.y(), p0.z(), p1.x(), p1.y(), p1.z() ); + + double xdist = fabs(p0.x() - p1.x()); + double ydist = fabs(p0.y() - p1.y()); + // cout << "xdist = " << xdist << " ydist = " << ydist << endl; + x_err_min = xdist + 1.0; + y_err_min = ydist + 1.0; + + if ( xdist > ydist ) { + // cout << "use y = mx + b" << endl; + + // sort these in a sensible order + Point3D p_min, p_max; + if ( p0.x() < p1.x() ) { + p_min = p0; + p_max = p1; + } else { + p_min = p1; + p_max = p0; + } + + m = (p_min.y() - p_max.y()) / (p_min.x() - p_max.x()); + b = p_max.y() - m * p_max.x(); + + // cout << "m = " << m << " b = " << b << endl; + + for ( int i = 0; i < (int)nodes.size(); ++i ) { + // cout << i << endl; + Point3D current = nodes[i]; + + if ( (current.x() > (p_min.x() + FG_EPSILON)) + && (current.x() < (p_max.x() - FG_EPSILON)) ) { + + // printf( "found a potential candidate %.7f %.7f %.7f\n", + // current.x(), current.y(), current.z() ); + + y_err = fabs(current.y() - (m * current.x() + b)); + // cout << "y_err = " << y_err << endl; + + if ( y_err < tgAirportEpsilon ) { + // cout << "FOUND EXTRA SEGMENT NODE (Y)" << endl; + // cout << p_min << " < " << current << " < " + // << p_max << endl; + found_node = true; + if ( y_err < y_err_min ) { + *result = current; + y_err_min = y_err; + } + } + } + } + } else { + // cout << "use x = m1 * y + b1" << endl; + + // sort these in a sensible order + Point3D p_min, p_max; + if ( p0.y() < p1.y() ) { + p_min = p0; + p_max = p1; + } else { + p_min = p1; + p_max = p0; + } + + m1 = (p_min.x() - p_max.x()) / (p_min.y() - p_max.y()); + b1 = p_max.x() - m1 * p_max.y(); + + // cout << " m1 = " << m1 << " b1 = " << b1 << endl; + // printf( " m = %.8f b = %.8f\n", 1/m1, -b1/m1); + + // cout << " should = 0 = " + // << fabs(p_min.x() - (m1 * p_min.y() + b1)) << endl; + // cout << " should = 0 = " + // << fabs(p_max.x() - (m1 * p_max.y() + b1)) << endl; + + for ( int i = 0; i < (int)nodes.size(); ++i ) { + Point3D current = nodes[i]; + + if ( (current.y() > (p_min.y() + FG_EPSILON)) + && (current.y() < (p_max.y() - FG_EPSILON)) ) { + + // printf( "found a potential candidate %.7f %.7f %.7f\n", + // current.x(), current.y(), current.z() ); + + x_err = fabs(current.x() - (m1 * current.y() + b1)); + // cout << "x_err = " << x_err << endl; + + // if ( temp ) { + // cout << " (" << counter << ") x_err = " << x_err << endl; + // } + + if ( x_err < tgAirportEpsilon ) { + // cout << "FOUND EXTRA SEGMENT NODE (X)" << endl; + // cout << p_min << " < " << current << " < " + // << p_max << endl; + found_node = true; + if ( x_err < x_err_min ) { + *result = current; + x_err_min = x_err; + } + } + } + } + } + + return found_node; +} + + +// Attempt to reduce degeneracies where a subsequent point of a +// polygon lies *on* a previous line segment. These artifacts are +// occasionally introduced by the gpc polygon clipper. +static point_list reduce_contour_degeneracy( const point_list& contour ) { + point_list result = contour; + + Point3D p0, p1, bad_node; + bool done = false; + + while ( !done ) { + // traverse the contour until we find the first bad node or + // hit the end of the contour + cout << " ... not done ... " << endl; + bool bad = false; + + int i = 0; + while ( i < (int)result.size() - 1 && !bad ) { + p0 = result[i]; + p1 = result[i+1]; + + bad = find_intermediate_node( p0, p1, result, &bad_node ); + ++i; + } + if ( !bad ) { + // do the end/start connecting segment + p0 = result[result.size() - 1]; + p1 = result[0]; + + bad = find_intermediate_node( p0, p1, result, &bad_node ); + } + if ( bad ) { + // remove bad node from contour + point_list tmp; tmp.clear(); + for ( int j = 0; j < (int)result.size(); ++j ) { + if ( result[j] == bad_node ) { + // skip + } else { + tmp.push_back( result[j] ); + } + } + result = tmp; + } else { + done = true; + } + } + + return result; +} + +// Search each segment of each contour for degenerate points (i.e. out +// of order points that lie coincident on other segments +FGPolygon reduce_degeneracy( const FGPolygon& poly ) { + FGPolygon result; + + for ( int i = 0; i < poly.contours(); ++i ) { + point_list contour = poly.get_contour(i); + contour = reduce_contour_degeneracy( contour ); + result.add_contour( contour, poly.get_hole_flag(i) ); + + // maintain original hole flag setting + // result.set_hole_flag( i, poly.get_hole_flag( i ) ); + } + + return result; +} + + +// remove any degenerate contours +FGPolygon remove_bad_contours( const FGPolygon &poly ) { + FGPolygon result; + result.erase(); + + for ( int i = 0; i < poly.contours(); ++i ) { + point_list contour = poly.get_contour( i ); + if ( contour.size() >= 3 ) { + // good + int flag = poly.get_hole_flag( i ); + result.add_contour( contour, flag ); + } + } + + return result; +} + + diff --git a/src/Lib/Geometry/poly_support.hxx b/src/Lib/Geometry/poly_support.hxx index e8ac2ebf..bf2ddda7 100644 --- a/src/Lib/Geometry/poly_support.hxx +++ b/src/Lib/Geometry/poly_support.hxx @@ -62,6 +62,24 @@ Point3D calc_point_inside_old( const FGPolygon& p, const int contour, // calculate some "arbitrary" point inside each of the polygons contours void calc_points_inside( FGPolygon& p ); +// remove duplicate nodes in a polygon should they exist. Returns the +// fixed polygon +FGPolygon remove_dups( const FGPolygon &poly ); + + +// Search each segment of each contour for degenerate points (i.e. out +// of order points that lie coincident on other segments +FGPolygon reduce_degeneracy( const FGPolygon& poly ); + + +// Find a point in the given node list that lies between start and +// end, return true if something found, false if nothing found. +bool find_intermediate_node( const Point3D& start, const Point3D& end, + const point_list& nodes, Point3D *result ); + +// remove any degenerate contours +FGPolygon remove_bad_contours( const FGPolygon &poly ); + #endif // _POLY_SUPPORT_HXX diff --git a/src/Lib/Polygon/Makefile.am b/src/Lib/Polygon/Makefile.am index 4164233b..a7d42338 100644 --- a/src/Lib/Polygon/Makefile.am +++ b/src/Lib/Polygon/Makefile.am @@ -4,7 +4,8 @@ libPolygon_a_SOURCES = \ index.cxx index.hxx \ names.cxx names.hxx \ polygon.cxx polygon.hxx \ - split.cxx split.hxx \ + simple_clip.cxx simple_clip.hxx \ + split-bin.cxx split.hxx \ superpoly.cxx superpoly.hxx INCLUDES += -I$(top_srcdir)/src -I$(top_srcdir)/src/Lib diff --git a/src/Lib/Polygon/names.hxx b/src/Lib/Polygon/names.hxx index 7e270bd3..6debc3c8 100644 --- a/src/Lib/Polygon/names.hxx +++ b/src/Lib/Polygon/names.hxx @@ -49,35 +49,36 @@ enum AreaType { GlacierArea = 10, OceanArea = 11, UrbanArea = 12, - MarshArea = 13, // USGS Land Covers // These are low-priority, since known polygons should always win. - BuiltUpCover = 14, // Urban and Built-Up Land - DryCropPastureCover = 15, // Dryland Cropland and Pasture - IrrCropPastureCover = 16, // Irrigated Cropland and Pasture - MixedCropPastureCover = 17, // Mixed Dryland/Irrigated Cropland and Pasture - CropGrassCover = 18, // Cropland/Grassland Mosaic - CropWoodCover = 19, // Cropland/Woodland Mosaic - GrassCover = 20, // Grassland - ShrubCover = 21, // Shrubland - ShrubGrassCover = 22, // Mixed Shrubland/Grassland - SavannaCover = 23, // Savanna - DeciduousBroadCover = 24, // Deciduous Broadleaf Forest - DeciduousNeedleCover = 25, // Deciduous Needleleaf Forest - EvergreenBroadCover = 26, // Evergreen Broadleaf Forest - EvergreenNeedleCover = 27, // Evergreen Needleleaf Forest - MixedForestCover = 28, // Mixed Forest - WaterBodyCover = 29, // Water Bodies - HerbWetlandCover = 30, // Herbaceous Wetland - WoodedWetlandCover = 31, // Wooded Wetland - BarrenCover = 32, // Barren or Sparsely Vegetated - HerbTundraCover = 33, // Herbaceous Tundra - WoodedTundraCover = 34, // Wooded Tundra - MixedTundraCover = 35, // Mixed Tundra - BareTundraCover = 36, // Bare Ground Tundra - SnowCover = 37, // Snow or Ice + BuiltUpCover = 13, // Urban and Built-Up Land + DryCropPastureCover = 14, // Dryland Cropland and Pasture + IrrCropPastureCover = 15, // Irrigated Cropland and Pasture + MixedCropPastureCover = 16, // Mixed Dryland/Irrigated Cropland and Pasture + CropGrassCover = 17, // Cropland/Grassland Mosaic + CropWoodCover = 18, // Cropland/Woodland Mosaic + GrassCover = 19, // Grassland + ShrubCover = 20, // Shrubland + ShrubGrassCover = 21, // Mixed Shrubland/Grassland + SavannaCover = 22, // Savanna + DeciduousBroadCover = 23, // Deciduous Broadleaf Forest + DeciduousNeedleCover = 24, // Deciduous Needleleaf Forest + EvergreenBroadCover = 25, // Evergreen Broadleaf Forest + EvergreenNeedleCover = 26, // Evergreen Needleleaf Forest + MixedForestCover = 27, // Mixed Forest + WaterBodyCover = 28, // Water Bodies + HerbWetlandCover = 29, // Herbaceous Wetland + WoodedWetlandCover = 30, // Wooded Wetland + BarrenCover = 31, // Barren or Sparsely Vegetated + HerbTundraCover = 32, // Herbaceous Tundra + WoodedTundraCover = 33, // Wooded Tundra + MixedTundraCover = 34, // Mixed Tundra + BareTundraCover = 35, // Bare Ground Tundra + SnowCover = 36, // Snow or Ice + + MarshArea = 37, IslandArea = 38, DefaultArea = 39, diff --git a/src/Lib/Polygon/simple_clip.cxx b/src/Lib/Polygon/simple_clip.cxx new file mode 100644 index 00000000..baf6a999 --- /dev/null +++ b/src/Lib/Polygon/simple_clip.cxx @@ -0,0 +1,567 @@ +// simple_clip.cxx -- simplistic polygon clipping routine. Returns +// the portion of a polygon that is above or below +// a horizontal line of y = a. Only works with +// single contour polygons (no holes.) +// +// Written by Curtis Olson, started August 1999. +// +// Copyright (C) 1999 Curtis L. Olson - curt@flightgear.org +// +// 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., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ + + +#include + +#include "simple_clip.hxx" + +FG_USING_STD(cout); +FG_USING_STD(endl); + +#define CLIP_EPSILON 0.000000000001 + + +// Given a line segment specified by two endpoints p1 and p2, return +// the x value of a point on the line that intersects with the +// horizontal line through y. Return true if an intersection is found, +// false otherwise. +static bool intersects_y( Point3D p0, Point3D p1, double y, Point3D *result ) { + // sort the end points + if ( p0.y() > p1.y() ) { + Point3D tmp = p0; + p0 = p1; + p1 = tmp; + } + + if ( (y < p0.y()) || (y > p1.y()) ) { + // out of range of line segment, bail right away + return false; + } + + // equation of a line through (x0,y0) and (x1,y1): + // + // y = y1 + (x - x1) * (y0 - y1) / (x0 - x1) + // x = x1 + (y - y1) * (x0 - x1) / (y0 - y1) + + double x; + + if ( fabs(p0.y() - p1.y()) > CLIP_EPSILON ) { + x = p1.x() + (y - p1.y()) * (p0.x() - p1.x()) / (p0.y() - p1.y()); + } else { + return false; + } + result->setx(x); + result->sety(y); + + if ( p0.x() <= p1.x() ) { + if ( (p0.x() <= x) && (x <= p1.x()) ) { + return true; + } + } else { + if ( (p0.x() >= x) && (x >= p1.x()) ) { + return true; + } + } + + return false; +} + + +// find the smallest point in contour 0 of poly such that x > min_x +// and y = y. Returns index of the found point, -1 if no match found. +static int find_point( const FGPolygon& poly, double min_x, double y ) { + Point3D p, save; + int index = -1; + + save.setx( 361.0 ); + + for ( int i = 0; i < poly.contour_size( 0 ); ++i ) { + p = poly.get_pt( 0, i ); + if ( p.y() == y ) { + // printf("(%d) p.y() = %.12f y = %.12f\n", i, p.y(), y); + // cout << " " << p << endl; + if ( p.x() > min_x ) { + if ( p.x() < save.x() ) { + save = p; + index = i; + } + } + } + } + + return index; +} + + +// return if interesection is valid (not in the ignore list) +static bool valid_intersection( int intersection, const int_list& ignore_ints ) +{ + for ( int i = 0; i < (int)ignore_ints.size(); ++i ) { + if ( intersection == ignore_ints[i] ) { + return false; + } + } + return true; +} + + +// return index of next valid intersection +static int next_intersection( const int_list& keep_ints, + const int_list& ignore_ints, + const int beginning_at ) +{ + // cout << "[ni] start_int = " << beginning_at << endl; + int i = beginning_at; + if ( i < 0 ) { i = 0; } + while ( i < (int)keep_ints.size() ) { + // cout << " i = " << i << endl; + if ( keep_ints[i] != -1 ) { + if ( valid_intersection(keep_ints[i], ignore_ints) ) { + return i; + } + } + ++i; + } + + return -1; +} + + +// return true if point.y() touches or is inside of line, else return false +inline bool is_on_or_inside( double line, Point3D p, fgSideType side ) { + if ( side == Above ) { + if ( p.y() >= line ) { + return true; + } + } else if ( side == Below ) { + if ( p.y() <= line ) { + return true; + } + } + + return false; +} + + +// return true if point.y() is inside of line, else return false +inline bool is_inside( double line, Point3D p, fgSideType side ) { + if ( side == Above ) { + if ( p.y() > line ) { + return true; + } + } else if ( side == Below ) { + if ( p.y() < line ) { + return true; + } + } + + return false; +} + + +// Walk through the input polygon and split it into the +// portion that is inside the clip region +static bool simple_clip( const FGPolygon& in, const double y, + const fgSideType side, + FGPolygon& result ) +{ + Point3D p, p_last, p_int; + int i; + + result.erase(); + + cout << "input poly size = " << in.total_size() << endl; + + p_last = in.get_pt( 0, in.contour_size(0)-1 ); + + for ( i = 0; i < (int)in.contour_size(0); ++i ) { + p = in.get_pt( 0, i ); + + if ( (fabs(p.x() - p_last.x()) < CLIP_EPSILON) && + (fabs(p.y() - p_last.y()) < CLIP_EPSILON) && + (i > 0) ) { + // cout << "WARNING: p and p_last are identical at index = " + // << i << endl; + } + + if ( is_on_or_inside(y, p, side) ) { + if ( is_on_or_inside(y, p_last, side) ) { + // cout << "inside & inside " << i << " " << p << endl; + result.add_node( 0, p ); + } else { + if ( !intersects_y(p, p_last, y, &p_int) ) { + cout << "Huh, this should have intersected!" << endl; + return false; + } else { + // cout << "intersection outside to inside " << i << " " + // << p_int << endl; + // cout << " i - 1 = " << in.get_pt( 0, i-1 ) << endl; + // cout << " i = " << in.get_pt( 0, i ) << endl; + // cout << " i + 1 = " << in.get_pt( 0, i+1 ) << endl; + result.add_node( 0, p_int ); + if ( (fabs(p.x() - p_int.x()) < CLIP_EPSILON) && + (fabs(p.y() - p_int.y()) < CLIP_EPSILON) ) + { + // cout << "WARNING: p and p_int are identical, "; + // cout << "omitting p" << endl; + } else { + cout << "adding intersection" << i << " " << p << endl; + result.add_node( 0, p ); + } + } + } + } else { + if ( is_inside(y, p_last, side) ) { + if ( !intersects_y(p, p_last, y, &p_int) ) { + cout << "Huh, this should have intersected!" << endl; + return false; + } else { + // cout << "intersection inside to outside " << i << " " + // << p_int << endl; + if ( (fabs(p.x() - p_int.x()) < CLIP_EPSILON) && + (fabs(p.y() - p_int.y()) < CLIP_EPSILON) ) + { + cout << "WARNING: p and p_int are identical, "; + cout << "omitting p" << endl; + } else { + result.add_node( 0, p_int ); + } + } + } + } + + p_last = p; + } + + return true; +} + + +// build the list of intersections +static bool build_intersections( const FGPolygon& arcs, double line, + fgSideType side, + int_list& keep_ints, + int_list& ignore_ints ) +{ + keep_ints.clear(); + ignore_ints.clear(); + + int index = 0; + double current_x = -181.0; + + while ( index >= 0 ) { + index = find_point( arcs, current_x, line ); + if ( index >= 0 ) { + cout << "intersection at " << index << " = " + << arcs.get_pt( 0, index ) << endl; + keep_ints.push_back( index ); + current_x = arcs.get_pt( 0, index ).x(); + + int before = index - 1; + if ( before < 0 ) { before += arcs.contour_size(0); } + int after = (index + 1) % arcs.contour_size(0); + cout << endl; + cout << " before = " << arcs.get_pt(0, before) << endl; + cout << " int = " << arcs.get_pt(0, index) << endl; + cout << " after = " << arcs.get_pt(0, after) << endl; + if ( side == Above ) { + if ( (arcs.get_pt(0, before).y() > line) && + (arcs.get_pt(0, after).y() > line) ) + { + cout << "side = above" << endl; + cout << "V intersection with clip line from above" << endl; + cout << "Adding intersection to ignore_ints" << endl; + ignore_ints.push_back( index ); + } + if ( (arcs.get_pt(0, before).y() <= line) && + (arcs.get_pt(0, after).y() <= line) ) + { + cout << "side = above" << endl; + cout << "V intersection with clip line from BELOW" << endl; + cout << "or an extra in-clip-line intersection." << endl; + cout << "Adding intersection to ignore_ints" << endl; + ignore_ints.push_back( index ); + } + } else if ( side == Below ) { + if ( (arcs.get_pt(0, before).y() >= line) && + (arcs.get_pt(0, after).y() >= line) ) + { + cout << "side = below" << endl; + cout << "V intersection with clip line from above" << endl; + cout << "Adding intersection to ignore_ints" << endl; + ignore_ints.push_back( index ); + } + if ( (arcs.get_pt(0, before).y() < line) && + (arcs.get_pt(0, after).y() < line) ) + { + cout << "side = below" << endl; + cout << "V intersection with clip line from BELOW" << endl; + cout << "or an extra in-clip-line intersection." << endl; + cout << "Adding intersection to ignore_ints" << endl; + ignore_ints.push_back( index ); + } + } + } + } + + return true; +} + + +// test for duplicate nodes +FGPolygon fix_dups( FGPolygon& in ) { + FGPolygon result; + + double x_last = -20000.0; + double y_last = -20000.0; + double x, y; + + for ( int i = 0; i < (int)in.contour_size(0); ++i ) { + x = in.get_pt(0, i).x(); + y = in.get_pt(0, i).y(); + if ( (x == x_last) && (y == y_last) ) { + // ignore + } else { + result.add_node(0, in.get_pt(0, i)); + } + x_last = x; + y_last = y; + } + + return result; +} + + +// simple polygon clipping routine. Returns the portion of a polygon +// that is above and below a horizontal line of y = a. Only works +// with single contour polygons (no holes.) Returns true if routine +// thinks it was successful. + +static bool clip_contour( const FGPolygon& in, const double y, + const fgSideType side, + FGPolygon& result ) +{ + FGPolygon result_arcs, arcs; + int i, i1, i2, index; + + + // Step 1: sanity checks + if ( (int)in.contours() != 1 ) { + cout << "we only handle single contour polygons" << endl; + return false; + } + + if ( (int)in.contour_size( 0 ) < 3 ) { + cout << "we must have at least three vertices to work" << endl; + return false; + } + + + // Step 2: walk through the input polygon and split it into the + // portion that is on or inside the clip line + + if ( simple_clip( in, y, side, result_arcs ) ) { + if ( result_arcs.contours() > 0 ) { + cout << "result_arcs size = " + << result_arcs.total_size() << endl; + } else { + cout << "empty result" << endl; + } + } else { + cout << "simple_clip_above() failed!" << endl; + exit(-1); + } + + + // Step 3: check for trivial cases + + result.erase(); + + // trivial -- nothing inside of clip line + if ( result_arcs.contours() == 0 ) { + cout << "trivially empty" << endl; + return true; + } + + // trivial -- everything inside of clip line + i1 = find_point( result_arcs, -181.0, y ); + if ( i1 < 0 ) { + cout << "trivially full" << endl; + result = result_arcs; + return true; + } + + // trivial -- single clip line intersection (polygon just nicks + // it) -- everything inside + i2 = find_point( result_arcs, result_arcs.get_pt(0,i1).x(), y ); + if ( i2 < 0 ) { + cout << "trivially full (clip line nicks edge)" << endl; + result = result_arcs; + return true; + } + + + // Step 4: If we are finding the "below" clip, reverse the points + // before extracting the cycles. (and remove duplicates) + + FGPolygon tmp; + tmp.erase(); + + if ( side == Below ) { + for ( i = result_arcs.contour_size(0) - 1; i >= 0; --i ) { + Point3D p = result_arcs.get_pt( 0, i ); + tmp.add_node( 0, p ); + } + } else { + tmp = result_arcs; + } + + arcs = fix_dups( tmp ); + + // Step 5: Build the intersections lists + + int_list keep_ints, ignore_ints; + build_intersections( arcs, y, side, keep_ints, ignore_ints ); + cout << "total keep_ints = " << keep_ints.size() << endl; + cout << "total ignore_ints = " << ignore_ints.size() << endl; + + + // Step 6: Walk through the result_arcs and extract the cycles (or + // individual contours.) + + int start_int = next_intersection( keep_ints, ignore_ints, 0 ); + int next_int = next_intersection( keep_ints, ignore_ints, start_int+1 ); + cout << "start_int = " << start_int << endl; + cout << "next_int = " << next_int << endl; + + int count = 0; + + // while we have keep_ints left to process + while ( start_int >= 0 ) { + point_list contour; + contour.clear(); + + index = keep_ints[next_int]; + keep_ints[next_int] = -1; + cout << endl << "starting at point = " << arcs.get_pt(0,index) << endl; + + while ( index != keep_ints[start_int] ) { + cout << "index = " << index << " start_int = " << start_int + << " keep_ints[start_int] = " << keep_ints[start_int] + << endl; + + // start with the 2nd item in the intersection list and + // traverse until we find another intersection + contour.push_back( arcs.get_pt(0,index) ); + index = (index + 1) % arcs.contour_size(0); + + while ( (arcs.get_pt(0,index).y() != y) || + ! valid_intersection(index, ignore_ints) ) + { + contour.push_back( arcs.get_pt(0,index) ); + index = (index + 1) % arcs.contour_size(0); + } + contour.push_back( arcs.get_pt(0,index) ); + cout << "exited at poly index = " << index << " " + << arcs.get_pt(0,index) << endl; + + // find which intersection we came out on in our list + cout << "finding exit intersection for " << index << endl; + i = 0; + while ( i < (int)keep_ints.size() ) { + // cout << " keep_int[" << i << "] = " << keep_ints[i] << endl; + if ( index == keep_ints[i] ) { + cout << " intersection index = " << i << endl; + if ( index != keep_ints[start_int] ) { + cout << " not start index so keep going" << endl; + keep_ints[i] = -1; + next_int = next_intersection( keep_ints, ignore_ints, + i+1 ); + index = keep_ints[next_int]; + keep_ints[next_int] = -1; + cout << " next_int = " << next_int << " index = " + << index << endl; + } + break; + } + ++i; + } + if ( i == (int)keep_ints.size() ) { + cout << "oops, didn't find that intersection, you are screwed" + << endl; + exit(-1); + } + } + keep_ints[start_int] = -1; + result.add_contour( contour, count ); + count++; + + // find next keep_ints + start_int = next_intersection( keep_ints, ignore_ints, -1 ); + next_int = next_intersection( keep_ints, ignore_ints, start_int+1 ); + + cout << "start_int = " << start_int << endl; + cout << "next_int = " << next_int << endl; + } + + return true; +} + + +// simple polygon clipping routine. Returns the portion of a polygon +// that is above and below a horizontal line of y = a. Clips +// multi-contour polygons individually and then reassembles the +// results. Doesn't work with holes. Returns true if routine thinks +// it was successful. + +FGPolygon horizontal_clip( const FGPolygon& in, const double y, + const fgSideType side ) +{ + FGPolygon result; + result.erase(); + + // Step 1: sanity checks + if ( (int)in.contours() == 0 ) { + cout << "Error: 0 contour polygon" << endl; + return result; + } + + // clip each contour individually + FGPolygon tmp_poly, clip_poly; + point_list contour; + + for ( int i = 0; i < in.contours(); ++i ) { + if ( (int)in.contour_size( i ) < 3 ) { + cout << "we must have at least three vertices to work" << endl; + return result; + } + + tmp_poly.erase(); + + contour = in.get_contour( i ); + tmp_poly.add_contour( contour, 0 ); + + clip_contour( tmp_poly, y, side, clip_poly ); + + // add each clip_poly contour to the result poly + for ( int j = 0; j < clip_poly.contours(); ++j ) { + contour = clip_poly.get_contour( j ); + result.add_contour( contour, 0 ); + } + } + + return result; +} diff --git a/src/Lib/Polygon/simple_clip.hxx b/src/Lib/Polygon/simple_clip.hxx new file mode 100644 index 00000000..2934158f --- /dev/null +++ b/src/Lib/Polygon/simple_clip.hxx @@ -0,0 +1,53 @@ +// simple_clip.hxx -- simplistic polygon clipping routine. Returns +// the portion of a polygon that is above or below +// a horizontal line of y = a. Only works with +// single contour polygons (no holes.) +// +// Written by Curtis Olson, started August 1999. +// +// Copyright (C) 1999 Curtis L. Olson - curt@flightgear.org +// +// 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., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ + + +#ifndef _SIMPLE_CLIP_HXX +#define _SIMPLE_CLIP_HXX + + +#include "polygon.hxx" + + +// side type +enum fgSideType { + Above = 0, + Below = 1 +}; + + +// simple polygon clipping routine. Returns the portion of a polygon +// that is above and below a horizontal line of y = a. Clips +// multi-contour polygons individually and then reassembles the +// results. Doesn't work with holes. Returns true if routine thinks +// it was successful. + +FGPolygon horizontal_clip( const FGPolygon& in, const double y, + const fgSideType side ); + + +#endif // _SIMPLE_CLIP_HXX + + diff --git a/src/Lib/Polygon/split.cxx b/src/Lib/Polygon/split.cxx index 646fc4bb..cc17a51f 100644 --- a/src/Lib/Polygon/split.cxx +++ b/src/Lib/Polygon/split.cxx @@ -181,7 +181,7 @@ void split_polygon(const string& path, AreaType area, const FGPolygon& shape) { if ( (dx > 2880) || (dy > 1440) ) { FG_LOG( FG_GENERAL, FG_ALERT, - "somethings really wrong!!!!" ); + "something is really wrong in split_polygon()!!!!" ); exit(-1); } @@ -238,6 +238,7 @@ void split_polygon(const string& path, AreaType area, const FGPolygon& shape) { b_cur = fgBucketOffset(min.x(), min.y(), i, j); clip_and_write_poly( path, index, area, b_cur, clip_row ); } + cout << " (done)" << endl; } // string answer; cin >> answer; } diff --git a/src/Prep/GSHHS/Makefile.am b/src/Prep/GSHHS/Makefile.am index 49d440ba..f2fd5591 100644 --- a/src/Prep/GSHHS/Makefile.am +++ b/src/Prep/GSHHS/Makefile.am @@ -3,8 +3,7 @@ bin_PROGRAMS = gshhs debug gshhs_SOURCES = \ main.cxx \ gshhs.h \ - gshhs_split.cxx gshhs_split.hxx \ - simple_clip.hxx simple_clip.cxx + gshhs_split.cxx gshhs_split.hxx gshhs_LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ diff --git a/src/Prep/GSHHS/gshhs_split.cxx b/src/Prep/GSHHS/gshhs_split.cxx index ca7fc312..afe4b0db 100644 --- a/src/Prep/GSHHS/gshhs_split.cxx +++ b/src/Prep/GSHHS/gshhs_split.cxx @@ -32,9 +32,9 @@ #include #include #include +#include #include "gshhs_split.hxx" -#include "simple_clip.hxx" FG_USING_STD(cout); FG_USING_STD(string); diff --git a/src/Prep/ShapeFile/Makefile.am b/src/Prep/ShapeFile/Makefile.am index f301389b..da613c52 100644 --- a/src/Prep/ShapeFile/Makefile.am +++ b/src/Prep/ShapeFile/Makefile.am @@ -1,6 +1,8 @@ -bin_PROGRAMS = shape-decode +bin_PROGRAMS = shape-decode noaa-decode -shape_decode_SOURCES = main.cxx +shape_decode_SOURCES = shape-decode.cxx + +noaa_decode_SOURCES = noaa-decode.cxx shape_decode_LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ @@ -8,4 +10,10 @@ shape_decode_LDADD = \ $(top_builddir)/src/Lib/shapelib/libshape.a \ -lsgdebug -lsgbucket -lsgmisc -lz -lgpc +noaa_decode_LDADD = \ + $(top_builddir)/src/Lib/Polygon/libPolygon.a \ + $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ + $(top_builddir)/src/Lib/shapelib/libshape.a \ + -lsgdebug -lsgbucket -lsgmisc -lz -lgpc + INCLUDES += -I$(top_srcdir)/src/Lib diff --git a/src/Prep/ShapeFile/noaa-decode.cxx b/src/Prep/ShapeFile/noaa-decode.cxx new file mode 100644 index 00000000..d5edded1 --- /dev/null +++ b/src/Prep/ShapeFile/noaa-decode.cxx @@ -0,0 +1,461 @@ +// main.cxx -- process shapefiles and extract polygon outlines, +// clipping against and sorting them into the revelant +// tiles. +// +// Written by Curtis Olson, started February 1999. +// +// Copyright (C) 1999 Curtis L. Olson - curt@flightgear.org +// +// 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., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ + + +#include + +#include STL_STRING + +#include + +#include +#include +#include +#include +#include + +#ifdef _MSC_VER +# include +#endif + + +// return the type of the shapefile record +AreaType get_shapefile_type(DBFHandle& hDBF, int rec) { + int *panWidth, i, iRecord; + char szFormat[32]; + int nWidth, nDecimals; + int bMultiLine = 0; + char szTitle[12]; + +#if 0 + // grab the meta-information for all the fields + // this applies to all the records in the DBF file. + for ( i = 0; i < DBFGetFieldCount(hDBF); i++ ) { + DBFFieldType eType; + string pszTypeName; + + eType = DBFGetFieldInfo( hDBF, i, szTitle, &nWidth, &nDecimals ); + if( eType == FTString ) + pszTypeName = "String"; + else if( eType == FTInteger ) + pszTypeName = "Integer"; + else if( eType == FTDouble ) + pszTypeName = "Double"; + else if( eType == FTInvalid ) + pszTypeName = "Invalid"; + + // printf( "Field %d: Type=%s, Title=`%s', Width=%d, Decimals=%d\n", + // i, pszTypeName.c_str(), szTitle, nWidth, nDecimals ); + } + + // Compute offsets to use when printing each of the field + // values. We make each field as wide as the field title+1, or the + // field value + 1. + + panWidth = (int *) malloc( DBFGetFieldCount( hDBF ) * sizeof(int) ); + for ( i = 0; i < DBFGetFieldCount(hDBF) && !bMultiLine; i++ ) { + DBFFieldType eType; + + eType = DBFGetFieldInfo( hDBF, i, szTitle, &nWidth, &nDecimals ); + if( (int)strlen(szTitle) > nWidth ) { + panWidth[i] = strlen(szTitle); + } else { + panWidth[i] = nWidth; + } + + if( eType == FTString ) { + sprintf( szFormat, "%%-%ds ", panWidth[i] ); + } else { + sprintf( szFormat, "%%%ds ", panWidth[i] ); + } + printf( szFormat, szTitle ); + } + printf( "\n" ); + + for( iRecord = 0; iRecord < DBFGetRecordCount(hDBF); iRecord++ ) { + for( i = 0; i < DBFGetFieldCount(hDBF); i++ ) { + DBFFieldType eType; + + eType = DBFGetFieldInfo( hDBF, i, szTitle, &nWidth, &nDecimals ); + + switch( eType ) + { + case FTString: + sprintf( szFormat, "%%-%ds", nWidth ); + printf( szFormat, + DBFReadStringAttribute( hDBF, iRecord, i ) ); + break; + + case FTInteger: + sprintf( szFormat, "%%%dd", nWidth ); + printf( szFormat, + DBFReadIntegerAttribute( hDBF, iRecord, i ) ); + break; + + case FTDouble: + sprintf( szFormat, "%%%d.%dlf", nWidth, nDecimals ); + printf( szFormat, + DBFReadDoubleAttribute( hDBF, iRecord, i ) ); + break; + } + } + } + printf( "\n" ); + +#endif + + string area = DBFReadStringAttribute( hDBF, rec, 2 ); + string code = DBFReadStringAttribute( hDBF, rec, 1 ); + // cout << "next record = " << test << endl; + + // strip leading spaces + while ( area[0] == ' ' ) { + area = area.substr(1, area.length() - 1); + } + // strip trailing spaces + while ( area[area.length() - 1] == ' ' ) { + area = area.substr(0, area.length() - 1); + } + // strip other junk encountered + while ( (int)area[area.length() - 1] == 9 ) { + area = area.substr(0, area.length() - 1); + } + + FG_LOG( FG_GENERAL, FG_INFO, " raw area = " << area ); + FG_LOG( FG_GENERAL, FG_INFO, " code = " << code ); + + // there may be a better way to handle various area types, but + // this is what we are trying for starters + if ( area == "Urban (1990 Enhanced)" ) { + area = "Urban"; + } else if ( area == "Residential" ) { + area = "Urban"; + } else if ( area == "Commercial and Services" ) { + area = "Urban"; + } else if ( area == "Industrial" ) { + area = "Urban"; + } else if ( area == "Transportation, Communications and Utilities" ) { + area = "Urban"; + } else if ( area == "Transportation, Communications, and Utilities" ) { + area = "Urban"; + } else if ( area == "Industrial and Commercial Complexes" ) { + area = "Urban"; + } else if ( area == "Mixed Urban or Built-up Land" ) { + area = "BuiltUpCover"; + } else if ( area == "Mixed Urban or Built up Land" ) { + area = "BuiltUpCover"; + } else if ( area == "Other Urban or Built-up Land" ) { + area = "BuiltUpCover"; + } else if ( area == "Other Urban or Built up Land" ) { + area = "BuiltUpCover"; + + } else if ( area == "Cropland and Pasture" ) { + area = "MixedCropPastureCover"; + } else if ( area == "Orchards, Groves, Vineyards, Nurseries, and Ornamental Horticultural Areas" ) { + area = "IrrCropPastureCover"; + } else if ( area == "Orchards, Groves, Vineyards, Nurseries and Ornament" ) { + area = "IrrCropPastureCover"; + } else if ( area == "Orchards, Groves, Vineyards, Nurseries and Ornamental Hort" ) { + area = "IrrCropPastureCover"; + } else if ( area == "Confined Feeding Operations" ) { + area = "MixedCropPastureCover"; + } else if ( area == "Other Agricultural Land" ) { + area = "IrrCropPastureCover"; + + } else if ( area == "Herbaceous Rangeland" ) { + area = "GrassCover"; + } else if ( area == "Shrub and Brush Rangeland" ) { + area = "ShrubCover"; + } else if ( area == "Mixed Rangeland" ) { + area = "ShrubGrassCover"; + + } else if ( area == "Deciduous Forest Land" ) { + area = "DeciduousBroadCover"; + } else if ( area == "Evergreen Forest Land" ) { + area = "EvergreenNeedleCover"; + } else if ( area == "Mixed Forest Land" ) { + area = "MixedForestCover"; + + } else if ( area == "Streams and Canals" ) { + area = "Stream"; + } else if ( area == "Lakes" ) { + area = "Lake"; + } else if ( area == "Reservoirs" ) { + area = "Reservoir"; + } else if ( area == "Bays and Estuaries" ) { + area = "Ocean"; + + } else if ( area == "Forested Wetland" ) { + area = "WoodedWetlandCover"; + } else if ( area == "Nonforested Wetland" ) { + area = "HerbWetlandCover"; + + } else if ( area == "Dry Salt Flats" ) { + area = "DryLake"; + } else if ( area == "Beaches" ) { + area = "DryLake"; + } else if ( area == "Sandy Areas Other than Beaches" ) { + area = "DryLake"; + } else if ( area == "Bare Exposed Rock" ) { + area = "BarrenCover"; + } else if ( area == "Strip Mines, Quarries, and Gravel Pits" ) { + area = "BarrenCover"; + } else if ( area == "Transitional Areas" ) { + area = "BarrenCover"; + } else if ( area == "Mixed Barren Land" ) { + area = "BarrenCover"; + + } else if ( area == "Shrub and Brush Tundra" ) { + area = "WoodedTundraCover"; + } else if ( area == "Herbaceous Tundra" ) { + area = "HerbTundraCover"; + } else if ( area == "Bare Ground" ) { + area = "BareTundraCover"; + } else if ( area == "Wet Tundra" ) { + area = "HerbTundraCover"; + } else if ( area == "Mixed Tundra" ) { + area = "MixedTundraCover"; + + } else if ( area == "Perennial Snowfields" ) { + area = "Glacier"; + } else if ( area == "Glaciers" ) { + area = "Glacier"; + + } else if ( area == "No Data" ) { + area = "Null"; + } + + return get_area_type( area ); +} + + +int main( int argc, char **argv ) { + FGPolygon shape; + int i, j; + + fglog().setLogLevels( FG_ALL, FG_DEBUG ); + + if ( argc < 3 ) { + FG_LOG( FG_GENERAL, FG_ALERT, "Usage: " << argv[0] + << " [ area_string ]" ); + exit(-1); + } + + FG_LOG( FG_GENERAL, FG_DEBUG, "Opening " << argv[1] << " for reading." ); + + // make work directory + string work_dir = argv[2]; + +#ifdef _MSC_VER + fg_mkdir( work_dir.c_str() ); +#else + string command = "mkdir -p " + work_dir; + system( command.c_str() ); +#endif + + // allow us to override the area type from the command line. All + // polygons in the processed shape file will be assigned this area + // type + string force_area_type = ""; + if ( argc == 4 ) { + force_area_type = argv[3]; + } + + // initialize persistant polygon counter + string counter_file = work_dir + "/../poly_counter"; + poly_index_init( counter_file ); + + string dbffile = argv[1]; + dbffile += ".dbf"; + DBFHandle hDBF = DBFOpen( dbffile.c_str(), "rb" ); + if( hDBF == NULL ) { + FG_LOG( FG_GENERAL, FG_ALERT, "DBFOpen(" << dbffile + << ",\"rb\") failed." ); + exit( -1 ); + } + + string shpfile = argv[1]; + shpfile += ".shp"; + SHPHandle hSHP = SHPOpen( shpfile.c_str(), "rb" ); + if( hSHP == NULL ) { + FG_LOG( FG_GENERAL, FG_ALERT, "SHPOpen(" << shpfile + << ",\"rb\") failed." ); + exit( -1 ); + } + + int nShapeType, nEntities; + double adfMinBound[4], adfMaxBound[4]; + SHPGetInfo( hSHP, &nEntities, &nShapeType, adfMinBound, adfMaxBound ); + + FG_LOG( FG_GENERAL, FG_INFO, "shape file records = " << nEntities << endl ); + + string shapetype = SHPTypeName( nShapeType ); + + if ( shapetype != "Polygon" ) { + FG_LOG( FG_GENERAL, FG_ALERT, "Can't handle non-polygon shape files" ); + exit(-1); + } + + int iPart; + const char *pszPlus; + for ( i = 0; i < nEntities; i++ ) { + // fetch i-th record (shape) + SHPObject *psShape; + + shape.erase(); + + psShape = SHPReadObject( hSHP, i ); + + FG_LOG( FG_GENERAL, FG_DEBUG, "Processing record = " << i + << " of " << nEntities + << " rings = " << psShape->nParts + << " total vertices = " << psShape->nVertices ); + + AreaType area = DefaultArea; + if ( force_area_type.length() == 0 ) { + area = get_shapefile_type(hDBF, i); + FG_LOG( FG_GENERAL, FG_DEBUG, " area type = " + << get_area_name(area) << " (" << (int)area << ")" ); + } else { + area = get_area_type( force_area_type ); + } + + FG_LOG( FG_GENERAL, FG_INFO, " record type = " + << SHPTypeName(psShape->nSHPType) ); + FG_LOG( FG_GENERAL, FG_INFO, " bounds = (" + << psShape->dfXMin << "," << psShape->dfYMin << ") " + << psShape->dfZMin << "," << psShape->dfMMin + << " to (" << psShape->dfXMax << "," << psShape->dfYMax << ") " + << psShape->dfZMax << "," << psShape->dfMMax ); + +#if 0 + printf( "\nShape:%d (%s) nVertices=%d, nParts=%d\n" + " Bounds:(%12.3f,%12.3f, %g, %g)\n" + " to (%12.3f,%12.3f, %g, %g)\n", + i, SHPTypeName(psShape->nSHPType), + psShape->nVertices, psShape->nParts, + psShape->dfXMin, psShape->dfYMin, + psShape->dfZMin, psShape->dfMMin, + psShape->dfXMax, psShape->dfYMax, + psShape->dfZMax, psShape->dfMMax ); +#endif + + for ( j = 0, iPart = 1; j < psShape->nVertices; j++ ) { + const char *pszPartType = ""; + + if ( j == 0 && psShape->nParts > 0 ) { + pszPartType = SHPPartTypeName( psShape->panPartType[0] ); + } + + if( iPart < psShape->nParts + && psShape->panPartStart[iPart] == j ) + { + pszPartType = SHPPartTypeName( psShape->panPartType[iPart] ); + iPart++; + pszPlus = "+"; + } else { + pszPlus = " "; + } + + shape.add_node( iPart - 1, + Point3D(psShape->padfX[j], psShape->padfY[j], 0) + ); +#if 0 + printf("%d %d %s (%12.3f,%12.3f, %g, %g) %s \n", + iPart, j, + pszPlus, + psShape->padfX[j], + psShape->padfY[j], + psShape->padfZ[j], + psShape->padfM[j], + pszPartType ); +#endif + } + + SHPDestroyObject( psShape ); + + // check/set hole status for each contour. negative area + // means counter clockwise winding indicating the ring/contour + // is a hole. + for ( int i = 0; i < shape.contours(); ++i ) { + double area = shape.area_contour( i ); + if ( area > 0 ) { + cout << "contour " << i << " = area" << endl; + shape.set_hole_flag( i, false ); + } else { + cout << "contour " << i << " = hole" << endl; + shape.set_hole_flag( i, true ); + } + } + + if ( force_area_type.length() > 0 ) { + // interior of polygon is assigned to force_area_type, + // holes are preserved + + area = get_area_type( force_area_type ); + split_polygon(work_dir, area, shape); + } else if ( area == OceanArea ) { + // interior of polygon is ocean, holes are islands + + FG_LOG( FG_GENERAL, FG_ALERT, "Ocean area ... SKIPPING!" ); + + // Ocean data now comes from GSHHS so we want to ignore + // all other ocean data + // split_polygon(work_dir, area, shape); + } else if ( area == VoidArea ) { + // interior is ???? + + // skip for now + FG_LOG( FG_GENERAL, FG_ALERT, "Void area ... SKIPPING!" ); + + if ( shape.contours() > 1 ) { + FG_LOG( FG_GENERAL, FG_ALERT, " Void area with holes!" ); + // exit(-1); + } + + // split_polygon(work_dir, area, shape); + } else if ( area == NullArea ) { + // interior is ???? + + // skip for now + FG_LOG( FG_GENERAL, FG_ALERT, "Null area ... SKIPPING!" ); + + if ( shape.contours() > 1 ) { + FG_LOG( FG_GENERAL, FG_ALERT, " Null area with holes!" ); + // exit(-1); + } + + // split_polygon(work_dir, area, shape); + } else { + split_polygon(work_dir, area, shape); + } + } + + DBFClose( hDBF ); + SHPClose( hSHP ); + + return 0; +} + + diff --git a/src/Prep/ShapeFile/main.cxx b/src/Prep/ShapeFile/shape-decode.cxx similarity index 88% rename from src/Prep/ShapeFile/main.cxx rename to src/Prep/ShapeFile/shape-decode.cxx index 852f0394..74f7d57d 100644 --- a/src/Prep/ShapeFile/main.cxx +++ b/src/Prep/ShapeFile/shape-decode.cxx @@ -91,6 +91,7 @@ AreaType get_shapefile_type(DBFHandle& hDBF, int rec) { } printf( szFormat, szTitle ); } + printf( "\n" ); for( iRecord = 0; iRecord < DBFGetRecordCount(hDBF); iRecord++ ) { for( i = 0; i < DBFGetFieldCount(hDBF); i++ ) { @@ -120,9 +121,13 @@ AreaType get_shapefile_type(DBFHandle& hDBF, int rec) { } } } + printf( "\n" ); + #endif string area = DBFReadStringAttribute( hDBF, rec, 4 ); + string test = DBFReadStringAttribute( hDBF, rec, 3 ); + cout << "next record = " << test << endl; // strip leading spaces while ( area[0] == ' ' ) { @@ -292,10 +297,6 @@ int main( int argc, char **argv ) { // holes are preserved area = get_area_type( force_area_type ); - split_polygon(work_dir, area, shape); - } else if ( area == MarshArea ) { - // interior of polygon is marsh, holes are preserved - split_polygon(work_dir, area, shape); } else if ( area == OceanArea ) { // interior of polygon is ocean, holes are islands @@ -305,38 +306,6 @@ int main( int argc, char **argv ) { // Ocean data now comes from GSHHS so we want to ignore // all other ocean data // split_polygon(work_dir, area, shape); - } else if ( area == LakeArea ) { - // interior of polygon is lake, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == DryLakeArea ) { - // interior of polygon is dry lake, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == IntLakeArea ) { - // interior of polygon is intermittent lake, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == ReservoirArea ) { - // interior of polygon is reservoir, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == IntReservoirArea ) { - // interior of polygon is intermittent reservoir, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == StreamArea ) { - // interior of polygon is stream, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == CanalArea ) { - // interior of polygon is canal, holes are islands - - split_polygon(work_dir, area, shape); - } else if ( area == GlacierArea ) { - // interior of polygon is glacier, holes are dry land - - split_polygon(work_dir, area, shape); } else if ( area == VoidArea ) { // interior is ???? @@ -362,8 +331,7 @@ int main( int argc, char **argv ) { // split_polygon(work_dir, area, shape); } else { - FG_LOG( FG_GENERAL, FG_ALERT, "Uknown area!" ); - exit(-1); + split_polygon(work_dir, area, shape); } }