[Gdal-dev] Cookie cutting shapefiles in OGR

Alexandre Leroux alexandre.leroux at ec.gc.ca
Tue Jul 4 13:31:20 EDT 2006


Thanks a lot Martin & Frank,

It's funny because we coded a solution using GPC. We have a few bugs 
left (in our bindings?) that should be corrected very soon. Thanks a lot 
for your C code Martin, we'll compare it with ours right away! We are 
planning to share our resulting code with the ogr crowd :-)

Cheers!

Alex

Chapman, Martin wrote:
> 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,

-- 
Alexandre Leroux, M.Sc., Ing.
Environnement Canada / Environment Canada
Centre météorologique canadien / Canadian Meteorological Centre
Division de la réponse aux urgences environnementales /
Environmental Emergency Response Division
alexandre.leroux at ec.gc.ca



More information about the Gdal-dev mailing list