[Gdal-dev] Cookie cutting shapefiles in OGR
Gauthier,Jean-Philippe [CMC]
Jean-Philippe.Gauthier at ec.gc.ca
Tue Jul 4 13:56:17 EDT 2006
There is only one twist with the GPC library, it sometimes return the results of the operations with no specific hole/extern ring order. In specific case when the result of the operation returns multiple polygon with holes, this might be a problem when converting back to OGR structures.
This is something I'm looking in right now.
================================================================
Jean-Philippe Gauthier Bilodeau Direction des opérations
Analyste Programmeur Scientifique Centre Météorologique Canadien
Tel : (514) 421-4642 2121 Route Transcanadienne
Fax: (514) 421-4679 Dorval,Québec
Email: Jean-Philippe.Gauthier at ec.gc.ca H9P 1J3
-----Original Message-----
From: Leroux,Alexandre [CMC]
Sent: 04 July, 2006 13:31
To: Chapman, Martin
Cc: Frank Warmerdam; gdal-dev at lists.maptools.org; Gauthier,Jean-Philippe [CMC]
Subject: Re: [Gdal-dev] Cookie cutting shapefiles in OGR
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())
> pOgrPolygon->;
>
> 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