Konvertieren Sie den Polygon-Feature-Schwerpunkt mithilfe von Python in Punkte

9

Ich möchte einige polygonbasierte shp-Dateien mit mehreren Polygon-Features in Punkte für jedes Feature konvertieren, die im Wesentlichen den Mittelpunkt jedes Polygon-Features darstellen. Ich weiß, dass ich in der ArcGIS-Welt das Feature-To-Point-Tool verwenden kann, aber ich möchte dies in einem Skript behalten, das auf PCs ausgeführt werden kann, auf denen kein Arcpy vorhanden ist. Daher suche ich nach einer Open-Source-Alternative dazu. Ist jemandem eine Bibliothek bekannt, die ich dafür verwenden könnte, zusammen mit einer Anleitung, wie man sie nutzen kann, um dies zu erreichen?

wilbev
quelle
Ich habe immer noch mehrere Probleme mit der Antwort, die Gene unten gegeben hat. Die Probleme sind, wie die Attribute von ihrer ursprünglichen Reihenfolge in alphabetische Reihenfolge umgeordnet werden, was ein Problem darstellt. Zweitens wird die Formdatei möglicherweise aufgrund der Datei, die ich konvertieren möchte, mit über 250 Attributen beschädigt.
Wilbev
In QGIS gibt es ein Standardwerkzeug namens "Polygon Centroids", das genau dies tut - benötigen Sie ein Skript? Ich denke, es wäre einfach genug, Skripte mit PyQGIS zu schreiben.
Dùn Caan
Es muss ein Skript sein und auf PCs funktionieren, auf denen kein QGIS installiert ist.
Wilbev

Antworten:

9

Sie können einen ogr2ogrBefehl ausführen (z. B. über eine OSGeo4w-Shell). ZB auf einem Shapefile von Ländern:

cd path/to/shapefiles
ogr2ogr -sql "SELECT ST_Centroid(geometry), * FROM countries" -dialect sqlite countries_centroid.shp countries.shp

Das neue Shapefile countries_centroid.shpsollte der Eingabe ähnlich sein, jedoch nur einen Punkt pro [Multi] -Polygon enthalten.

@PEL zeigt auch ein gutes Beispiel mit ST_PointOnSurface, das in diesem Befehl einfach zu ersetzen ist.


Ähnliches kann bei Bedarf in Python getan werden, es kann jedoch einige Codezeilen länger dauern:

import os
from osgeo import ogr

ogr.UseExceptions()
os.chdir('path/to/shapefiles')

ds = ogr.Open('countries.shp')
ly = ds.ExecuteSQL('SELECT ST_Centroid(geometry), * FROM countries', dialect='sqlite')
drv = ogr.GetDriverByName('Esri shapefile')
ds2 = drv.CreateDataSource('countries_centroid.shp')
ds2.CopyLayer(ly, '')
ly = ds = ds2 = None  # save, close
Mike T.
quelle
Ich denke, Sie haben eine einfachste Lösung mit OGR und SQL vorgeschlagen. Ich denke, es ist sicherer, einen Parameter zu OGR mit -nlt Point
PEL
Leider kann ich das nicht zum Laufen bringen. Ich erhalte eine Fehlermeldung, dass ST_Centroid nicht funktioniert.
Wilbev
1
Es benötigt die in GDAL integrierte SQLite-Dialektoption (wie gezeigt) und Spatialite, was nicht immer garantiert ist. OSGeo4W hat eine gute Version von GDAL, die diesen Befehl korrekt ausführt.
Mike T
Ich kann Ihr Top-Skript innerhalb von ogr2ogr empfehlen, ohne Probleme zu arbeiten. Ich muss dies jedoch in einem eigenständigen Python-Skript tun, damit ich versuche, Ihren zweiten Codesatz zum Laufen zu bringen. Dort, wo ich mit ST_Centroid fortfahre, ist kein Funktionsfehler. Mein Code ist identisch mit dem, was Sie oben haben, einschließlich des SQLite-Dialekts.
Wilbev
1
Der von Ihnen beschriebene Fehler tritt auf, wenn GDAL ohne Spatialite-Unterstützung erstellt wurde. Einige gdal-python-Pakete unterstützen dies, aber nicht alle. Versuchen Sie, eine OSGeo4W-Shell zu öffnen und das Python-Skript in dieser Umgebung auszuführen. Ich denke, die Standardpakete beziehen sich auf gdal-bindiese Unterstützung.
Mike T
9

