Ralf Gerlich:

This patch for TerraGear extends the shape-decode and noaa-decode tools to also handle shapefiles containing point- and line-data.

src/Prep/DemChop/Makefile.am

@@ -23,7 +23,7 @@
 #---------------------------------------------------------------------------
 
 
-bin_PROGRAMS = demchop hgtchop testassem
+bin_PROGRAMS = demchop hgtchop fillvoids testassem
 
 demchop_SOURCES = \
 	demchop.cxx point2d.hxx
@@ -41,6 +41,14 @@ hgtchop_LDADD = \
 	-lsgbucket -lsgmisc -lsgdebug -lsgxml -lz
 	$(base_LIBS)
 
+fillvoids_SOURCES = \
+	fillvoids.cxx
+
+fillvoids_LDADD = \
+	$(top_builddir)/src/Lib/Array/libArray.a \
+	-lsgbucket -lsgmisc -lsgdebug -lsgxml -lz
+	$(base_LIBS)
+
 testassem_SOURCES = \
 	testassem.cxx
 

src/Prep/ShapeFile/Makefile.am

@@ -2,18 +2,9 @@ bin_PROGRAMS = shape-decode noaa-decode
 
 shape_decode_SOURCES = shape-decode.cxx
 
-noaa_decode_SOURCES = noaa-decode.cxx
+noaa_decode_SOURCES = shape-decode.cxx
 
