1
0
Fork 0

fix ogr-decode double save issue

- added some experimental (and unused) iterative tests
This commit is contained in:
Peter Sadrozinski 2013-03-11 09:12:42 -04:00
parent 08f5c17221
commit a74edcd91b
5 changed files with 118 additions and 113 deletions

View file

@ -9,7 +9,7 @@
#include "tg_chopper.hxx"
void tgChopper::Clip( const tgPolygon& subject,
tgPolygon tgChopper::Clip( const tgPolygon& subject,
const std::string& type,
SGBucket& b )
{
@ -42,24 +42,12 @@ void tgChopper::Clip( const tgPolygon& subject,
base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), max.getLatitudeDeg()) );
base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), max.getLatitudeDeg()) );
SG_LOG(SG_GENERAL, SG_DEBUG, "shape contours = " << subject.Contours() );
for ( unsigned int ii = 0; ii < subject.Contours(); ii++ ) {
SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << subject.GetContour(ii).GetHole() );
}
result = tgPolygon::Intersect( subject, base );
SG_LOG(SG_GENERAL, SG_DEBUG, "result contours = " << result.Contours() );
for ( unsigned int ii = 0; ii < result.Contours(); ii++ ) {
SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << result.GetContour(ii).GetHole() );
}
if ( subject.GetPreserve3D() ) {
result.InheritElevations( subject );
}
if ( result.Contours() > 0 ) {
result.SetPreserve3D( subject.GetPreserve3D() );
if ( subject.GetPreserve3D() ) {
result.InheritElevations( subject );
result.SetPreserve3D( true );
}
result.SetTexParams( subject.GetTexParams() );
if ( subject.GetTexMethod() == TG_TEX_BY_GEODE ) {
// need to set center latitude for geodetic texturing
@ -71,6 +59,8 @@ void tgChopper::Clip( const tgPolygon& subject,
bp_map[b.gen_index()].push_back( result );
lock.unlock();
}
return result;
}
void tgChopper::Add( const tgPolygon& subject, const std::string& type )
@ -97,98 +87,96 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type )
if ( b_min == b_max ) {
// shape entirely contained in a single bucket, write and bail
Clip( subject, type, b_min );
return;
}
} else {
SGBucket b_cur;
int dx, dy;
SGBucket b_cur;
int dx, dy;
sgBucketDiff(b_min, b_max, &dx, &dy);
SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" );
SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy );
sgBucketDiff(b_min, b_max, &dx, &dy);
SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" );
SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy );
if ( (dx > 2880) || (dy > 1440) )
throw sg_exception("something is really wrong in split_polygon()!!!!");
if ( (dx > 2880) || (dy > 1440) )
throw sg_exception("something is really wrong in split_polygon()!!!!");
if ( dy <= 1 ) {
// we are down to at most two rows, write each column and then bail
double min_center_lat = b_min.get_center_lat();
double min_center_lon = b_min.get_center_lon();
for ( int j = 0; j <= dy; ++j ) {
for ( int i = 0; i <= dx; ++i ) {
b_cur = sgBucketOffset(min_center_lon, min_center_lat, i, j);
Clip( subject, type, b_cur );
}
}
return;
}
if ( dy <= 1 ) {
// we are down to at most two rows, write each column and then bail
double min_center_lat = b_min.get_center_lat();
double min_center_lon = b_min.get_center_lon();
for ( int j = 0; j <= dy; ++j ) {
for ( int i = 0; i <= dx; ++i ) {
b_cur = sgBucketOffset(min_center_lon, min_center_lat, i, j);
Clip( subject, type, b_cur );
// we have two or more rows left, split in half (along a
// horizontal dividing line) and recurse with each half
// find mid point (integer math)
int mid = (dy + 1) / 2 - 1;
// determine horizontal clip line
SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, mid);
double clip_line = b_clip.get_center_lat();
if ( (clip_line >= -90.0 + SG_HALF_BUCKET_SPAN)
&& (clip_line < 90.0 - SG_HALF_BUCKET_SPAN) )
clip_line += SG_HALF_BUCKET_SPAN;
else if ( clip_line < -89.0 )
clip_line = -89.0;
else if ( clip_line >= 89.0 )
clip_line = 90.0;
else {
SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << clip_line );
}
{
//
// Crop bottom area (hopefully by putting this in it's own
// scope we can shorten the life of some really large data
// structures to reduce memory use)
//
SG_LOG( SG_GENERAL, SG_DEBUG, "Generating bottom half (" << bb.getMin().getLatitudeDeg() << "-" << clip_line << ")" );
tgPolygon bottom, bottom_clip;
bottom.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMin().getLatitudeDeg()) );
bottom.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMin().getLatitudeDeg()) );
bottom.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) );
bottom.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) );
bottom_clip = tgPolygon::Intersect( subject, bottom );
if ( (bottom_clip.TotalNodes() > 0) && (bottom_clip.TotalNodes() != subject.TotalNodes()) ) {
Add( bottom_clip, type );
}
}
return;
}
// we have two or more rows left, split in half (along a
// horizontal dividing line) and recurse with each half
{
//
// Crop top area (hopefully by putting this in it's own scope
// we can shorten the life of some really large data
// structures to reduce memory use)
//
// find mid point (integer math)
int mid = (dy + 1) / 2 - 1;
SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half (" << clip_line << "-" << bb.getMax().getLatitudeDeg() << ")" );
// determine horizontal clip line
SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, mid);
double clip_line = b_clip.get_center_lat();
if ( (clip_line >= -90.0 + SG_HALF_BUCKET_SPAN)
&& (clip_line < 90.0 - SG_HALF_BUCKET_SPAN) )
clip_line += SG_HALF_BUCKET_SPAN;
else if ( clip_line < -89.0 )
clip_line = -89.0;
else if ( clip_line >= 89.0 )
clip_line = 90.0;
else {
SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << clip_line );
}
tgPolygon top, top_clip;
{
//
// Crop bottom area (hopefully by putting this in it's own
// scope we can shorten the life of some really large data
// structures to reduce memory use)
//
top.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) );
top.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) );
top.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMax().getLatitudeDeg()) );
top.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMax().getLatitudeDeg()) );
SG_LOG( SG_GENERAL, SG_DEBUG, "Generating bottom half (" << bb.getMin().getLatitudeDeg() << "-" << clip_line << ")" );
top_clip = tgPolygon::Intersect( subject, top );
tgPolygon bottom, bottom_clip;
bottom.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMin().getLatitudeDeg()) );
bottom.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMin().getLatitudeDeg()) );
bottom.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) );
bottom.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) );
bottom_clip = tgPolygon::Intersect( subject, bottom );
if ( (bottom_clip.TotalNodes() > 0) && (bottom_clip.TotalNodes() != subject.TotalNodes()) ) {
Add( bottom_clip, type );
if ( (top_clip.TotalNodes() > 0) && (top_clip.TotalNodes() != subject.TotalNodes()) ) {
Add( top_clip, type );
}
}
}
{
//
// Crop top area (hopefully by putting this in it's own scope
// we can shorten the life of some really large data
// structures to reduce memory use)
//
SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half (" << clip_line << "-" << bb.getMax().getLatitudeDeg() << ")" );
tgPolygon top, top_clip;
top.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) );
top.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) );
top.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMax().getLatitudeDeg()) );
top.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMax().getLatitudeDeg()) );
top_clip = tgPolygon::Intersect( subject, top );
if ( (top_clip.TotalNodes() > 0) && (top_clip.TotalNodes() != subject.TotalNodes()) ) {
Add( top_clip, type );
} else
return;
}
}
long int tgChopper::GenerateIndex( std::string path )
@ -265,4 +253,4 @@ void tgChopper::Save( void )
gzclose( fp );
}
}
}