Verwenden Sie einfach Fiona oder GeoPandas (Python 2.7.x und 3.x)

Einige Polygone

Geben Sie hier die Bildbeschreibung ein

import geopandas as gpd
# GeoDataFrame creation
poly = gpd.read_file("geoch_poly.shp")
poly.head()

Geben Sie hier die Bildbeschreibung ein

Transformation zu Punkten (Zentroiden)

# copy poly to new GeoDataFrame
points = poly.copy()
# change the geometry
points.geometry = points['geometry'].centroid
# same crs
points.crs =poly.crs
points.head()

Geben Sie hier die Bildbeschreibung ein

# save the shapefile
points.to_file('geoch_centroid.shp')

Ergebnis

Geben Sie hier die Bildbeschreibung ein

Gen
quelle
Danke für das Antwortgen. Ich denke, Sie haben möglicherweise einen Tippfehler oben, bei dem die Variable 'gdf' im Code points.crs = gdf.crs 'poly' sein sollte. Ich habe auch einige andere Probleme, bei denen die PRJ-Datei nicht abgerufen wird erstellt, wird es leer angezeigt und die Reihenfolge der Attributfelder ändert ihre Reihenfolge gegenüber den Polygondaten, da sie jetzt alphabetisch zu sein scheinen. Es ist wichtig, dass sie in derselben Reihenfolge bleiben. Sie kennen eine Möglichkeit, das Feldattribut in derselben Reihenfolge zu halten bestellen?
wilbev
Danke, korrigiert. Für die Reihenfolge der Attributfelder ändern Sie einfach die Reihenfolge des GeoPandas GeoDataFrame (= Pandas DataFrame)
Gens
Danke Gene, aber ich bin mir nicht sicher, wo ich diese Reihenfolge im Code hier ändern würde. Ich stoße auch auf zwei andere Probleme. Zuerst ist die * .prj-Datei in der neuen shp-Datei leer. Zweitens, wenn ich versuche, die shp-Datei im shp-Reader zu öffnen, tritt beim Öffnen der Datei ein Fehler auf, als wäre sie beschädigt. Es scheint zu funktionieren, ohne beschädigt zu werden, wenn die shp-Datei nur eine einzige Funktion hat, aber bei mehreren scheint es Probleme zu geben.
Wilbev
Entschuldigung, aber Sie müssen Pandas dafür kennen und ich habe kein Problem mit dem Skript mit meinen Daten (ich verwende die neuesten Versionen von GeoPandas, Fiona und Numpy)
Gen
Ich kann Ihnen die shp-Datei senden, damit Sie selbst sehen können, aber diese shp-Datei enthält über 250 Datenspalten, von denen ich mir vorstelle, dass sie Probleme verursachen. Ich habe dies an einer shp-Datei mit viel weniger Attributen versucht und es scheint keine Probleme zu geben.
Wilbev
5

Ein anderer, vielleicht "niedrigerer" Weg wäre die direkte Verwendung fionaund shapelyfür die E / A- und Geometrieverarbeitung.

import fiona
from shapely.geometry import shape, mapping

with fiona.open('input_shapefile.shp') as src:
    meta = src.meta
    meta['schema']['geometry'] = 'Point'
    with fiona.open('output_shapefile.shp', 'w', **meta) as dst:
        for f in src:
            centroid = shape(f['geometry']).centroid
            f['geometry'] = mapping(centroid)
            dst.write(f)
