diff --git a/src/Lib/Polygon/chop-bin.cxx b/src/Lib/Polygon/chop-bin.cxx index 9f93870a..cbc2a955 100644 --- a/src/Lib/Polygon/chop-bin.cxx +++ b/src/Lib/Polygon/chop-bin.cxx @@ -39,6 +39,8 @@ # include #endif +#include + #include "index.hxx" #include "names.hxx" #include "simple_clip.hxx" @@ -94,7 +96,25 @@ static void clip_and_write_poly( string root, long int p_index, AreaType area, fclose(bfp); */ + cout << "shape contours = " << shape.contours() << " "; + for ( int ii = 0; ii < shape.contours(); ii++ ) { + cout << "hole = " << shape.get_hole_flag(ii) << " "; + } + cout << endl; + result = tgPolygonInt( base, shape ); + + // write_polygon(shape, "shape"); + // write_polygon(result, "result"); + + cout << "result contours = " << result.contours() << " "; + for ( int ii = 0; ii < result.contours(); ii++ ) { + cout << "hole = " << result.get_hole_flag(ii) << " "; + //if ( result.get_hole_flag(ii) ) { + // exit(0); + //} + } + cout << endl; if ( preserve3d ) { result.inherit_elevations( shape ); } diff --git a/src/Prep/E00Lines/Makefile.am b/src/Prep/E00Lines/Makefile.am index 5a3e0d28..7e8392b0 100644 --- a/src/Prep/E00Lines/Makefile.am +++ b/src/Prep/E00Lines/Makefile.am @@ -5,6 +5,7 @@ e00lines_SOURCES = main.cxx e00lines_LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ $(top_builddir)/src/Lib/Geometry/libGeometry.a \ + $(top_builddir)/src/Lib/Output/libOutput.a \ $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ $(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \ $(top_builddir)/src/Lib/e00/libe00.a \ diff --git a/src/Prep/GSHHS/Makefile.am b/src/Prep/GSHHS/Makefile.am index 12a525d0..db961331 100644 --- a/src/Prep/GSHHS/Makefile.am +++ b/src/Prep/GSHHS/Makefile.am @@ -8,6 +8,7 @@ gshhs_SOURCES = \ gshhs_LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ $(top_builddir)/src/Lib/Geometry/libGeometry.a \ + $(top_builddir)/src/Lib/Output/libOutput.a \ $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ -lsgdebug -lsgbucket -lsgmath -lsgmisc -lsgstructure -lsgxml \ -lgenpolyclip -lz diff --git a/src/Prep/ShapeFile/Makefile.am b/src/Prep/ShapeFile/Makefile.am index 1470f003..fc208f3f 100644 --- a/src/Prep/ShapeFile/Makefile.am +++ b/src/Prep/ShapeFile/Makefile.am @@ -7,6 +7,7 @@ noaa_decode_SOURCES = shape-decode.cxx LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ $(top_builddir)/src/Lib/Geometry/libGeometry.a \ + $(top_builddir)/src/Lib/Output/libOutput.a \ $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ $(top_builddir)/src/Lib/shapelib/libshape.a \ $(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \ diff --git a/src/Prep/ShapeFile/shape-decode.cxx b/src/Prep/ShapeFile/shape-decode.cxx index 72b68d62..3e287ee6 100644 --- a/src/Prep/ShapeFile/shape-decode.cxx +++ b/src/Prep/ShapeFile/shape-decode.cxx @@ -33,6 +33,7 @@ #include #include +#include #include #include @@ -52,11 +53,12 @@ SG_USING_STD( string ); SG_USING_STD( map ); map area_code_map; -int use_area_code_map=0; -int continue_on_errors=0; -int area_column=4,code_column=3; +bool use_area_code_map = false; +int area_column = 4, code_column = 3; double max_segment = 0.0; // zero = no splitting int width_column = -1; +bool ivmap = false; +double line_width_pad = 0.0; void load_noaa_area_codes() { area_code_map["Urban (1990 Enhanced)"]="Urban"; @@ -112,7 +114,7 @@ void load_noaa_area_codes() { area_code_map["Glaciers"]="Glacier"; area_code_map["No Data"]="Null"; - use_area_code_map=1; + use_area_code_map = true; } // return the type of the shapefile record @@ -220,19 +222,52 @@ AreaType get_shapefile_type(DBFHandle& hDBF, int rec) { SG_LOG( SG_GENERAL, SG_INFO, " raw area = " << area ); SG_LOG( SG_GENERAL, SG_INFO, " code = " << code ); - if (use_area_code_map) { - area=area_code_map[area]; + if ( use_area_code_map ) { + area = area_code_map[area]; SG_LOG( SG_GENERAL, SG_INFO, "cooked area = " << area ); } + if ( ivmap ) { + if ( area == "40" ) { + area = "WhiteLine"; + } else if ( area == "156" ) { + area = "YellowLine"; + } else { + SG_LOG( SG_GENERAL, SG_ALERT, + "Unknown ivmap area type = " << area ); + } + } + return get_area_type( area ); } // get attribute value of 'width' from shapefile -int get_shapefile_width(DBFHandle& hDBF, int rec) { +float get_shapefile_width(DBFHandle& hDBF, int rec) { string wstring = DBFReadStringAttribute( hDBF, rec, width_column ); - SG_LOG( SG_GENERAL, SG_DEBUG, "wstring = " << wstring ); - return atoi(wstring.c_str()); + SG_LOG( SG_GENERAL, SG_DEBUG, "width string = " << wstring ); + double width = atof(wstring.c_str()); + return width + line_width_pad; +} + +// return an arbitrary shape file attribute as a string +string get_attribute(DBFHandle& hDBF, int rec, int column) { + + string attrib = DBFReadStringAttribute( hDBF, rec, column ); + + // strip leading spaces + while ( attrib[0] == ' ' ) { + attrib = attrib.substr(1, attrib.length() - 1); + } + // strip trailing spaces + while ( attrib[attrib.length() - 1] == ' ' ) { + attrib = attrib.substr(0, attrib.length() - 1); + } + // strip other junk encountered + while ( (int)attrib[attrib.length() - 1] == 9 ) { + attrib = attrib.substr(0, attrib.length() - 1); + } + + return attrib; } void processPolygon(SHPObject* psShape, @@ -298,23 +333,172 @@ void processPolygon(SHPObject* psShape, tgChopNormalPolygon(work_dir, area, shape, preserve3D); } + +// Calculate theta of angle (a, b, c) +double calc_angle(Point3D a, Point3D b, Point3D c) { + Point3D u, v; + double udist, vdist, uv_dot, tmp; + + // u . v = ||u|| * ||v|| * cos(theta) + + u.setx( b.x() - a.x() ); + u.sety( b.y() - a.y() ); + udist = sqrt( u.x() * u.x() + u.y() * u.y() ); + // printf("udist = %.6f\n", udist); + + v.setx( b.x() - c.x() ); + v.sety( b.y() - c.y() ); + vdist = sqrt( v.x() * v.x() + v.y() * v.y() ); + // printf("vdist = %.6f\n", vdist); + + uv_dot = u.x() * v.x() + u.y() * v.y(); + // printf("uv_dot = %.6f\n", uv_dot); + + tmp = uv_dot / (udist * vdist); + // printf("tmp = %.6f\n", tmp); + + return acos(tmp); +} + + +// given an input line, try to identify locations in the line where +// the data doubles back on itself and remove those (IVLab data glitch +// work around) +tg::Line fixDegenerateLine1( tg::Line line ) { + tg::Line result; + + int size = line.getPointCount(); + if ( size < 3 ) { + return line; + } + + // prime the pump, and copy over the starting point + Point3D a = line.getPoint(0); + Point3D b = line.getPoint(1); + result.addPoint( a ); + + for ( int i = 1; i < size - 1; i++ ) { + Point3D c = line.getPoint(i+1); + double angle = calc_angle( a, b, c ); + if ( angle < 0 ) { + cout << "OOPS, negative angle!!!! PLEASE CHECK ME!!!" << endl; + } + if ( angle > SGD_PI_2 ) { + // good angle, march forward (small angles indicate the + // line is doubling back on itself) + result.addPoint(b); + a = b; + b = c; + } else { + // degenerate angle, skip point c and try next point + printf("IVLab data glitch at = %.8f %.8f\n", b.lon(), b.lat()); + } + } + + // copy the final point. + result.addPoint( line.getPoint(size-1) ); + + return result; +} + + +// given an input line, try to identify locations in the line where +// the points are packed really tightly and remove those (IVLab data glitch +// work around) +tg::Line fixDegenerateLine2( tg::Line line ) { + const double min_segment_length_m = 10.0; + tg::Line result; + + int size = line.getPointCount(); + if ( size <= 2 ) { + return line; + } + + // prime the pump, and copy over the starting point + Point3D a = line.getPoint(0); + result.addPoint( a ); + + for ( int i = 1; i < size - 1; i++ ) { + Point3D b = line.getPoint(i); + double az1, az2, dist; + geo_inverse_wgs_84( 0.0, a.y(), a.x(), b.y(), b.x(), + &az1, &az2, &dist ); + if ( dist > min_segment_length_m ) { + // distance ok, march forward (small distances could + // indicate clusters of points that might cause problems. + result.addPoint(b); + a = b; + } else { + // distance too small, skip point b and try next point + printf("IVLab data glitch at = %.8f %.8f\n", b.lon(), b.lat()); + } + } + + // add the last point (we always want the first and last points of + // the line to minimize visual glitches + result.addPoint( line.getPoint(size-1) ); + + return result; +} + + +// given an input line, avoid a situation where there is a completely +// loop, i.e. after a long route, the road loops around and connects +// up with itself. Back off from the end of the line until there is a +// small gap of more than 1 meter. +tg::Line fixDegenerateLine3( tg::Line line ) { + tg::Line result; + + int size = line.getPointCount(); + if ( size <= 2 ) { + return line; + } + + // prime the pump + Point3D a = line.getPoint(0); + int end = size - 1; + + for ( int i = size - 1; i > 1; i-- ) { + end = i; + Point3D b = line.getPoint(end); + double az1, az2, dist; + geo_inverse_wgs_84( 0.0, a.y(), a.x(), b.y(), b.x(), + &az1, &az2, &dist ); + if ( dist > 2.0 ) { + // distance ok, go ahead and exit + break; + } else { + // distance too small, march backwards one step and try + // the previous point + printf("IVLab data loop at = %.8f %.8f\n", b.lon(), b.lat()); + } + } + + for ( int i = 0; i <= end; i++ ) { + result.addPoint( line.getPoint(i) ); + } + + return result; +} + + void processLine(SHPObject* psShape, const string& work_dir, AreaType area, - int linewidth) { + float linewidth) { int iPart,j=0,partEnd=psShape->nVertices; - + double minx = 200, miny = 200, maxx = -200, maxy = -200; for (iPart=psShape->nParts-1;iPart>=0;iPart--) { tg::Line line; TGPolygon shape; -#if 0 + // #if 0 SG_LOG( SG_GENERAL, SG_DEBUG, - "Part " << iPart + "In processLine() Part " << iPart << " type = " << SHPPartTypeName(psShape->panPartType[iPart]) << " start = " << psShape->panPartStart[iPart] << " end = " << partEnd << " # of points = " << partEnd-psShape->panPartStart[iPart]); -#endif + // #endif if (partEnd-psShape->panPartStart[iPart]<=1) { SG_LOG( SG_GENERAL, SG_WARN, @@ -323,20 +507,37 @@ void processLine(SHPObject* psShape, } shape.erase(); for (j=psShape->panPartStart[iPart];jpadfX[j] << ", " << psShape->padfY[j] << ")"); -#endif + if ( psShape->padfX[j] < minx ) { minx = psShape->padfX[j]; } + if ( psShape->padfX[j] > maxx ) { maxx = psShape->padfX[j]; } + if ( psShape->padfY[j] < miny ) { miny = psShape->padfY[j]; } + if ( psShape->padfY[j] > maxy ) { maxy = psShape->padfY[j]; } + //#endif line.addPoint(Point3D(psShape->padfX[j],psShape->padfY[j],0)); } + SG_LOG( SG_GENERAL, SG_DEBUG, "(" << minx << "," << miny << ") (" << maxx << "," << maxy << ")" ); + + if ( ivmap ) { + cout << "original line size = " << line.getPointCount() << endl; + line = fixDegenerateLine1( line ); + cout << "angle line size = " << line.getPointCount() << endl; + line = fixDegenerateLine2( line ); + cout << "dist line size = " << line.getPointCount() << endl; + line = fixDegenerateLine3( line ); + cout << "cycle line size = " << line.getPointCount() << endl; + } + partEnd=psShape->panPartStart[iPart]; tg::makePolygon(line,linewidth,shape); if ( max_segment > 1.0 ) { shape = tgPolygonSplitLongEdges( shape, max_segment ); } + cout << "hole flag = " << shape.get_hole_flag(0) << endl; tgChopNormalPolygon(work_dir, area, shape, false); } } @@ -406,14 +607,16 @@ void usage(char* progname) { } int main( int argc, char **argv ) { - int i, j; - int pointwidth = 500, linewidth = 50, force_linewidth = -1; + int i; + int pointwidth = 500; + float force_linewidth = -1, linewidth = 50.0; int start_record = 0; char* progname=argv[0]; SGPath programPath(progname); string force_area_type = ""; + bool continue_on_errors = false; - sglog().setLogLevels( SG_ALL, SG_INFO ); + sglog().setLogLevels( SG_ALL, SG_DEBUG ); if (programPath.file()=="noaa-decode") { SG_LOG( SG_GENERAL, SG_INFO, "Entering noaa-decode mode" ); @@ -426,7 +629,7 @@ int main( int argc, char **argv ) { if (!strcmp(argv[1],"--line-width")) { if (argc<3) usage(progname); - force_linewidth=atoi(argv[2]); + force_linewidth=atof(argv[2]); argv+=2; argc-=2; } else if (!strcmp(argv[1],"--point-width")) { @@ -456,7 +659,7 @@ int main( int argc, char **argv ) { } else if (!strcmp(argv[1],"--continue-on-errors")) { argv++; argc--; - continue_on_errors=1; + continue_on_errors=true; } else if (!strcmp(argv[1],"--line-width-column")) { if (argc<3) usage(progname); @@ -469,6 +672,20 @@ int main( int argc, char **argv ) { start_record=atoi(argv[2]); argv+=2; argc-=2; + } else if (!strcmp(argv[1],"--ivlanemarkings")) { + argv++; + argc--; + ivmap = true; + width_column = 6; + } else if (!strcmp(argv[1],"--ivlane")) { + argv++; + argc--; + ivmap = true; + width_column = 3; + } else if (!strcmp(argv[1],"--pad-width")) { + line_width_pad = atof(argv[2]); + argv+=2; + argc-=2; } else if (!strcmp(argv[1],"--help")) { usage(progname); } else @@ -533,7 +750,8 @@ int main( int argc, char **argv ) { nShapeType != SHPT_ARCM && nShapeType != SHPT_POLYGON && nShapeType != SHPT_POLYGONM && - nShapeType != SHPT_POLYGONZ) { + nShapeType != SHPT_POLYGONZ) + { SG_LOG( SG_GENERAL, SG_ALERT, "Can only handle 2D-Points, 2D-Multipoints, " "2D-Arcs and Polygons (2D and 3D)" ); @@ -591,11 +809,18 @@ int main( int argc, char **argv ) { area = get_area_type( force_area_type ); } - if ( force_linewidth == -1) { - if (width_column != -1) // line width from shape file - linewidth = get_shapefile_width(hDBF, i); - } else - linewidth = force_linewidth; + if ( force_linewidth < 0 ) { + if (width_column != -1) // line width from shape file + linewidth = get_shapefile_width(hDBF, i); + } else + linewidth = force_linewidth; + + if ( ivmap ) { + // get skip type + string skip = get_attribute(hDBF, i, 7); + + // FIXME: implement skip logic + } SG_LOG( SG_GENERAL, SG_DEBUG, " record type = " << SHPTypeName(psShape->nSHPType) ); diff --git a/src/Prep/TGVPF/Makefile.am b/src/Prep/TGVPF/Makefile.am index 490d68e7..720e1a48 100644 --- a/src/Prep/TGVPF/Makefile.am +++ b/src/Prep/TGVPF/Makefile.am @@ -5,6 +5,7 @@ tgvpf_SOURCES = tgvpf.cxx tgvpf_LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ $(top_builddir)/src/Lib/Geometry/libGeometry.a \ + $(top_builddir)/src/Lib/Output/libOutput.a \ $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ $(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \ $(top_builddir)/src/Lib/vpf/libvpf.a \ diff --git a/src/Prep/UserDef/Makefile.am b/src/Prep/UserDef/Makefile.am index 68deac25..ca0da508 100644 --- a/src/Prep/UserDef/Makefile.am +++ b/src/Prep/UserDef/Makefile.am @@ -5,6 +5,7 @@ tguserdef_SOURCES = tguserdef.cxx tguserdef_LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ $(top_builddir)/src/Lib/Geometry/libGeometry.a \ + $(top_builddir)/src/Lib/Output/libOutput.a \ $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ $(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \ -lsgbucket -lsgmisc -lsgmath -lsgprops -lsgio -lsgstructure -lsgxml \