View file

@ -18,7 +18,7 @@ public:
private:
long int GenerateIndex( std::string path );
void Clip( const tgPolygon& subject, const std::string& type, SGBucket& b );
tgPolygon Clip( const tgPolygon& subject, const std::string& type, SGBucket& b );
void Chop( const tgPolygon& subject, const std::string& type );
std::string root_path;

View file

@ -36,6 +36,23 @@ void tgRectangle::setMax (const SGGeod &p)
_max = p;
}
void tgRectangle::expandBy(const tgRectangle& r)
{
if ( r.getMin().getLongitudeDeg() < _min.getLongitudeDeg() ) {
_min.setLongitudeDeg( r.getMin().getLongitudeDeg() );
}
if ( r.getMin().getLatitudeDeg() < _min.getLatitudeDeg() ) {
_min.setLatitudeDeg( r.getMin().getLatitudeDeg() );
}
if ( r.getMax().getLongitudeDeg() > _max.getLongitudeDeg() ) {
_max.setLongitudeDeg( r.getMax().getLongitudeDeg() );
}
if ( r.getMax().getLatitudeDeg() > _max.getLatitudeDeg() ) {
_max.setLatitudeDeg( r.getMax().getLatitudeDeg() );
}
}
void tgRectangle::sanify ()
{
double tmp;

View file

@ -90,6 +90,13 @@ public:
*/
virtual void setMax (const SGGeod& p);
/**
* Expand the rectangle to encompass a given rectangle
*
* @param r The rectangle to expand by
*/
virtual void expandBy(const tgRectangle& r);
/**
* Make the rectangle sane.
*

View file

@ -212,6 +212,8 @@ void Decoder::run()
// as long as we have geometry to parse, do so
while (!global_workQueue.empty()) {
OGRFeature *poFeature = global_workQueue.pop();
// SG_LOG( SG_GENERAL, SG_INFO, " remaining features is " << global_workQueue.size() );
if ( poFeature ) {
OGRGeometry *poGeometry = poFeature->GetGeometryRef();
@ -315,9 +317,9 @@ void Decoder::run()
}
}
OGRMultiLineString* multils=(OGRMultiLineString*)poGeometry;
for (int i=0;i<multils->getNumGeometries();i++) {
processLineString((OGRLineString*)poGeometry, area_type_name, width, texture_lines);
OGRMultiLineString* multilines=(OGRMultiLineString*)poGeometry;
for (int i=0;i<multilines->getNumGeometries();i++) {
processLineString((OGRLineString*)(multilines->getGeometryRef(i)), area_type_name, width, texture_lines);
}
break;
}
@ -388,9 +390,7 @@ void processLayer(OGRLayer* poLayer, tgChopper& results )
/* 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 );
@ -441,6 +441,7 @@ void processLayer(OGRLayer* poLayer, tgChopper& results )
}
// Now process the workqueue with threads
// this just generates all the tgPolygons
std::vector<Decoder *> decoders;
for (int i=0; i<num_threads; i++) {
Decoder* decoder = new Decoder( poCT, area_type_field, point_width_field, line_width_field, results );
@ -448,19 +449,11 @@ void processLayer(OGRLayer* poLayer, tgChopper& results )
decoders.push_back( decoder );
}
#if 0
while (!global_workQueue.empty()) {
sleep(1);
}
#endif
// Then wait until they are finished
for (unsigned int i=0; i<decoders.size(); i++) {
decoders[i]->join();
}
results.Save();
OCTDestroyCoordinateTransformation ( poCT );
}
@ -659,7 +652,6 @@ int main( int argc, char **argv ) {
SG_LOG( SG_GENERAL, SG_ALERT, "Processing datasource " << datasource );
OGRLayer *poLayer;
if (argc>3) {
for (int i=3;i<argc;i++) {
poLayer = poDS->GetLayerByName( argv[i] );
@ -681,9 +673,10 @@ int main( int argc, char **argv ) {
}
}
results.Save();
OGRDataSource::DestroyDataSource( poDS );
SG_LOG(SG_GENERAL, SG_ALERT, "Saving to buckets");
results.Save();
return 0;
}