Loïc Dutrieux
quelle
Danke Loic. Das behebt definitiv das Sortierproblem, das ich hatte, aber es behebt das Problem nicht mit so vielen Attributen, dass die Datei beschädigt wird. Haben Sie weitere Ideen, wie Sie dieses Problem lösen können? Ich denke, ich muss möglicherweise Attribute löschen. Ich kann Ihnen eine Beispieldatei senden, wenn es helfen würde.
Wilbev
@wilbev Senden Sie einen Download-Link zu Ihren Daten, wenn Sie können. Sonst sehe ich nicht, was falsch sein könnte.
Loïc Dutrieux
Loic, ich habe dir eine Beispieldatei per E-Mail geschickt. Hoffentlich gibt Ihnen das eine gute Vorstellung von dem Problem, auf das ich stoße.
Wilbev
@wilbev was meinst du mit 'die Datei wird beschädigt'? Mit der von Ihnen gesendeten Datei konnte ich die Schwerpunkte erstellen und das Ausgabe-Shapefile in QGIS problemlos öffnen. Die Attributtabelle bleibt zwischen beiden Dateien unverändert.
Loïc Dutrieux
Mit "korrupt" meine ich, dass es sich im Wesentlichen um eine leere DBF-Datei handelt, da nach dem Ausführen des Skripts eine DBF-Datei mit einer Größe von 1 KB erstellt wird und diese beim Öffnen vollständig leer ist. Wenn ich dasselbe Skript wie das oben aufgeführte in einer Datei mit weniger Attributen ausführe, funktioniert es ohne Probleme. Ich habe sogar einen zweiten PC ausprobiert und das gleiche Ergebnis erzielt. Ich verstehe es nicht
Wilbev
2

Ich denke, der einfachste Weg ist die Verwendung des virtuellen Formats gdal / ogr. ( http://www.gdal.org/drv_vrt.html ) und SQL / SQLITE-Dialekt ( http://www.gdal.org/ogr_sql.html und https://www.gaia-gis.it/spatialite-3.0) .0-BETA / Spatialite-SQL-3.0.0.html )

Mein Polygon-Shapefile heißt poly.shp. Dann erstelle ich diese XML-ähnliche Datei mit dem Namen vrt.vrt. In dieser Datei (vrt.vrt) finden Sie hier den Inhalt, der in Punkte konvertiert werden soll

<OGRVRTDataSource>
    <OGRVRTLayer name="poly">
        <SrcDataSource relativeToVRT="1">poly.shp</SrcDataSource>
        <SrcSQL dialect="sqlite">SELECT ST_PointOnSurface(geometry) as geom_point, poly.* from poly</SrcSQL>
        <GeometryType>wkbPoints</GeometryType> 
        <GeometryField name="geom_point" />
    </OGRVRTLayer>
</OGRVRTDataSource>

Zu diesem Zeitpunkt können Sie diese Datei zur Validierung in Qgis integrieren. Das Rendering ist sicher langsamer als die Rohquelle, da jedes Feature bei jeder Rendering-Abfrage als Punkt gewertet wird.

Konvertieren Sie anschließend diese Datei (vrt.vrt) mit gdal / ogr-Utils aus einer Python-Shell / einem Python-Skript in eine andere Datei

os.system("ogr2ogr point_from_vrt.shp vrt.vrt poly")

Sie erhalten ein Punkt-Shapefile mit dem Namen point_from_vrt.shp.

PEL
quelle
Ich konnte dies zum Laufen bringen, aber ich möchte all dies in einem Python-Skript behalten, da ich Hunderte von Dateien mit unterschiedlichen Dateinamen konvertieren muss. Ich möchte die @ Mike T-Lösung verwenden, erhalte jedoch die Meldung "Keine solche Funktion: ST_Centroid, wenn ich diese verwende, und ich habe auch ST_PointOnSurface ausprobiert, das auch besagt, dass es keine solche Funktion gibt. Eine Idee, warum angegeben wird, dass dies keine Funktionen von ExecuteSQL sind." ()?
Wilbev
Ich bekomme'wkbPoints' is not a valid value of the atomic type
Ben Sinclair