Calcular el área común entre polígonos usando python?

Question

Estoy tratando de calcular el común de "vacío" de la zona entre los vecinos de los polígonos con el fin de calcular el volumen de la porosidad (la zona roja en la imagen)

Mi idea era, primero para uso Convertir FeatureToVertices entonces, para calcular la distancia más corta entre un punto al siguiente polígono y, a continuación, tratando de generar los nuevos puntos (puntos azules) y, finalmente, para generar el común de polígonos a partir de estos puntos. Yo no podía encontrar la manera de escribir el código y que sólo han logrado hacer los dos primeros pasos. Alguna idea de cómo manejar este problema?

Answer 1

Tal vez usted podría crear un nuevo polígono que encierra a todos los demás (o en un sub-grupo de los polígonos existentes, por ejemplo, los de un lado de la carretera) utilizando el convexo, o mejor, el casco cóncavo y, a continuación, la diferencia de la original de los polígonos de este derivado del polígono. Esto le puede dar una medida de la porosidad.

EDITAR Bueno, he escrito un pequeño script en el libre y open-source GIS Whitebox las Herramientas de Análisis Geoespacial (disponible para su descarga aquí) que identificará el espacio vacío entre la construcción de los polígonos. Funciona mediante la ampliación de la línea de los segmentos que componen la construcción de polígonos, entonces polygonizing el grupo de extender los segmentos de línea. Entonces se disuelve la resultante de múltiples polígono y, finalmente, las diferencias en el resultado de los edificios originales. No es perfecto, pero como mencioné en mi comentario, no creo que hay un consistente conjunto de reglas que se pueden replicar lo que han hecho manualmente en la pregunta anterior. Aún así, creo que el resultado puede ser satisfactorio para el propósito.

Nota, el guión está escrito en Groovy y se puede ejecutar desde el Whitebox GAT scripter.

import whitebox.geospatialfiles.ShapeFile
import whitebox.geospatialfiles.shapefile.*
import whitebox.geospatialfiles.shapefile.attributes.*
import whitebox.utilities.FileUtilities
import whitebox.utilities.Topology
import com.vividsolutions.jts.operation.polygonize.Polygonizer
import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.LineString
import com.vividsolutions.jts.geom.Geometry
import com.vividsolutions.jts.geom.GeometryCollection
import com.vividsolutions.jts.geom.GeometryFactory
import com.vividsolutions.jts.geom.MultiPolygon
import com.vividsolutions.jts.geom.Polygon
import com.vividsolutions.jts.operation.polygonize.Polygonizer
import com.vividsolutions.jts.geom.impl.CoordinateArraySequence;

com.vividsolutions.jts.geom.Geometry g
com.vividsolutions.jts.geom.Geometry buildings
double x1, x2, stX, endX, y1, y2, stY, endY, slope
GeometryFactory factory = new GeometryFactory();
Polygonizer polygonizer = new Polygonizer();

// Change the following three values as needed
String inputFile = pluginHost.getWorkingDirectory() + "building footprints.shp"
String outputFile = pluginHost.getWorkingDirectory() + "tmp5.shp"
int extensionDistance = 70

ShapeFile input = new ShapeFile(inputFile)
ArrayList<com.vividsolutions.jts.geom.Polygon> buildingPolygonList = new ArrayList<>();
ArrayList<com.vividsolutions.jts.geom.Geometry> inputGeometryList = new ArrayList<>();
boolean firstFeature = true
for (ShapeFileRecord record : input.records) {
    double[][] points = record.getGeometry().getPoints()
    recJTSGeometries = record.getGeometry().getJTSGeometries()
    for (Geometry k : recJTSGeometries) {
        buildingPolygonList.add((com.vividsolutions.jts.geom.Polygon)k)
    }
    for (int p = 0; p < points.length - 1; p++) {
        x1 = points[p][0]
        y1 = points[p][1]
        x2 = points[p + 1][0]
        y2 = points[p + 1][1]

        // Extend this line segment on both sides by extensionDistance
        if (x1 - x2 != 0) {
            slope = Math.atan2((y2 - y1) , (x2 - x1))
            xSt = x1 - extensionDistance * Math.cos(slope)
            ySt = y1 - extensionDistance * Math.sin(slope)
            xEnd = x2 + extensionDistance * Math.cos(slope)
            yEnd = y2 + extensionDistance * Math.sin(slope)
        } else {
            xSt = x1
            xEnd = x2
            if (y2 > y1) {
                ySt = y1 - extensionDistance
            } else {
                ySt = y1 + extensionDistance
            }
        }
        CoordinateArraySequence coordArray = new CoordinateArraySequence(2);
        coordArray.setOrdinate(0, 0, xSt);
        coordArray.setOrdinate(0, 1, ySt);
        coordArray.setOrdinate(1, 0, xEnd);
        coordArray.setOrdinate(1, 1, yEnd);

        if (!firstFeature) {
            inputGeometryList.add(factory.createLineString(coordArray))
        } else {
            g = factory.createLineString(coordArray)
            firstFeature = false
        }
    }
}