-shape_decode_LDADD = \
+LDADD = \
-	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
-	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
-	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \
-	$(top_builddir)/src/Lib/shapelib/libshape.a \
-	$(top_builddir)/src/Lib/TriangleJRS/libTriangleJRS.a \
-	-lsgdebug -lsgbucket -lsgmath -lsgmisc -lsgstructure -lsgxml \
-	-lgenpolyclip -lz
-
-noaa_decode_LDADD = \
 	$(top_builddir)/src/Lib/Polygon/libPolygon.a \
 	$(top_builddir)/src/Lib/Geometry/libGeometry.a \
 	$(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \

src/Prep/ShapeFile/shape-decode.cxx

@@ -1,10 +1,12 @@
 // main.cxx -- process shapefiles and extract polygon outlines,
-//             clipping against and sorting them into the revelant
+//             lines and points, clipping against and sorting
-//             tiles.
+//             them into the revelant tiles.
 //
 // Written by Curtis Olson, started February 1999.
+// Modified by Ralf Gerlich, August 2005.
 //
 // Copyright (C) 1999  Curtis L. Olson  - http://www.flightgear.org/~curt
+// Copyright (C) 2005  Ralf Gerlich     - http://home.easylink.de/rgerlich
 //
 // 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
@@ -26,9 +28,13 @@
 #include <simgear/compiler.h>
 
 #include STL_STRING
+#include <map>
 
 #include <simgear/debug/logstream.hxx>
+#include <simgear/misc/sg_path.hxx>
 
+#include <Geometry/line.hxx>
+#include <Geometry/util.hxx>
 #include <Polygon/chop.hxx>
 #include <Polygon/index.hxx>
 #include <Polygon/names.hxx>
@@ -41,6 +47,69 @@
 
 SG_USING_STD( cout );
 SG_USING_STD( string );
+SG_USING_STD( map );
+
+map<string,string> area_code_map;
+int use_area_code_map=0;
+int continue_on_errors=0;
+int area_column=4,code_column=3;
+
+void load_noaa_area_codes() {
+    area_code_map["Urban (1990 Enhanced)"]="Urban";
+    area_code_map["Residential"]="Urban";
+    area_code_map["Commercial and Services"]="Urban";
+    area_code_map["Industrial"]="Urban";
+    area_code_map["Transportation, Communications and Utilities"]="Urban";
+    area_code_map["Transportation, Communications, and Utilities"]="Urban";
+    area_code_map["Industrial and Commercial Complexes"]="Urban";
+    area_code_map["Mixed Urban or Built-up Land"]="BuiltUpCover";
+    area_code_map["Mixed Urban or Built up Land"]="BuiltUpCover";
+    area_code_map["Other Urban or Built-up Land"]="BuiltUpCover";
+    area_code_map["Other Urban or Built up Land"]="BuiltUpCover";
+    
+    area_code_map["Cropland and Pasture"]="MixedCropPastureCover";
+    area_code_map["Orchards, Groves, Vineyards, Nurseries, and Ornamental Horticultural Areas"]="IrrCropPastureCover";
+    area_code_map["Orchards, Groves, Vineyards, Nurseries and Ornament"]="IrrCropPastureCover";
+    area_code_map["Orchards, Groves, Vineyards, Nurseries and Ornamental Hort"]="IrrCropPastureCover";
+    area_code_map["Confined Feeding Operations"]="MixedCropPastureCover";
+    area_code_map["Other Agricultural Land"]="IrrCropPastureCover";
+    
+    area_code_map["Herbaceous Rangeland"]="GrassCover";
+    area_code_map["Shrub and Brush Rangeland"]="ShrubCover";
+    area_code_map["Mixed Rangeland"]="ShrubGrassCover";
+    
+    area_code_map["Deciduous Forest Land"]="DeciduousBroadCover";
+    area_code_map["Evergreen Forest Land"]="EvergreenNeedleCover";
+    area_code_map["Mixed Forest Land"]="MixedForestCover";
+    
+    area_code_map["Streams and Canals"]="Stream";
+    area_code_map["Lakes"]="Lake";
+    area_code_map["Reservoirs"]="Reservoir";
+    area_code_map["Bays and Estuaries"]="Ocean";
+    
+    area_code_map["Forested Wetland"]="WoodedWetlandCover";
+    area_code_map["Nonforested Wetland"]="HerbWetlandCover";
+    
+    area_code_map["Dry Salt Flats"]="DryLake";
+    area_code_map["Beaches"]="DryLake";
+    area_code_map["Sandy Areas Other than Beaches"]="DryLake";
+    area_code_map["Bare Exposed Rock"]="BarrenCover";
+    area_code_map["Strip Mines, Quarries, and Gravel Pits"]="BarrenCover";
+    area_code_map["Transitional Areas"]="BarrenCover";
+    area_code_map["Mixed Barren Land"]="BarrenCover";
+    
+    area_code_map["Shrub and Brush Tundra"]="WoodedTundraCover";
+    area_code_map["Herbaceous Tundra"]="HerbTundraCover";
+    area_code_map["Bare Ground"]="BareTundraCover";
+    area_code_map["Wet Tundra"]="HerbTundraCover";
+    area_code_map["Mixed Tundra"]="MixedTundraCover";
+    
+    area_code_map["Perennial Snowfields"]="Glacier";
+    area_code_map["Glaciers"]="Glacier";
+    
+    area_code_map["No Data"]="Null";
+    use_area_code_map=1;
+}
 
 // return the type of the shapefile record
 AreaType get_shapefile_type(DBFHandle& hDBF, int rec) {
@@ -127,9 +196,9 @@ AreaType get_shapefile_type(DBFHandle& hDBF, int rec) {
 
 #endif
  
-    string area = DBFReadStringAttribute( hDBF, rec, 4 );
+    string area = DBFReadStringAttribute( hDBF, rec, area_column );
-    string test = DBFReadStringAttribute( hDBF, rec, 3 );
+    string code = DBFReadStringAttribute( hDBF, rec, code_column );
-    cout << "next record = " << test << endl;
+    cout << "next record = " << code << endl;
 
     // strip leading spaces
     while ( area[0] == ' ' ) {
@@ -144,24 +213,199 @@ AreaType get_shapefile_type(DBFHandle& hDBF, int rec) {
 	area = area.substr(0, area.length() - 1);
     }
 
-    SG_LOG( SG_GENERAL, SG_INFO, "  raw area = " << area );
+    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];
+	SG_LOG( SG_GENERAL, SG_INFO, "cooked area = " << area );
+    }
 
     return get_area_type( area );
 }
 
+void processPolygon(SHPObject* psShape,
+		    const string& work_dir,
+		    AreaType area,
+		    bool preserve3D) {
+    int iPart,j;
+    TGPolygon shape;
+#if 0
+    const char *pszPlus;
+    const char *pszPartType = SHPPartTypeName( psShape->panPartType[0] );
+#endif
+
+    shape.erase();
+    for ( j = 0, iPart = 1; j < psShape->nVertices; j++ ) {
+
+	if( iPart < psShape->nParts
+	    && psShape->panPartStart[iPart] == j )
+	{
+	    iPart++;
+#if 0
+	    pszPartType = SHPPartTypeName( psShape->panPartType[iPart] );
+	    pszPlus = "+";
+	} else {
+	    pszPlus = " ";
+#endif
+	}
+
+	shape.add_node( iPart - 1, 
+			Point3D(psShape->padfX[j],
+				psShape->padfY[j],
+				psShape->padfZ[j])
+			);
+#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
+    }
+    
+    // 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 );
+	}
+    }
+
+    tgChopNormalPolygon(work_dir, area, shape, preserve3D);
+}
+
+void processLine(SHPObject* psShape,
+		 const string& work_dir,
+		 AreaType area,
+		 int linewidth) {
+    int iPart,j=0,partEnd=psShape->nVertices;
+
+    for (iPart=psShape->nParts-1;iPart>=0;iPart--) {
+	tg::Line line;
+	TGPolygon shape;
+#if 0
+	SG_LOG( SG_GENERAL, SG_DEBUG,
+		   "Part " << iPart
+		<< " type = " << SHPPartTypeName(psShape->panPartType[iPart])
+		<< " start = " << psShape->panPartStart[iPart]
+		<< " end = " << partEnd
+		<< " # of points = " << partEnd-psShape->panPartStart[iPart]);
+#endif
+
+	if (partEnd-psShape->panPartStart[iPart]<=1) {
+	    SG_LOG( SG_GENERAL, SG_WARN,
+	    	    "Skipping line with less than two points" );
+	    continue;
+	}
+	shape.erase();
+	for (j=psShape->panPartStart[iPart];j<partEnd;j++) {
+#if 0
+	    SG_LOG( SG_GENERAL, SG_DEBUG,
+	    	       "   Point ("
+		    << psShape->padfX[j] << ", "
+		    << psShape->padfY[j] << ")");
+#endif
+	    line.addPoint(Point3D(psShape->padfX[j],psShape->padfY[j],0));
+	}
+	partEnd=psShape->panPartStart[iPart];
+	tg::makePolygon(line,linewidth,shape);
+	tgChopNormalPolygon(work_dir, area, shape, false);
+    }
+}
+
+void processPoints(SHPObject* psShape,
+		   const string& work_dir,
+		   AreaType area,
+		   int pointwidth) {
+    TGPolygon shape;
+    int j;
+
+    for (j=0;j<psShape->nVertices;j++) {
+#if 0
+	SG_LOG( SG_GENERAL, SG_DEBUG,
+		   "   Point ("
+		<< psShape->padfX[j] << ", "
+		<< psShape->padfY[j] << ")");
+#endif
+	shape.erase();
+	tg::makePolygon(Point3D(psShape->padfX[j],psShape->padfY[j],0),
+			pointwidth,
+			shape);
+	tgChopNormalPolygon(work_dir, area, shape, false);
+    }
+}
+
+void usage(char* progname) {
+    SG_LOG( SG_GENERAL, SG_ALERT, "Usage: " << progname 
+	    << " [--line-width width] [--point-width width]"
+	       " [--area-column col] [--code-col col]"
+	       " [--continue-on-errors]"
+	       " <shape_file> <work_dir> [ area_string ]" );
+    exit(-1);
+}
 
 int main( int argc, char **argv ) {
-    TGPolygon shape;
     int i, j;
+    int pointwidth=500,linewidth=50;
+    char* progname=argv[0];
+    SGPath programPath(progname);
+    string force_area_type = "";
 
     sglog().setLogLevels( SG_ALL, SG_DEBUG );
 
-    if ( argc < 3 ) {
+    if (programPath.file()=="noaa-decode") {
-	SG_LOG( SG_GENERAL, SG_ALERT, "Usage: " << argv[0] 
+	SG_LOG( SG_GENERAL, SG_INFO, "Entering noaa-decode mode" );
-		<< " <shape_file> <work_dir> [ area_string ]" );
+	area_column=2;
-	exit(-1);
+	code_column=1;
+	load_noaa_area_codes();
     }
 
+    while (argc>1) {
+	if (!strcmp(argv[1],"--line-width")) {
+	    if (argc<3)
+		usage(progname);
+	    linewidth=atoi(argv[2]);
+	    argv+=2;
+	    argc-=2;
+	} else if (!strcmp(argv[1],"--point-width")) {
+	    if (argc<3)
+		usage(progname);
+	    pointwidth=atoi(argv[2]);
+	    argv+=2;
+	    argc-=2;
+	} else if (!strcmp(argv[1],"--area-column")) {
+	    if (argc<3)
+		usage(progname);
+	    area_column=atoi(argv[2]);
+	    argv+=2;
+	    argc-=2;
+	} else if (!strcmp(argv[1],"--code-column")) {
+	    if (argc<3)
+		usage(progname);
+	    code_column=atoi(argv[2]);
+	    argv+=2;
+	    argc-=2;
+	} else if (!strcmp(argv[1],"--continue-on-errors")) {
+	    continue_on_errors=1;
+	} else if (!strcmp(argv[1],"--help")) {
+	    usage(progname);
+	} else
+		break;
+    }
+    
+    if (argc<3)
+	usage(progname);
+    
     SG_LOG( SG_GENERAL, SG_DEBUG, "Opening " << argv[1] << " for reading." );
 
     // make work directory
@@ -177,7 +421,6 @@ int main( int argc, char **argv ) {
     // 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];
     }
@@ -210,27 +453,60 @@ int main( int argc, char **argv ) {
 
     SG_LOG( SG_GENERAL, SG_INFO, "shape file records = " << nEntities << endl );
 
-    string shapetype = SHPTypeName( nShapeType );
+    if ( nShapeType != SHPT_POINT &&
-
+         nShapeType != SHPT_POINTM &&
-    if ( shapetype != string("Polygon") ) {
+         nShapeType != SHPT_MULTIPOINT &&
-	SG_LOG( SG_GENERAL, SG_ALERT, "Can't handle non-polygon shape files" );
+         nShapeType != SHPT_MULTIPOINTM &&
+         nShapeType != SHPT_ARC &&
+         nShapeType != SHPT_ARCM &&
+         nShapeType != SHPT_POLYGON &&
+         nShapeType != SHPT_POLYGONM &&
+         nShapeType != SHPT_POLYGONZ) {
+	SG_LOG( SG_GENERAL, SG_ALERT,
+		"Can only handle 2D-Points, 2D-Multipoints, "
+		"2D-Arcs and Polygons (2D and 3D)" );
 	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 );
 
+	if (psShape->nSHPType==0) {
+	    SG_LOG( SG_GENERAL, SG_DEBUG, "Skipping NULL record " << i);
+	    continue;
+	}
+	
 	SG_LOG( SG_GENERAL, SG_DEBUG, "Processing record = " << i 
 		<< "  rings = " << psShape->nParts
 		<< "  total vertices = " << psShape->nVertices );
 
+	if (psShape->nParts<0) {
+	    SG_LOG( SG_GENERAL, SG_ALERT,
+	    	    "Record with negative part count, "
+		    "file may be corrupted!");
+	    if (continue_on_errors)
+		continue;
+	    else
+		exit(-1);
+	}
+	
+	if (psShape->nParts==0 &&
+	    psShape->nSHPType!=SHPT_POINT &&
+	    psShape->nSHPType!=SHPT_POINTM &&
+	    psShape->nSHPType!=SHPT_MULTIPOINT &&
+	    psShape->nSHPType!=SHPT_MULTIPOINTM) {
+	    SG_LOG( SG_GENERAL, SG_ALERT,
+	    	    "Non-point record with zero part count, "
+		    "file may be corrupted!");
+	    if (continue_on_errors)
+		continue;
+	    else
+		exit(-1);
+	}
+	
 	AreaType area = DefaultArea;
 	if ( force_area_type.length() == 0 ) {
 	    area = get_shapefile_type(hDBF, i);
@@ -260,81 +536,79 @@ int main( int argc, char **argv ) {
                 psShape->dfZMax, psShape->dfMMax );
 #endif
 
-        for ( j = 0, iPart = 1; j < psShape->nVertices; j++ ) {
+	if ( area == OceanArea ) {
-            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 );
-
-	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 );
-	    tgChopNormalPolygon(work_dir, area, shape, false);
-	} else if ( area == OceanArea ) {
 	    // interior of polygon is ocean, holes are islands
 
 	    SG_LOG(  SG_GENERAL, SG_ALERT, "Ocean area ... SKIPPING!" );
 
 	    // Ocean data now comes from GSHHS so we want to ignore
 	    // all other ocean data
-	    // tgChopPolygon(work_dir, area, shape, false);
+	    continue;
 	} else if ( area == VoidArea ) {
 	    // interior is ????
 
 	    // skip for now
 	    SG_LOG(  SG_GENERAL, SG_ALERT, "Void area ... SKIPPING!" );
 
-	    if ( shape.contours() > 1 ) {
+	    if ( psShape->nParts > 1 ) {
 		SG_LOG(  SG_GENERAL, SG_ALERT, "  Void area with holes!" );
 		// exit(-1);
 	    }
-
+	    
-	    // tgChopPolygon(work_dir, area, shape, false);
+	    continue;
 	} else if ( area == NullArea ) {
 	    // interior is ????
 
 	    // skip for now
 	    SG_LOG(  SG_GENERAL, SG_ALERT, "Null area ... SKIPPING!" );
 
-	    if ( shape.contours() > 1 ) {
+	    if ( psShape->nParts > 1 ) {
 		SG_LOG(  SG_GENERAL, SG_ALERT, "  Null area with holes!" );
 		// exit(-1);
 	    }
 
-	    // tgChopPolygon(work_dir, area, shape, false);
+	    continue;
-	} else {
-	    tgChopNormalPolygon(work_dir, area, shape, false);
 	}
+
+	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 );
+	}
+	
+	switch (psShape->nSHPType) {
+	case SHPT_POLYGON:
+	case SHPT_POLYGONM:
+		/* 2D Polygons */
+		processPolygon(psShape,work_dir,area,false);
+		break;
+	case SHPT_POLYGONZ:
+		/* 3D Polygons */
+		processPolygon(psShape,work_dir,area,true);
+		break;
+	case SHPT_ARC:
+	case SHPT_ARCM:
+		/* 2D Lines */
+		processLine(psShape,work_dir,area,linewidth);
+		break;
+	case SHPT_POINT:
+	case SHPT_POINTM:
+	case SHPT_MULTIPOINT:
+	case SHPT_MULTIPOINTM:
+		/* 2D Points */
+		processPoints(psShape,work_dir,area,pointwidth);
+		break;
+	default:
+		SG_LOG( SG_GENERAL, SG_WARN,
+			"Can't handle "
+			<< SHPTypeName(psShape->nSHPType)
+			<< " yet - skipping" );
+		break;
+	}
+	
+        SHPDestroyObject( psShape );
+
     }
 
     DBFClose( hDBF );