From ed22267dd9d5dc04ea5073e0e33b087a0a9737c0 Mon Sep 17 00:00:00 2001 From: curt Date: Tue, 8 Nov 2005 16:27:36 +0000 Subject: [PATCH] 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 | 10 +- src/Prep/ShapeFile/Makefile.am | 13 +- src/Prep/ShapeFile/shape-decode.cxx | 412 +++++++++++++++++++++++----- 3 files changed, 354 insertions(+), 81 deletions(-) diff --git a/src/Prep/DemChop/Makefile.am b/src/Prep/DemChop/Makefile.am index a7b1c307..a66dd35f 100644 --- a/src/Prep/DemChop/Makefile.am +++ b/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 diff --git a/src/Prep/ShapeFile/Makefile.am b/src/Prep/ShapeFile/Makefile.am index 34ad1f6b..1470f003 100644 --- a/src/Prep/ShapeFile/Makefile.am +++ b/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 = \ - $(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 = \ +LDADD = \ $(top_builddir)/src/Lib/Polygon/libPolygon.a \ $(top_builddir)/src/Lib/Geometry/libGeometry.a \ $(top_builddir)/src/Lib/poly2tri/libpoly2tri.a \ diff --git a/src/Prep/ShapeFile/shape-decode.cxx b/src/Prep/ShapeFile/shape-decode.cxx index 5296f09a..7b7921f9 100644 --- a/src/Prep/ShapeFile/shape-decode.cxx +++ b/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 -// tiles. +// lines and points, clipping against and sorting +// 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 #include STL_STRING +#include #include +#include +#include +#include #include #include #include @@ -41,6 +47,69 @@ SG_USING_STD( cout ); 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; + +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 test = DBFReadStringAttribute( hDBF, rec, 3 ); - cout << "next record = " << test << endl; + string area = DBFReadStringAttribute( hDBF, rec, area_column ); + string code = DBFReadStringAttribute( hDBF, rec, code_column ); + 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];jpadfX[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;jnVertices;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]" + " [ 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 ) { - SG_LOG( SG_GENERAL, SG_ALERT, "Usage: " << argv[0] - << " [ area_string ]" ); - exit(-1); + if (programPath.file()=="noaa-decode") { + SG_LOG( SG_GENERAL, SG_INFO, "Entering noaa-decode mode" ); + area_column=2; + 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 ( shapetype != string("Polygon") ) { - SG_LOG( SG_GENERAL, SG_ALERT, "Can't handle non-polygon shape files" ); + if ( nShapeType != SHPT_POINT && + nShapeType != SHPT_POINTM && + nShapeType != SHPT_MULTIPOINT && + 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++ ) { - 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 ) { + 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); - } else { - tgChopNormalPolygon(work_dir, area, shape, false); + continue; } + + 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 );