fly/terragear
1
0
Fork 0
Code Activity

Added ogrdecode, which decodes all OGR-supported vector input formats to TerraGear polygons. Might replace shape-decode one day.

Browse source
This commit is contained in:
Ralf Gerlich 2007-11-09 13:16:47 +01:00
parent 2a4560f312
commit 01f5b80457
4 changed files with 566 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
configure.ac
View file

@ -138,11 +138,15 @@ else
fi
GDAL_LIBS="$GDAL_LIBS $GDAL_DEP_LIBS"
if test "$USE_OGR" = "1"; then
OGR_DECODE=OGRDecode
fi
AC_SUBST(GDAL_LIBS)
AC_SUBST(GDAL_CFLAGS)
AC_SUBST(USE_GDAL)
AC_SUBST(USE_OGR)
AC_SUBST(OGR_DECODE)
fi
dnl Check for MS Windows environment
@ -460,6 +464,7 @@ AC_CONFIG_FILES([ \
src/Prep/GSHHS/Makefile \
src/Prep/MergerClipper/Makefile \
src/Prep/Photo/Makefile \
src/Prep/OGRDecode/Makefile \
src/Prep/ShapeFile/Makefile \
src/Prep/TGVPF/Makefile \
src/Prep/Terra/Makefile \

5
src/Prep/Makefile.am
View file

@ -12,4 +12,7 @@ SUBDIRS = \
Terra \
TerraFit \
Tower \
UserDef
UserDef \
$(OGR_DECODE)
EXTRA_SUBDIRS=OGRDecode

16
src/Prep/OGRDecode/Makefile.am Normal file
View file

@ -0,0 +1,16 @@
bin_PROGRAMS = ogr-decode
ogr_decode_SOURCES = ogr-decode.cxx
ogr_decode_CXXFLAGS=$(GDAL_CFLAGS)
ogr_decode_LDADD = \
$(GDAL_LIBS) \
$(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/TriangleJRS/libTriangleJRS.a \
-lsgdebug -lsgbucket -lsgmath -lsgmisc -lsgstructure -lsgxml \
-lgenpolyclip -lz
INCLUDES = -I$(top_srcdir)/src/Lib

541
src/Prep/OGRDecode/ogr-decode.cxx Normal file
View file

@ -0,0 +1,541 @@
// ogr-decode.cxx -- process OGR layers and extract polygon outlines,
// lines and points, clipping against and sorting
// them into the relevant tiles.
//
// Written by Ralf Gerlich, started February 2007.
//
// Copyright (C) 2007 Ralf Gerlich - ralf.gerlich@custom-scenery.org
// Loosely based on shape-decode.cxx by Curtis L. Olson
//
// 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 <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>
#include <Polygon/polygon.hxx>
#include <ogrsf_frmts.h>
#ifdef _MSC_VER
# include <Win32/mkdir.hpp>
#endif
SG_USING_STD( cout );
SG_USING_STD( string );
SG_USING_STD( map );
int line_width=50;
string line_width_col;
int point_width=500;
string point_width_col;
string area_type="Default";
string area_type_col;
int continue_on_errors=0;
int max_segment_length=0; // ==0 => don't split
int start_record=0;
bool use_attribute_query=false;
string attribute_query;
bool use_spatial_query=false;
double spat_min_x, spat_min_y, spat_max_x, spat_max_y;
void processPoint(OGRPoint* poGeometry,
const string& work_dir,
const string& area_type,
int width)
{
TGPolygon shape;
shape.erase();
tg::makePolygon(Point3D(poGeometry->getX(),poGeometry->getY(),0),
width,
shape);
if ( max_segment_length > 0 ) {
shape = tgPolygonSplitLongEdges( shape, max_segment_length );
}
tgChopNormalPolygon(work_dir, area_type, shape, false);
}
void processLineString(OGRLineString* poGeometry,
const string& work_dir,
const string& area_type,
int width)
{
tg::Line line;
TGPolygon shape;
if (poGeometry->getNumPoints()<2) {
SG_LOG( SG_GENERAL, SG_WARN,
"Skipping line with less than two points" );
return;
}
shape.erase();
for (int i=0;i<poGeometry->getNumPoints();i++) {
line.addPoint(Point3D(poGeometry->getX(i),poGeometry->getY(i),0));
}
tg::makePolygon(line,width,shape);
if ( max_segment_length > 0 ) {
shape = tgPolygonSplitLongEdges( shape, max_segment_length );
}
tgChopNormalPolygon(work_dir, area_type, shape, false);
}
void processPolygon(OGRPolygon* poGeometry,
const string& work_dir,
const string& area_type)
{
bool preserve3D = ((poGeometry->getGeometryType()&wkb25DBit)==wkb25DBit);
TGPolygon shape;
shape.erase();
OGRLinearRing *ring;
// first add the outer ring
ring=poGeometry->getExteriorRing();
for (int i=0;i<ring->getNumPoints();i++) {
shape.add_node(0,
Point3D(ring->getX(i),
ring->getY(i),
ring->getZ(i)));
}
shape.set_hole_flag( 0, 0 );
// then add the inner rings
for ( int j = 0 ; j < poGeometry->getNumInteriorRings() ; j ++ ) {
ring=poGeometry->getInteriorRing( j );
for (int i=0;i<ring->getNumPoints();i++) {
shape.add_node(j+1,
Point3D(ring->getX(i),
ring->getY(i),
ring->getZ(i)));
}
shape.set_hole_flag( j+1 , 1 );
}
if ( max_segment_length > 0 ) {
shape = tgPolygonSplitLongEdges( shape, max_segment_length );
}
tgChopNormalPolygon(work_dir, area_type, shape, preserve3D);
}
void processLayer(OGRLayer* poLayer,
const string& work_dir) {
int feature_count=poLayer->GetFeatureCount();
if (feature_count!=-1 && start_record>0 && start_record>=feature_count) {
SG_LOG( SG_GENERAL, SG_ALERT, "Layer has only " << feature_count << " records, but start record is set to " << start_record );
exit( 1 );
}
/* determine the indices of the required columns */
OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
string layername=poFDefn->GetName();
int point_width_field=-1, line_width_field=-1, area_type_field=-1;
if (!point_width_col.empty()) {
point_width_field=poFDefn->GetFieldIndex(point_width_col.c_str());
if (point_width_field==-1) {
SG_LOG( SG_GENERAL, SG_ALERT, "Field " << point_width_col << " for point-width not found in layer" );
if (!continue_on_errors)
exit( 1 );
}
}
if (!line_width_col.empty()) {
line_width_field=poFDefn->GetFieldIndex(line_width_col.c_str());
if (line_width_field==-1) {
SG_LOG( SG_GENERAL, SG_ALERT, "Field " << line_width_col << " for line-width not found in layer" );
if (!continue_on_errors)
exit( 1 );
}
}
if (!area_type_col.empty()) {
area_type_field=poFDefn->GetFieldIndex(area_type_col.c_str());
if (area_type_field==-1) {
SG_LOG( SG_GENERAL, SG_ALERT, "Field " << area_type_col << " for line-width not found in layer" );
if (!continue_on_errors)
exit( 1 );
}
}
/* setup a transformation to WGS84 */
OGRSpatialReference *oSourceSRS, oTargetSRS;
oSourceSRS=poLayer->GetSpatialRef();
if (oSourceSRS == NULL) {
SG_LOG( SG_GENERAL, SG_ALERT, "Layer " << layername << " has no defined spatial reference system" );
exit( 1 );
}
char* srsWkt;
oSourceSRS->exportToWkt(&srsWkt);
SG_LOG( SG_GENERAL, SG_DEBUG, "Source spatial reference system: " << srsWkt );
OGRFree(srsWkt);
oTargetSRS.SetWellKnownGeogCS( "WGS84" );
OGRCoordinateTransformation *poCT;
poCT = OGRCreateCoordinateTransformation(oSourceSRS, &oTargetSRS);
/* setup attribute and spatial queries */
if (use_spatial_query) {
double trans_min_x,trans_min_y,trans_max_x,trans_max_y;
/* do a simple reprojection of the source SRS */
OGRCoordinateTransformation *poCTinverse;
poCTinverse = OGRCreateCoordinateTransformation(&oTargetSRS, oSourceSRS);
trans_min_x=spat_min_x;
trans_min_y=spat_min_y;
trans_max_x=spat_max_x;
trans_max_y=spat_max_y;
poCTinverse->Transform(1,&trans_min_x,&trans_min_y);
poCTinverse->Transform(1,&trans_max_x,&trans_max_y);
poLayer->SetSpatialFilterRect(trans_min_x, trans_min_y,
trans_max_x, trans_max_y);
}
if (use_attribute_query) {
if (poLayer->SetAttributeFilter(attribute_query.c_str()) != OGRERR_NONE) {
SG_LOG( SG_GENERAL, SG_ALERT, "Error in query expression '" << attribute_query << "'" );
exit( 1 );
}
}
OGRFeature *poFeature;
poLayer->SetNextByIndex(start_record);
for ( ; (poFeature = poLayer->GetNextFeature()) != NULL; OGRFeature::DestroyFeature( poFeature ) )
{
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if (poGeometry==NULL) {
SG_LOG( SG_GENERAL, SG_INFO, "Found feature without geometry!" );
if (!continue_on_errors) {
SG_LOG( SG_GENERAL, SG_ALERT, "Aborting!" );
exit( 1 );
} else {
continue;
}
}
assert(poGeometry!=NULL);
OGRwkbGeometryType geoType=wkbFlatten(poGeometry->getGeometryType());
if (geoType!=wkbPoint && geoType!=wkbMultiPoint &&
geoType!=wkbLineString && geoType!=wkbMultiLineString &&
geoType!=wkbPolygon && geoType!=wkbMultiPolygon) {
SG_LOG( SG_GENERAL, SG_INFO, "Unknown feature" );
continue;
}
string area_type_name=area_type;
if (area_type_field!=-1) {
area_type_name=poFeature->GetFieldAsString(area_type_field);
}
if ( is_ocean_area(area_type_name) ) {
// 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
continue;
} else if ( is_void_area(area_type_name) ) {
// interior is ????
// skip for now
SG_LOG( SG_GENERAL, SG_ALERT, "Void area ... SKIPPING!" );
continue;
} else if ( is_null_area(area_type_name) ) {
// interior is ????
// skip for now
SG_LOG( SG_GENERAL, SG_ALERT, "Null area ... SKIPPING!" );
continue;
}
poGeometry->transform( poCT );
switch (geoType) {
case wkbPoint: {
SG_LOG( SG_GENERAL, SG_DEBUG, "Point feature" );
int width=point_width;
if (point_width_field!=-1) {
width=poFeature->GetFieldAsInteger(point_width_field);
}
processPoint((OGRPoint*)poGeometry,work_dir,area_type_name,width);
break;
}
case wkbMultiPoint: {
SG_LOG( SG_GENERAL, SG_DEBUG, "MultiPoint feature" );
int width=point_width;
if (point_width_field!=-1) {
width=poFeature->GetFieldAsInteger(point_width_field);
}
OGRMultiPoint* multipt=(OGRMultiPoint*)poGeometry;
for (int i=0;i<multipt->getNumGeometries();i++) {
processPoint((OGRPoint*)(multipt->getGeometryRef(i)),work_dir,area_type_name,width);
}
break;
}
case wkbLineString: {
SG_LOG( SG_GENERAL, SG_DEBUG, "LineString feature" );
int width=line_width;
if (line_width_field!=-1) {
width=poFeature->GetFieldAsInteger(line_width_field);
}
processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width);
break;
}
case wkbMultiLineString: {
SG_LOG( SG_GENERAL, SG_DEBUG, "MultiLineString feature" );
int width=line_width;
if (line_width_field!=-1) {
width=poFeature->GetFieldAsInteger(line_width_field);
}
OGRMultiLineString* multils=(OGRMultiLineString*)poGeometry;
for (int i=0;i<multils->getNumGeometries();i++) {
processLineString((OGRLineString*)(multils->getGeometryRef(i)),work_dir,area_type_name,width);
}
break;
}
case wkbPolygon: {
SG_LOG( SG_GENERAL, SG_DEBUG, "Polygon feature" );
processPolygon((OGRPolygon*)poGeometry,work_dir,area_type_name);
break;
}
case wkbMultiPolygon: {
SG_LOG( SG_GENERAL, SG_DEBUG, "MultiPolygon feature" );
OGRMultiPolygon* multipoly=(OGRMultiPolygon*)poGeometry;
for (int i=0;i<multipoly->getNumGeometries();i++) {
processPolygon((OGRPolygon*)(multipoly->getGeometryRef(i)),work_dir,area_type_name);
}
break;
}
}
}
OCTDestroyCoordinateTransformation ( poCT );
}
void usage(char* progname) {
SG_LOG( SG_GENERAL, SG_ALERT, "Usage: " <<
progname << " [options...] <work_dir> <datasource> [<layername>...]" );
SG_LOG( SG_GENERAL, SG_ALERT, "Options:" );
SG_LOG( SG_GENERAL, SG_ALERT, "--line-width width" );
SG_LOG( SG_GENERAL, SG_ALERT, " Width in meters for the lines" );
SG_LOG( SG_GENERAL, SG_ALERT, "--line-width-column colname" );
SG_LOG( SG_GENERAL, SG_ALERT, " Use value from colname as width for lines" );
SG_LOG( SG_GENERAL, SG_ALERT, " Overrides --line-width if present" );
SG_LOG( SG_GENERAL, SG_ALERT, "--point-width width" );
SG_LOG( SG_GENERAL, SG_ALERT, " Size in meters of the squares generated from points" );
SG_LOG( SG_GENERAL, SG_ALERT, "--point-width-column colname" );
SG_LOG( SG_GENERAL, SG_ALERT, " Use value from colname as width for points" );
SG_LOG( SG_GENERAL, SG_ALERT, " Overrides --point-width if present" );
SG_LOG( SG_GENERAL, SG_ALERT, "--area-type type" );
SG_LOG( SG_GENERAL, SG_ALERT, " Area type for all objects from file" );
SG_LOG( SG_GENERAL, SG_ALERT, "--area-type-column colname" );
SG_LOG( SG_GENERAL, SG_ALERT, " Use string from colname as area type" );
SG_LOG( SG_GENERAL, SG_ALERT, " Overrides --area-type if present" );
SG_LOG( SG_GENERAL, SG_ALERT, "--continue-on-errors" );
SG_LOG( SG_GENERAL, SG_ALERT, " Continue even if the file seems fishy" );
SG_LOG( SG_GENERAL, SG_ALERT, "--max-segment max_segment_length" );
SG_LOG( SG_GENERAL, SG_ALERT, " Maximum segment length in meters" );
SG_LOG( SG_GENERAL, SG_ALERT, "--start-record record_no" );
SG_LOG( SG_GENERAL, SG_ALERT, " Start processing at the specified record number (first record num=0)" );
SG_LOG( SG_GENERAL, SG_ALERT, "--where attrib_query" );
SG_LOG( SG_GENERAL, SG_ALERT, " Use an attribute query (like SQL WHERE)" );
SG_LOG( SG_GENERAL, SG_ALERT, "--spat xmin ymin xmax ymax" );
SG_LOG( SG_GENERAL, SG_ALERT, " spatial query extents" );
SG_LOG( SG_GENERAL, SG_ALERT, "" );
SG_LOG( SG_GENERAL, SG_ALERT, "<work_dir>" );
SG_LOG( SG_GENERAL, SG_ALERT, " Directory to put the polygon files in" );
SG_LOG( SG_GENERAL, SG_ALERT, "<datasource>" );
SG_LOG( SG_GENERAL, SG_ALERT, " The datasource from which to fetch the data" );
SG_LOG( SG_GENERAL, SG_ALERT, "<layername>..." );
SG_LOG( SG_GENERAL, SG_ALERT, " The layers to process." );
SG_LOG( SG_GENERAL, SG_ALERT, " If no layer is given, all layers in the datasource are used" );
exit(-1);
}
int main( int argc, char **argv ) {
char* progname=argv[0];
string datasource,work_dir;
sglog().setLogLevels( SG_ALL, SG_INFO );
while (argc>1) {
if (!strcmp(argv[1],"--line-width")) {
if (argc<3)
usage(progname);
line_width=atoi(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--line-width-column")) {
if (argc<3)
usage(progname);
line_width_col=argv[2];
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--point-width")) {
if (argc<3)
usage(progname);
point_width=atoi(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--point-width-column")) {
if (argc<3)
usage(progname);
point_width_col=argv[2];
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--area-type")) {
if (argc<3)
usage(progname);
area_type=argv[2];
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--area-type-column")) {
if (argc<3)
usage(progname);
area_type_col=argv[2];
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--continue-on-errors")) {
argv++;
argc--;
continue_on_errors=1;
} else if (!strcmp(argv[1],"--max-segment")) {
if (argc<3)
usage(progname);
max_segment_length=atoi(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--start-record")) {
if (argc<3)
usage(progname);
start_record=atoi(argv[2]);
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--where")) {
if (argc<3)
usage(progname);
use_attribute_query=true;
attribute_query=argv[2];
argv+=2;
argc-=2;
} else if (!strcmp(argv[1],"--spat")) {
if (argc<6)
usage(progname);
use_spatial_query=true;
spat_min_x=atof(argv[2]);
spat_min_y=atof(argv[3]);
spat_max_x=atof(argv[4]);
spat_max_y=atof(argv[5]);
argv+=5;
argc-=5;
} else if (!strcmp(argv[1],"--help")) {
usage(progname);
} else
break;
}
if (argc<3)
usage(progname);
work_dir=argv[1];
datasource=argv[2];
#ifdef _MSC_VER
fg_mkdir( work_dir.c_str() );
#else
string command = "mkdir -p " + work_dir;
system( command.c_str() );
#endif
// initialize persistant polygon counter
string counter_file = work_dir + "/../poly_counter";
poly_index_init( counter_file );
SG_LOG( SG_GENERAL, SG_DEBUG, "Opening datasource " << datasource << " for reading." );
OGRRegisterAll();
OGRDataSource *poDS;
poDS = OGRSFDriverRegistrar::Open( datasource.c_str(), FALSE );
if( poDS == NULL )
{
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening datasource " << datasource );
exit( 1 );
}
OGRLayer *poLayer;
if (argc>3) {
for (int i=3;i<argc;i++) {
poLayer = poDS->GetLayerByName( argv[i] );
if (poLayer == NULL )
{
SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource );
exit( 1 );
}
processLayer(poLayer, work_dir);
}
} else {
for (int i=0;i<poDS->GetLayerCount();i++) {
poLayer = poDS->GetLayer(i);
assert(poLayer != NULL);
processLayer(poLayer, work_dir);
}
}
OGRDataSource::DestroyDataSource( poDS );
return 0;
}