[Gdal-dev] Cookie cutting shapefiles in OGR
Chapman, Martin
MChapman at sanz.com
Tue Jul 4 13:09:12 EDT 2006
Alex,
There are two ways you can do this. I'm not sure what you want to do:
1. You can do the overlaps (bbox) test first using OGR (without GEOS)
as Frank suggested, to eliminate all geometries who's envelopes don't
overlap, then you do another pass through the remaining geometries and
do a strict intersection test to eliminate those geometries that don't
exactly intersect. I wrote the code at the bottom in java to do this
and I have the c++ port around somewhere if you want it, but this code
can be ported easily. Point and LineString intersections can use the
same...slightly modified isPolyInPoly()...maybe create a
isLineStringInPoly() and isPointInPoly() function and so on to handle
all your geometry type cases. By the way...it's very fast.
2. If you actually want to clip intersecting geometries and reshape
them to your cookie cutter geometry then you need to use a library like
GPC (http://www.cs.man.ac.uk/~toby/alan/software/ ). This open source
library is very simple to use and works great...and it's very fast. If
you have problems using it let me know and I will help you out. If
Frank W. is reading this posting I would urge him to consider including
GPC as part of OGR since it's small, fast, powerful, and easy to use.
LineString and Point clipping is easy and you should be able to find
algorithms on-line or come up with one yourself. GPC also supports
difference, intersection, exclusive-or and union clip operations. The
C++ code I use for GPC is listed below the java code.
Best regards,
Martin
//////////////////////////JAVA
private boolean isPolyInPoly(double[] X1Points,
double[] Y1Points,
double[] X2Points,
double[] Y2Points)
{
try
{
int X1Length = X1Points.length;
int X2Length = X2Points.length;
for (int i = 0; i < X1Length; i++)
{
if (isPointInPoly(X2Length, X2Points, Y2Points,
X1Points[i], Y1Points[i]))
{
return true;
}
}
for (int j = 0; j < X2Length; j++)
{
if (isPointInPoly(X1Length, X1Points, Y1Points,
X2Points[j], Y2Points[j]))
{
return true;
}
}
Line2D.Double lineSegment = new Line2D.Double();
for (int k = 0; k < X1Points.length - 1; k++)
{
lineSegment.x1 = X1Points[k];
lineSegment.y1 = Y1Points[k];
lineSegment.x2 = X1Points[k + 1];
lineSegment.y2 = Y1Points[k + 1];
if (isLineInPoly(lineSegment, X2Points, Y2Points))
{
return true;
}
}
}
catch (Exception e)
{
//
}
return false;
}
private boolean isPointInPoly(int polySize, double[] XPoints,
double[] YPoints, double x, double y)
{
try
{
int j = polySize - 1;
boolean isInside = false;
for (int i = 0; i < polySize; j = i++)
{
if ( ( ( (YPoints[i] <= y) && (y < YPoints[j])) ||
( (YPoints[j] <= y) && (y < YPoints[i]))) &&
(x < (XPoints[j] - XPoints[i]) *
(y - YPoints[i]) / (YPoints[j] - YPoints[i]) +
XPoints[i]))
{
isInside = !isInside;
}
}
return isInside;
}
catch (Exception e)
{
//
}
return false;
}
private boolean isLineInPoly(Line2D.Double pL, double[] XPoints,
double[] YPoints)
{
try
{
if (isPointInPoly(XPoints.length, XPoints, YPoints, pL.x1,
pL.y1))
{
return true;
}
if (isPointInPoly(XPoints.length, XPoints, YPoints, pL.x2,
pL.y2))
{
return true;
}
Line2D.Double lineSegment = new Line2D.Double();
for (int i = 0; i < XPoints.length - 1; i++)
{
lineSegment.x1 = XPoints[i];
lineSegment.y1 = YPoints[i];
lineSegment.x2 = XPoints[i + 1];
lineSegment.y2 = YPoints[i + 1];
if (isLineInLine(pL, lineSegment))
{
return true;
}
}
}
catch (Exception e)
{
//
}
return false;
}
private boolean isLineInLine(Line2D.Double pL1, Line2D.Double pL2)
{
try
{
if ( (pL1.x1 == pL2.x1 && pL1.y1 == pL2.y1) || (pL1.x1 ==
pL2.x2 && pL1.y1 == pL2.y2) ||
(pL1.x2 == pL2.x1 && pL1.y2 == pL2.y1) || (pL1.x2 ==
pL2.x2 && pL1.y2 == pL2.y2))
{
return true;
}
isAngleCounterClockWise( (Point2D.Double) pL1.getP1(),
(Point2D.Double) pL1.getP2(), (Point2D.Double) pL2.getP1());
return (boolean) ( (isAngleCounterClockWise(
(Point2D.Double) pL1.getP1(), (Point2D.Double) pL1.getP2(),
(Point2D.Double) pL2.getP1()) !=
isAngleCounterClockWise(
(Point2D.Double) pL1.getP1(), (Point2D.Double) pL1.getP2(),
(Point2D.Double) pL2.getP2())) &&
(isAngleCounterClockWise( (Point2D.Double)
pL2.getP1(), (Point2D.Double) pL2.getP2(), (Point2D.Double) pL1.getP1())
!=
isAngleCounterClockWise( (Point2D.Double)
pL2.getP1(), (Point2D.Double) pL2.getP2(), (Point2D.Double)
pL1.getP2())));
}
catch (Exception e)
{
//
}
return false;
}
private boolean isAngleCounterClockWise(Point2D.Double p1,
Point2D.Double p2, Point2D.Double p3)
{
try
{
double dx1 = 0.0, dx2 = 0.0, dy1 = 0.0, dy2 = 0.0;
dx1 = p2.x - p1.x;
dy1 = p2.y - p1.y;
dx2 = p3.x - p2.x;
dy2 = p3.y - p2.y;
if ( (dx1 == 0 && dx2 == 0) || (dy1 == 0 && dy2 == 0))
{
return false;
}
if (dy1 * dx2 < dy2 * dx1)
{
return true;
}
else
{
return false;
}
}
catch (Exception e)
{
//
}
return false;
}
//////////////////////////////////////// GPC
void CVPAlgorithms::GetIntersection(OGRGeometry* pGeometryA,
OGRGeometry* pGeometryB, OGRGeometry* pResult)
{
try
{
if (!pGeometryA || !pGeometryB || !pResult) return;
if (pGeometryA->Intersects(pGeometryB) == FALSE) return;
gpc_polygon* pPolygonA = NULL;
gpc_polygon* pPolygonB = NULL;
// convert geometry A
OGRwkbGeometryType geomTypeA =
pGeometryA->getGeometryType();
if (geomTypeA == wkbPoint)
{
}
else if (geomTypeA == wkbLineString)
{
}
else if (geomTypeA == wkbPolygon)
{
OGRLinearRing* pLinearRing = NULL;
OGRPolygon* pPolygon = (OGRPolygon*) pGeometryA;
pPolygonA = new gpc_polygon();
pPolygonA->num_contours =
pPolygon->getNumInteriorRings() + 1;
pPolygonA->hole = (int*) new
int[pPolygonA->num_contours];
pPolygonA->contour = (gpc_vertex_list*) new
gpc_vertex_list[pPolygonA->num_contours];
for (int c = 0; c < pPolygonA->num_contours;
c++)
{
if (c == 0)
pLinearRing =
pPolygon->getExteriorRing();
else
pLinearRing =
pPolygon->getInteriorRing(c - 1);
pPolygonA->contour[c].num_vertices =
pLinearRing->getNumPoints();
pPolygonA->hole[c] = c == 0 ? FALSE :
TRUE;
pPolygonA->contour[c].vertex = new
gpc_vertex[pPolygonA->contour[c].num_vertices];
if (c == 0)
{
for (int v = 0; v <
pPolygonA->contour[c].num_vertices; v++)
{
pPolygonA->contour[c].vertex[v].x = pLinearRing->getX(v);
pPolygonA->contour[c].vertex[v].y = pLinearRing->getY(v);
}
}
else
{
for (int v =
pPolygonA->contour[c].num_vertices - 1; v >= 0; v--)
{
pPolygonA->contour[c].vertex[v].x = pLinearRing->getX(v);
pPolygonA->contour[c].vertex[v].y = pLinearRing->getY(v);
}
}
}
}
else if (geomTypeA == wkbMultiPoint)
{
}
else if (geomTypeA == wkbMultiLineString)
{
}
else if (geomTypeA == wkbMultiPolygon)
{
}
else
throw "Unsupported geometry type for parameter
1.";
// convert geometry B
OGRwkbGeometryType geomTypeB =
pGeometryB->getGeometryType();
if (geomTypeB == wkbPoint)
{
}
else if (geomTypeB == wkbLineString)
{
}
else if (geomTypeB == wkbPolygon)
{
OGRLinearRing* pLinearRing = NULL;
OGRPolygon* pPolygon = (OGRPolygon*) pGeometryB;
pPolygonB = new gpc_polygon();
pPolygonB->num_contours =
pPolygon->getNumInteriorRings() + 1;
pPolygonB->hole = (int*) new
int[pPolygonB->num_contours];
pPolygonB->contour = (gpc_vertex_list*) new
gpc_vertex_list[pPolygonB->num_contours];
for (int c = 0; c < pPolygonB->num_contours;
c++)
{
if (c == 0)
pLinearRing =
pPolygon->getExteriorRing();
else
pLinearRing =
pPolygon->getInteriorRing(c - 1);
pPolygonB->contour[c].num_vertices =
pLinearRing->getNumPoints();
pPolygonB->hole[c] = c == 0 ? FALSE :
TRUE;
pPolygonB->contour[c].vertex = new
gpc_vertex[pPolygonB->contour[c].num_vertices];
if (c == 0)
{
for (int v = 0; v <
pPolygonB->contour[c].num_vertices; v++)
{
pPolygonB->contour[c].vertex[v].x = pLinearRing->getX(v);
pPolygonB->contour[c].vertex[v].y = pLinearRing->getY(v);
}
}
else
{
for (int v =
pPolygonB->contour[c].num_vertices - 1; v >= 0; v--)
{
pPolygonB->contour[c].vertex[v].x = pLinearRing->getX(v);
pPolygonB->contour[c].vertex[v].y = pLinearRing->getY(v);
}
}
}
}
else if (geomTypeB == wkbMultiPoint)
{
}
else if (geomTypeB == wkbMultiLineString)
{
}
else if (geomTypeB == wkbMultiPolygon)
{
}
else
throw "Unsupported geometry type for parameter
2.";
// polygon to polygon
if (geomTypeA == wkbPolygon && geomTypeB == wkbPolygon)
{
gpc_polygon* pGPCPolygon = new gpc_polygon();
gpc_polygon_clip(GPC_INT, pPolygonA, pPolygonB,
pGPCPolygon);
OGRPolygon* pOgrPolygon = (OGRPolygon*) pResult;
pOgrPolygon->assignSpatialReference(pGeometryA->getSpatialReference());
gpc_vertex_list* pVertexList = NULL;
for (int i = 0; i < pGPCPolygon->num_contours;
i++)
{
OGRLinearRing ogrLinearRing;
pVertexList = pGPCPolygon->contour;
int nNumPoints =
pVertexList->num_vertices;
for (int n = 0; n < nNumPoints; n++)
{
double nX =
pVertexList->vertex[n].x;
double nY =
pVertexList->vertex[n].y;
ogrLinearRing.addPoint(nX, nY);
}
pOgrPolygon->addRing(&ogrLinearRing);
}
if (pGPCPolygon)
{
for (int c = 0; c <
pGPCPolygon->num_contours; c++)
delete []
pGPCPolygon->contour[c].vertex;
delete [] pGPCPolygon->hole;
delete [] pGPCPolygon->contour;
delete pGPCPolygon;
pGPCPolygon = NULL;
}
}
// clean up
if (pPolygonB)
{
for (int c = 0; c < pPolygonB->num_contours;
c++)
delete [] pPolygonB->contour[c].vertex;
delete [] pPolygonB->hole;
delete [] pPolygonB->contour;
delete pPolygonB;
pPolygonB = NULL;
}
if (pPolygonA)
{
for (int c = 0; c < pPolygonA->num_contours;
c++)
delete [] pPolygonA->contour[c].vertex;
delete [] pPolygonA->hole;
delete [] pPolygonA->contour;
delete pPolygonA;
pPolygonA = NULL;
}
}
catch(exception e)
{throw e;}
catch(char* e)
{throw e;}
catch(...)
{throw "An unexpected error occurred in
CVPAlgorithms::GetIntersection().";}
}
-----Original Message-----
From: gdal-dev-bounces at lists.maptools.org
[mailto:gdal-dev-bounces at lists.maptools.org] On Behalf Of Frank
Warmerdam
Sent: Sunday, July 02, 2006 2:19 PM
To: Alexandre Leroux
Cc: gdal-dev at lists.maptools.org
Subject: Re: [Gdal-dev] Cookie cutting shapefiles in OGR
Alexandre Leroux wrote:
>
> Hi list,
>
> I haven't got any news from my cookie cutting challenge below. Anyone
> has a clue to help me get started?
>
> Thanks and cheers! :-)
>
> Alex
>
>
>> Can someone tell me if I can do the following with OGR: I have two
>> shapefiles and need to cookie cut one with the other while keeping
the
>> attributes of the first shapefile in the newly generated third
shapefile.
>>
>> I saw there is OGRGeometry::Intersection which would allow me to do
>> this, but I'd need to verify the intersection of all features from
>> shapefile #1 with all features of shapefile #2, which would take days
>> and days to compute with OGR (that's what my collegue which knows OGR
>> way better told me). Is there a way to cookie cut shapefiles directly
>> at the layer level?
Alex,
There is no built in support to do this conveniently. You could do an
intersection of all features in layer one against each feature in layer
two, with a pretest to avoid intersection operations on features whose
BBOXes don't even intersect. Even better would be to keep some sort of
spatial index in memory so you would only do bbox comparisons with
features that are relatively nearby. But I will also warn that the
Intersection() operator is pretty expensive, partly because each OGR
geometry must be turned into a GEOS geometry internally. So another
optimization would be to translate all your geometries into GEOS
geometries
once, and do the intersection directly with geos. That would avoid a
lot
of re-conversion compared to a more ogr based approach.
But basically, any approach with OGR and GEOS is going to involve quite
a bit of scripting or programming work on your part.
You might want to look at options for doing this in GRASS. I believe it
has some sort of layer overlay operation (perhaps done in the raster
domain?).
Best regards,
--
---------------------------------------+--------------------------------
------
I set the clouds in motion - turn up | Frank Warmerdam,
warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | President OSGF,
http://osgeo.org
_______________________________________________
Gdal-dev mailing list
Gdal-dev at lists.maptools.org
http://lists.maptools.org/mailman/listinfo/gdal-dev
More information about the Gdal-dev
mailing list