- fix ogr-decode missing landclass
sometimes, sgBucketSpan returned an extra row ( a sliver ) chop would pick this and save it as the landclass output - instead of the full bucket poly. construct wouldn't have any poly for the area, and mark it as ocean
This commit is contained in:
Peter Sadrozinski
parent
d42b7269cf
commit
a6098f3c4a
7 changed files with 273 additions and 181 deletions
|
|
@ -1,8 +1,8 @@
|
|||
/*******************************************************************************
|
||||
* *
|
||||
* Author : Angus Johnson *
|
||||
* Version : 5.1.2 *
|
||||
* Date : 25 February 2013 *
|
||||
* Version : 5.1.4 *
|
||||
* Date : 24 March 2013 *
|
||||
* Website : http://www.angusj.com *
|
||||
* Copyright : Angus Johnson 2010-2013 *
|
||||
* *
|
||||
|
|
@ -543,14 +543,18 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2,
|
|||
IntPoint &ip, bool UseFullInt64Range)
|
||||
{
|
||||
double b1, b2;
|
||||
if (SlopesEqual(edge1, edge2, UseFullInt64Range)) return false;
|
||||
if (SlopesEqual(edge1, edge2, UseFullInt64Range))
|
||||
{
|
||||
if (edge2.ybot > edge1.ybot) ip.Y = edge2.ybot;
|
||||
else ip.Y = edge1.ybot;
|
||||
return false;
|
||||
}
|
||||
else if (NEAR_ZERO(edge1.dx))
|
||||
{
|
||||
ip.X = edge1.xbot;
|
||||
if (NEAR_EQUAL(edge2.dx, HORIZONTAL))
|
||||
{
|
||||
ip.Y = edge2.ybot;
|
||||
} else
|
||||
else
|
||||
{
|
||||
b2 = edge2.ybot - (edge2.xbot / edge2.dx);
|
||||
ip.Y = Round(ip.X / edge2.dx + b2);
|
||||
|
|
@ -560,14 +564,14 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2,
|
|||
{
|
||||
ip.X = edge2.xbot;
|
||||
if (NEAR_EQUAL(edge1.dx, HORIZONTAL))
|
||||
{
|
||||
ip.Y = edge1.ybot;
|
||||
} else
|
||||
else
|
||||
{
|
||||
b1 = edge1.ybot - (edge1.xbot / edge1.dx);
|
||||
ip.Y = Round(ip.X / edge1.dx + b1);
|
||||
}
|
||||
} else
|
||||
}
|
||||
else
|
||||
{
|
||||
b1 = edge1.xbot - edge1.ybot * edge1.dx;
|
||||
b2 = edge2.xbot - edge2.ybot * edge2.dx;
|
||||
|
|
@ -586,7 +590,8 @@ bool IntersectPoint(TEdge &edge1, TEdge &edge2,
|
|||
ip.X = edge1.xtop;
|
||||
ip.Y = edge1.ytop;
|
||||
return TopX(edge2, edge1.ytop) < edge1.xtop;
|
||||
} else
|
||||
}
|
||||
else
|
||||
{
|
||||
ip.X = edge2.xtop;
|
||||
ip.Y = edge2.ytop;
|
||||
|
|
@ -1761,17 +1766,11 @@ void Clipper::IntersectEdges(TEdge *e1, TEdge *e2,
|
|||
}
|
||||
else if ( e1Contributing )
|
||||
{
|
||||
if ((e2Wc == 0 || e2Wc == 1) &&
|
||||
(m_ClipType != ctIntersection ||
|
||||
e2->polyType == ptSubject || (e2->windCnt2 != 0)))
|
||||
DoEdge1(e1, e2, pt);
|
||||
if (e2Wc == 0 || e2Wc == 1) DoEdge1(e1, e2, pt);
|
||||
}
|
||||
else if ( e2contributing )
|
||||
{
|
||||
if ((e1Wc == 0 || e1Wc == 1) &&
|
||||
(m_ClipType != ctIntersection ||
|
||||
e1->polyType == ptSubject || (e1->windCnt2 != 0)))
|
||||
DoEdge2(e1, e2, pt);
|
||||
if (e1Wc == 0 || e1Wc == 1) DoEdge2(e1, e2, pt);
|
||||
}
|
||||
else if ( (e1Wc == 0 || e1Wc == 1) &&
|
||||
(e2Wc == 0 || e2Wc == 1) && !e1stops && !e2stops )
|
||||
|
|
@ -2204,28 +2203,28 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge)
|
|||
TEdge* e = GetNextInAEL( horzEdge , dir );
|
||||
while( e )
|
||||
{
|
||||
if ( e->xcurr == horzEdge->xtop && !eMaxPair )
|
||||
{
|
||||
if (SlopesEqual(*e, *horzEdge->nextInLML, m_UseFullRange))
|
||||
{
|
||||
//if output polygons share an edge, they'll need joining later ...
|
||||
if (horzEdge->outIdx >= 0 && e->outIdx >= 0)
|
||||
AddJoin(horzEdge->nextInLML, e, horzEdge->outIdx);
|
||||
break; //we've reached the end of the horizontal line
|
||||
}
|
||||
else if (e->dx < horzEdge->nextInLML->dx)
|
||||
//we really have got to the end of the intermediate horz edge so quit.
|
||||
//nb: More -ve slopes follow more +ve slopes ABOVE the horizontal.
|
||||
break;
|
||||
}
|
||||
|
||||
TEdge* eNext = GetNextInAEL( e, dir );
|
||||
|
||||
if (eMaxPair ||
|
||||
((dir == dLeftToRight) && (e->xcurr <= horzRight)) ||
|
||||
((dir == dRightToLeft) && (e->xcurr >= horzLeft)))
|
||||
((dir == dLeftToRight) && (e->xcurr < horzRight)) ||
|
||||
((dir == dRightToLeft) && (e->xcurr > horzLeft)))
|
||||
{
|
||||
//ok, so far it looks like we're still in range of the horizontal edge
|
||||
if ( e->xcurr == horzEdge->xtop && !eMaxPair )
|
||||
{
|
||||
if (SlopesEqual(*e, *horzEdge->nextInLML, m_UseFullRange))
|
||||
{
|
||||
//if output polygons share an edge, they'll need joining later ...
|
||||
if (horzEdge->outIdx >= 0 && e->outIdx >= 0)
|
||||
AddJoin(horzEdge->nextInLML, e, horzEdge->outIdx);
|
||||
break; //we've reached the end of the horizontal line
|
||||
}
|
||||
else if (e->dx < horzEdge->nextInLML->dx)
|
||||
//we really have got to the end of the intermediate horz edge so quit.
|
||||
//nb: More -ve slopes follow more +ve slopes ABOVE the horizontal.
|
||||
break;
|
||||
}
|
||||
|
||||
//so far we're still in range of the horizontal edge
|
||||
if( e == eMaxPair )
|
||||
{
|
||||
//horzEdge is evidently a maxima horizontal and we've arrived at its end.
|
||||
|
|
@ -2261,8 +2260,8 @@ void Clipper::ProcessHorizontal(TEdge *horzEdge)
|
|||
}
|
||||
SwapPositionsInAEL( horzEdge, e );
|
||||
}
|
||||
else if( (dir == dLeftToRight && e->xcurr > horzRight && m_SortedEdges) ||
|
||||
(dir == dRightToLeft && e->xcurr < horzLeft && m_SortedEdges) ) break;
|
||||
else if( (dir == dLeftToRight && e->xcurr >= horzRight) ||
|
||||
(dir == dRightToLeft && e->xcurr <= horzLeft) ) break;
|
||||
e = eNext;
|
||||
} //end while
|
||||
|
||||
|
|
@ -2345,7 +2344,7 @@ void Clipper::BuildIntersectList(const long64 botY, const long64 topY)
|
|||
{
|
||||
e->prevInSEL = e->prevInAEL;
|
||||
e->nextInSEL = e->nextInAEL;
|
||||
e->tmpX = TopX( *e, topY );
|
||||
e->xcurr = TopX( *e, topY );
|
||||
e = e->nextInAEL;
|
||||
}
|
||||
|
||||
|
|
@ -2359,9 +2358,10 @@ void Clipper::BuildIntersectList(const long64 botY, const long64 topY)
|
|||
{
|
||||
TEdge *eNext = e->nextInSEL;
|
||||
IntPoint pt;
|
||||
if(e->tmpX > eNext->tmpX &&
|
||||
IntersectPoint(*e, *eNext, pt, m_UseFullRange))
|
||||
if(e->xcurr > eNext->xcurr)
|
||||
{
|
||||
if (!IntersectPoint(*e, *eNext, pt, m_UseFullRange) && e->xcurr > eNext->xcurr +1)
|
||||
throw clipperException("Intersection error");
|
||||
if (pt.Y > botY)
|
||||
{
|
||||
pt.Y = botY;
|
||||
|
|
@ -2452,7 +2452,7 @@ void Clipper::DoMaxima(TEdge *e, long64 topY)
|
|||
if (!eNext) throw clipperException("DoMaxima error");
|
||||
IntersectEdges( e, eNext, IntPoint(X, topY), ipBoth );
|
||||
SwapPositionsInAEL(e, eNext);
|
||||
eNext = eNext->nextInAEL;
|
||||
eNext = e->nextInAEL;
|
||||
}
|
||||
if( e->outIdx < 0 && eMaxPair->outIdx < 0 )
|
||||
{
|
||||
|
|
@ -2600,25 +2600,22 @@ void Clipper::FixupOutPolygon(OutRec &outRec)
|
|||
|
||||
void Clipper::BuildResult(Polygons &polys)
|
||||
{
|
||||
int k = 0;
|
||||
polys.resize(m_PolyOuts.size());
|
||||
polys.reserve(m_PolyOuts.size());
|
||||
for (PolyOutList::size_type i = 0; i < m_PolyOuts.size(); ++i)
|
||||
{
|
||||
if (m_PolyOuts[i]->pts)
|
||||
{
|
||||
Polygon* pg = &polys[k];
|
||||
pg->clear();
|
||||
Polygon pg;
|
||||
OutPt* p = m_PolyOuts[i]->pts;
|
||||
do
|
||||
{
|
||||
pg->push_back(p->pt);
|
||||
pg.push_back(p->pt);
|
||||
p = p->prev;
|
||||
} while (p != m_PolyOuts[i]->pts);
|
||||
//make sure each polygon has at least 3 vertices ...
|
||||
if (pg->size() < 3) pg->clear(); else k++;
|
||||
if (pg.size() > 2)
|
||||
polys.push_back(pg);
|
||||
}
|
||||
}
|
||||
polys.resize(k);
|
||||
}
|
||||
//------------------------------------------------------------------------------
|
||||
|
||||
|
|
@ -3178,14 +3175,14 @@ PolyOffsetBuilder(const Polygons& in_polys, Polygons& out_polys,
|
|||
|
||||
switch (jointype)
|
||||
{
|
||||
case jtRound: //limit defaults to 0.125
|
||||
if (limit <= 0) limit = 0.125;
|
||||
case jtRound:
|
||||
if (limit <= 0) limit = 0.25;
|
||||
else if (limit > std::fabs(delta)) limit = std::fabs(delta);
|
||||
break;
|
||||
case jtMiter: //limit defaults to twice delta's width ...
|
||||
case jtMiter:
|
||||
if (limit < 2) limit = 2;
|
||||
break;
|
||||
default: //otherwise limit is unused
|
||||
default: //unused
|
||||
limit = 1;
|
||||
}
|
||||
m_RMin = 2.0/(limit*limit);
|
||||
|
|
@ -3275,9 +3272,8 @@ private:
|
|||
|
||||
void AddPoint(const IntPoint& pt)
|
||||
{
|
||||
Polygon::size_type len = m_curr_poly->size();
|
||||
if (len == m_curr_poly->capacity())
|
||||
m_curr_poly->reserve(len + buffLength);
|
||||
if (m_curr_poly->size() == m_curr_poly->capacity())
|
||||
m_curr_poly->reserve(m_curr_poly->capacity() + buffLength);
|
||||
m_curr_poly->push_back(pt);
|
||||
}
|
||||
//------------------------------------------------------------------------------
|
||||
|
|
@ -3410,38 +3406,62 @@ void SimplifyPolygons(Polygons &polys, PolyFillType fillType)
|
|||
}
|
||||
//------------------------------------------------------------------------------
|
||||
|
||||
void CleanPolygon(Polygon& in_poly, Polygon& out_poly, double distance)
|
||||
bool PointsAreClose(IntPoint pt1, IntPoint pt2, long64 distSqrd)
|
||||
{
|
||||
//delta = proximity in units/pixels below which vertices
|
||||
//will be stripped. Default ~= sqrt(2) so when adjacent
|
||||
//vertices have both x & y coords within 1 unit, then
|
||||
//the second vertex will be stripped.
|
||||
int len = in_poly.size();
|
||||
if (len < 3)
|
||||
out_poly.resize(0);
|
||||
else
|
||||
out_poly.resize(in_poly.size());
|
||||
|
||||
int d = (int)(distance * distance);
|
||||
IntPoint p = in_poly[0];
|
||||
int j = 1;
|
||||
for (int i = 1; i < len; i++)
|
||||
{
|
||||
if ((in_poly[i].X - p.X) * (in_poly[i].X - p.X) +
|
||||
(in_poly[i].Y - p.Y) * (in_poly[i].Y - p.Y) <= d)
|
||||
continue;
|
||||
out_poly[j] = in_poly[i];
|
||||
p = in_poly[i];
|
||||
j++;
|
||||
}
|
||||
p = in_poly[j - 1];
|
||||
if ((in_poly[0].X - p.X) * (in_poly[0].X - p.X) +
|
||||
(in_poly[0].Y - p.Y) * (in_poly[0].Y - p.Y) <= d)
|
||||
j--;
|
||||
if (j < len)
|
||||
out_poly.resize(j);
|
||||
long64 dx = pt1.X - pt2.X;
|
||||
long64 dy = pt1.Y - pt2.Y;
|
||||
return ((dx * dx) + (dy * dy) <= distSqrd);
|
||||
}
|
||||
//------------------------------------------------------------------------------
|
||||
|
||||
void CleanPolygon(Polygon& in_poly, Polygon& out_poly, double distance)
|
||||
{
|
||||
//distance = proximity in units/pixels below which vertices
|
||||
//will be stripped. Default ~= sqrt(2).
|
||||
int highI = in_poly.size() -1;
|
||||
long64 d = (int)(distance * distance);
|
||||
while (highI > 0 && PointsAreClose(in_poly[highI], in_poly[0], d)) highI--;
|
||||
if (highI < 2)
|
||||
{
|
||||
out_poly.clear();
|
||||
return;
|
||||
}
|
||||
out_poly.resize(highI + 1);
|
||||
bool UseFullRange = FullRangeNeeded(in_poly);
|
||||
IntPoint pt = in_poly[highI];
|
||||
int i = 0;
|
||||
int k = 0;
|
||||
for (;;)
|
||||
{
|
||||
if (i >= highI) break;
|
||||
int j = i + 1;
|
||||
|
||||
if (PointsAreClose(pt, in_poly[j], d))
|
||||
{
|
||||
i = j + 1;
|
||||
while (i <= highI && PointsAreClose(pt, in_poly[i], d)) i++;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (PointsAreClose(in_poly[i], in_poly[j], d) ||
|
||||
SlopesEqual(pt, in_poly[i], in_poly[j], UseFullRange))
|
||||
{
|
||||
i = j;
|
||||
continue;
|
||||
}
|
||||
|
||||
pt = in_poly[i++];
|
||||
out_poly[k++] = pt;
|
||||
}
|
||||
|
||||
if (i <= highI) out_poly[k++] = in_poly[i];
|
||||
if (k > 2 && SlopesEqual(out_poly[k -2], out_poly[k -1], out_poly[0], UseFullRange))
|
||||
k--;
|
||||
if (k < 3) out_poly.clear();
|
||||
else if (k <= highI) out_poly.resize(k);
|
||||
}
|
||||
//------------------------------------------------------------------------------
|
||||
|
||||
void CleanPolygons(Polygons& in_polys, Polygons& out_polys, double distance)
|
||||
{
|
||||
for (Polygons::size_type i = 0; i < in_polys.size(); ++i)
|
||||
|
|
|
|||
|
|
@ -1,8 +1,8 @@
|
|||
/*******************************************************************************
|
||||
* *
|
||||
* Author : Angus Johnson *
|
||||
* Version : 5.1.2 *
|
||||
* Date : 25 February 2013 *
|
||||
* Version : 5.1.4 *
|
||||
* Date : 24 March 2013 *
|
||||
* Website : http://www.angusj.com *
|
||||
* Copyright : Angus Johnson 2010-2013 *
|
||||
* *
|
||||
|
|
@ -134,7 +134,6 @@ struct TEdge {
|
|||
double dx;
|
||||
long64 deltaX;
|
||||
long64 deltaY;
|
||||
long64 tmpX;
|
||||
PolyType polyType;
|
||||
EdgeSide side;
|
||||
int windDelta; //1 or -1 depending on winding direction
|
||||
|
|
|
|||
|
|
@ -8,6 +8,9 @@
|
|||
|
||||
|
||||
#include "tg_chopper.hxx"
|
||||
#include "tg_shapefile.hxx"
|
||||
|
||||
unsigned int strip_clip = 0;
|
||||
|
||||
tgPolygon tgChopper::Clip( const tgPolygon& subject,
|
||||
const std::string& type,
|
||||
|
|
@ -42,7 +45,7 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
|
|||
base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), max.getLatitudeDeg()) );
|
||||
base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), max.getLatitudeDeg()) );
|
||||
|
||||
result = tgPolygon::Intersect( subject, base );
|
||||
result = tgPolygon::Intersect( base, subject );
|
||||
if ( result.Contours() > 0 ) {
|
||||
if ( subject.GetPreserve3D() ) {
|
||||
result.InheritElevations( subject );
|
||||
|
|
@ -63,17 +66,33 @@ tgPolygon tgChopper::Clip( const tgPolygon& subject,
|
|||
return result;
|
||||
}
|
||||
|
||||
// Pass in the center lat for clipping buckets from the row.
|
||||
// We can't rely on sgBucketOffset, as rounding error sometimes causes it to look like there are 2 rows
|
||||
// (the first being a sliver)
|
||||
// This leads to using that poly as the subject - which leads to having no usable polygon for this row.
|
||||
void tgChopper::ClipRow( const tgPolygon& subject, const double& center_lat, const std::string& type )
|
||||
{
|
||||
tgRectangle bb = subject.GetBoundingBox();
|
||||
SGBucket b_min( bb.getMin() );
|
||||
SGBucket b_max( bb.getMax() );
|
||||
double min_center_lon = b_min.get_center_lon();
|
||||
int dx, dy;
|
||||
|
||||
sgBucketDiff(b_min, b_max, &dx, &dy);
|
||||
|
||||
for ( int i = 0; i <= dx; ++i ) {
|
||||
SGBucket b_cur = sgBucketOffset(min_center_lon, center_lat, i, 0);
|
||||
Clip( subject, type, b_cur );
|
||||
}
|
||||
}
|
||||
|
||||
void tgChopper::Add( const tgPolygon& subject, const std::string& type )
|
||||
{
|
||||
tgRectangle bb;
|
||||
SGGeod p;
|
||||
|
||||
// bail out immediately if polygon is empty
|
||||
if ( subject.Contours() == 0 )
|
||||
return;
|
||||
|
||||
bb = subject.GetBoundingBox();
|
||||
|
||||
tgRectangle bb = subject.GetBoundingBox();
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " min = " << bb.getMin() << " max = " << bb.getMax() );
|
||||
|
||||
// find buckets for min, and max points of convex hull.
|
||||
|
|
@ -81,99 +100,58 @@ void tgChopper::Add( const tgPolygon& subject, const std::string& type )
|
|||
// polygons that span the date line
|
||||
SGBucket b_min( bb.getMin() );
|
||||
SGBucket b_max( bb.getMax() );
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket min = " << b_min );
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket max = " << b_max );
|
||||
SGBucket b_cur;
|
||||
int dx, dy;
|
||||
|
||||
if ( b_min == b_max ) {
|
||||
// shape entirely contained in a single bucket, write and bail
|
||||
Clip( subject, type, b_min );
|
||||
} else {
|
||||
SGBucket b_cur;
|
||||
int dx, dy;
|
||||
sgBucketDiff(b_min, b_max, &dx, &dy);
|
||||
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 == 0 )
|
||||
{
|
||||
// We just have a single row - no need to intersect first
|
||||
ClipRow( subject, b_min.get_center_lat(), type );
|
||||
}
|
||||
else
|
||||
{
|
||||
// Multiple rows - perform row intersection to reduce the number of bucket clips we need
|
||||
// since many shapes are narraw in some places, wide in others - bb will be at the widest part
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, "subject spans tile rows: bb is from lat " << bb.getMin().getLatitudeDeg() << " to " << bb.getMax().getLatitudeDeg() << " dy is " << dy );
|
||||
|
||||
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 );
|
||||
for ( int row = 0; row <= dy; row++ )
|
||||
{
|
||||
// Generate a clip rectangle - add some buffer on top and bottom, so we don't clip directly on an edge when we
|
||||
// clip the individual buckets
|
||||
// TODO : May no longer be necessary
|
||||
SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, row );
|
||||
double clip_bottom = b_clip.get_center_lat() - SG_HALF_BUCKET_SPAN; // + 0.01);
|
||||
double clip_top = b_clip.get_center_lat() + SG_HALF_BUCKET_SPAN; // + 0.01);
|
||||
tgPolygon clip_row, clipped;
|
||||
|
||||
SG_LOG( SG_GENERAL, SG_DEBUG, " row " << row << " center lat is " << b_clip.get_center_lat() << " clip_botton is " << clip_bottom << " clip_top is " << clip_top );
|
||||
|
||||
clip_row.AddNode( 0, SGGeod::fromDeg(-180.0, clip_bottom) );
|
||||
clip_row.AddNode( 0, SGGeod::fromDeg( 180.0, clip_bottom) );
|
||||
clip_row.AddNode( 0, SGGeod::fromDeg( 180.0, clip_top) );
|
||||
clip_row.AddNode( 0, SGGeod::fromDeg(-180.0, clip_top) );
|
||||
|
||||
clipped = tgPolygon::Intersect( clip_row, subject );
|
||||
if ( clipped.TotalNodes() > 0 ) {
|
||||
ClipRow( clipped, b_clip.get_center_lat(), type );
|
||||
|
||||
#if 0
|
||||
{
|
||||
char layer[32];
|
||||
char ds_name[64];
|
||||
sprintf(layer, "clipped_%d", strip_clip++ );
|
||||
sprintf(ds_name, "./stripped_%s", type.c_str() );
|
||||
|
||||
tgShapefile::FromPolygon( clipped, ds_name, layer, "poly" );
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
#endif
|
||||
|
||||
// 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 );
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
//
|
||||
// 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 );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -18,6 +18,7 @@ public:
|
|||
|
||||
private:
|
||||
long int GenerateIndex( std::string path );
|
||||
void ClipRow( const tgPolygon& subject, const double& center_lat, const std::string& type );
|
||||
tgPolygon Clip( const tgPolygon& subject, const std::string& type, SGBucket& b );
|
||||
void Chop( const tgPolygon& subject, const std::string& type );
|
||||
|
||||
|
|
|
|||
|
|
@ -73,6 +73,23 @@ double tgContour::GetArea( void ) const
|
|||
return fabs(area * 0.5);
|
||||
}
|
||||
|
||||
bool tgContour::IsInside( const tgContour& inside, const tgContour& outside )
|
||||
{
|
||||
// first contour is inside second if the intersection of first with second is == first
|
||||
// Intersection returns a polygon...
|
||||
tgPolygon result;
|
||||
bool isInside = false;
|
||||
result = Intersect( inside, outside );
|
||||
|
||||
if ( result.Contours() == 1 ) {
|
||||
if ( result.GetContour(0) == inside ) {
|
||||
isInside = true;
|
||||
}
|
||||
}
|
||||
|
||||
return isInside;
|
||||
}
|
||||
|
||||
bool tgContour::RemoveCycles( const tgContour& subject, tgcontour_list& result )
|
||||
{
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "remove cycles : contour has " << subject.GetSize() << " points" );
|
||||
|
|
@ -116,6 +133,27 @@ bool tgContour::RemoveCycles( const tgContour& subject, tgcontour_list& result )
|
|||
second.AddNode( subject.GetNode(n) );
|
||||
}
|
||||
|
||||
// determine hole vs boundary
|
||||
if ( IsInside( first, second ) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "first contur is within second contour " );
|
||||
|
||||
// first contour is inside second : mark first contour as opposite of subject
|
||||
first.SetHole( !subject.GetHole() );
|
||||
second.SetHole( subject.GetHole() );
|
||||
} else if ( IsInside( second, first ) ) {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "second contur is within first contour " );
|
||||
|
||||
// second contour is inside first : mark second contour as opposite of subject
|
||||
first.SetHole( subject.GetHole() );
|
||||
second.SetHole( !subject.GetHole() );
|
||||
} else {
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "conturs are (mostly) disjoint " );
|
||||
|
||||
// neither contour is inside - bots are the same as subject
|
||||
first.SetHole( subject.GetHole() );
|
||||
second.SetHole( subject.GetHole() );
|
||||
}
|
||||
|
||||
SG_LOG(SG_GENERAL, SG_DEBUG, "remove first: size " << first.GetSize() );
|
||||
first.SetHole( subject.GetHole() );
|
||||
RemoveCycles( first, result );
|
||||
|
|
@ -460,6 +498,32 @@ tgPolygon tgContour::Diff( const tgContour& subject, tgPolygon& clip )
|
|||
return result;
|
||||
}
|
||||
|
||||
tgPolygon tgContour::Intersect( const tgContour& subject, const tgContour& clip )
|
||||
{
|
||||
tgPolygon result;
|
||||
UniqueSGGeodSet all_nodes;
|
||||
|
||||
/* before diff - gather all nodes */
|
||||
for ( unsigned int i = 0; i < subject.GetSize(); ++i ) {
|
||||
all_nodes.add( subject.GetNode(i) );
|
||||
}
|
||||
|
||||
ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject );
|
||||
ClipperLib::Polygon clipper_clip = tgContour::ToClipper( clip );
|
||||
ClipperLib::Polygons clipper_result;
|
||||
|
||||
ClipperLib::Clipper c;
|
||||
c.Clear();
|
||||
c.AddPolygon(clipper_subject, ClipperLib::ptSubject);
|
||||
c.AddPolygon(clipper_clip, ClipperLib::ptClip);
|
||||
c.Execute(ClipperLib::ctIntersection, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd);
|
||||
|
||||
result = tgPolygon::FromClipper( clipper_result );
|
||||
result = tgPolygon::AddColinearNodes( result, all_nodes );
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
static bool FindIntermediateNode( const SGGeod& start, const SGGeod& end,
|
||||
const std::vector<SGGeod>& nodes, SGGeod& result,
|
||||
double bbEpsilon, double errEpsilon )
|
||||
|
|
|
|||
|
|
@ -77,6 +77,25 @@ public:
|
|||
double GetMinimumAngle( void ) const;
|
||||
double GetArea( void ) const;
|
||||
|
||||
bool operator==(const tgContour& other ) {
|
||||
bool isEqual = true;
|
||||
|
||||
if ( GetSize() == other.GetSize() )
|
||||
{
|
||||
for (unsigned int i=0; i<GetSize(); i++) {
|
||||
if ( GetNode(i) != other.GetNode(i) ) {
|
||||
isEqual = false;
|
||||
break;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
isEqual = false;
|
||||
}
|
||||
|
||||
return isEqual;
|
||||
}
|
||||
|
||||
|
||||
static tgContour Snap( const tgContour& subject, double snap );
|
||||
static tgContour RemoveDups( const tgContour& subject );
|
||||
static tgContour SplitLongEdges( const tgContour& subject, double dist );
|
||||
|
|
@ -85,7 +104,9 @@ public:
|
|||
|
||||
static tgPolygon Union( const tgContour& subject, tgPolygon& clip );
|
||||
static tgPolygon Diff( const tgContour& subject, tgPolygon& clip );
|
||||
static tgPolygon Intersect( const tgContour& subject, const tgContour& clip );
|
||||
|
||||
static bool IsInside( const tgContour& inside, const tgContour& outside );
|
||||
static tgContour AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes );
|
||||
static tgContour AddColinearNodes( const tgContour& subject, std::vector<SGGeod>& nodes );
|
||||
static bool FindColinearLine( const tgContour& subject, const SGGeod& node, SGGeod& start, SGGeod& end );
|
||||
|
|
|
|||
|
|
@ -63,6 +63,10 @@ string attribute_query;
|
|||
bool use_spatial_query=false;
|
||||
double spat_min_x, spat_min_y, spat_max_x, spat_max_y;
|
||||
int num_threads = 1;
|
||||
bool save_shapefiles=false;
|
||||
std::string ds_name=".";
|
||||
|
||||
const double gSnap = 0.00000001; // approx 1 mm
|
||||
|
||||
SGLockedQueue<OGRFeature *> global_workQueue;
|
||||
|
||||
|
|
@ -197,6 +201,7 @@ void Decoder::processPolygon(OGRPolygon* poGeometry, const string& area_type )
|
|||
|
||||
// first add the outer ring
|
||||
tgPolygon shape = tgShapefile::ToPolygon( poGeometry );
|
||||
shape = tgPolygon::Simplify( shape );
|
||||
|
||||
if ( max_segment_length > 0 ) {
|
||||
shape = tgPolygon::SplitLongEdges( shape, max_segment_length );
|
||||
|
|
@ -608,6 +613,10 @@ int main( int argc, char **argv ) {
|
|||
num_threads=boost::thread::hardware_concurrency();
|
||||
argv+=1;
|
||||
argc-=1;
|
||||
} else if (!strcmp(argv[1],"--debug")) {
|
||||
argv++;
|
||||
argc--;
|
||||
save_shapefiles=true;
|
||||
} else if (!strcmp(argv[1],"--help")) {
|
||||
usage(progname);
|
||||
} else {
|
||||
|
|
|
|||
Loading…
Reference in a new issue