// make the buildings geometry
com.vividsolutions.jts.geom.Polygon[] polygonArray = new com.vividsolutions.jts.geom.Polygon[buildingPolygonList.size()];
for (i = 0; i < buildingPolygonList.size(); i++) {
    polygonArray[i] = buildingPolygonList.get(i);
}

buildings = factory.createMultiPolygon(polygonArray);
buildings = buildings.buffer(0) //corrects the topology in layer

g = g.union(factory.buildGeometry(inputGeometryList))
polygonizer.add(g)
ArrayList<com.vividsolutions.jts.geom.Polygon> outputPoly = (ArrayList)polygonizer.getPolygons()

int FID = 0
if (outputPoly.size() > 0) {
    // perform a dissolve on the resulting polygons
    polygonArray = new com.vividsolutions.jts.geom.Polygon[outputPoly.size()];
    for (i = 0; i < outputPoly.size(); i++) {
        polygonArray[i] = outputPoly.get(i);
    }

    geometriesToBuffer = factory.createMultiPolygon(polygonArray);

    Geometry buffer = geometriesToBuffer.buffer(0)

    buffer = buffer.difference(buildings)

    // output the file
    DBFField[] fields = new DBFField[1];

    fields[0] = new DBFField();
    fields[0].setName("FID");
    fields[0].setDataType(DBFField.DBFDataType.NUMERIC);
    fields[0].setFieldLength(10);
    fields[0].setDecimalCount(0);

    ShapeFile output = new ShapeFile(outputFile, ShapeType.POLYGON, fields);

    if (buffer instanceof com.vividsolutions.jts.geom.MultiPolygon) {
        MultiPolygon mpBuffer = (MultiPolygon) buffer;
        FID = 0;
        n = 0;
        for (int a = 0; a < mpBuffer.getNumGeometries(); a++) {
            g = mpBuffer.getGeometryN(a);
            if (g instanceof com.vividsolutions.jts.geom.Polygon) {
                Polygon poly = (Polygon)(g)
                ArrayList<ShapefilePoint> pnts = new ArrayList<>();

                int[] parts = new int[poly.getNumInteriorRing() + 1];

                Coordinate[] coords = poly.getExteriorRing().getCoordinates();
                if (Topology.isClockwisePolygon(coords)) {
                    for (i = 0; i < coords.length; i++) {
                        pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
                    }
                } else {
                    for (i = coords.length - 1; i >= 0; i--) {
                        pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
                    }
                }

                for (int b = 0; b < poly.getNumInteriorRing(); b++) {
                    parts[b + 1] = pnts.size();
                    coords = poly.getInteriorRingN(b).getCoordinates();
                    if (Topology.isClockwisePolygon(coords)) {
                        for (i = coords.length - 1; i >= 0; i--) {
                            pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
                        }
                    } else {
                        for (i = 0; i < coords.length; i++) {
                            pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
                        }
                    }
                }

                PointsList pl = new PointsList(pnts);
                whitebox.geospatialfiles.shapefile.Geometry wbGeometry = new whitebox.geospatialfiles.shapefile.Polygon(parts, pl.getPointsArray());
                FID++;
                Object[] rowData = new Object[1];
                rowData[0] = new Double(FID);
                output.addRecord(wbGeometry, rowData);
            }
        }
    } else if (buffer instanceof com.vividsolutions.jts.geom.Polygon) {
        Polygon poly = (Polygon)(buffer)
        ArrayList<ShapefilePoint> pnts = new ArrayList<>();

        int[] parts = new int[poly.getNumInteriorRing() + 1];

        Coordinate[] coords = poly.getExteriorRing().getCoordinates();
        if (Topology.isClockwisePolygon(coords)) {
            for (i = 0; i < coords.length; i++) {
                pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
            }
        } else {
            for (i = coords.length - 1; i >= 0; i--) {
                pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
            }
        }

        for (int b = 0; b < poly.getNumInteriorRing(); b++) {
            parts[b + 1] = pnts.size();
            coords = poly.getInteriorRingN(b).getCoordinates();
            if (Topology.isClockwisePolygon(coords)) {
                for (i = coords.length - 1; i >= 0; i--) {
                    pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
                }
            } else {
                for (i = 0; i < coords.length; i++) {
                    pnts.add(new ShapefilePoint(coords[i].x, coords[i].y));
                }
            }
        }

        PointsList pl = new PointsList(pnts);
        whitebox.geospatialfiles.shapefile.Geometry wbGeometry = new whitebox.geospatialfiles.shapefile.Polygon(parts, pl.getPointsArray());
        FID++;
        Object[] rowData = new Object[1];
        rowData[0] = new Double(FID);
        output.addRecord(wbGeometry, rowData);
    }

    output.write()

    // display the output image
    pluginHost.returnData(outputFile)

    println("Operation complete")
} else {
    println("Something went wrong.")
}

El guión no es bonito, y estoy seguro de que podría haber limpiado bastante, pero funciona.