Convertir le centre de gravité d'une entité polygonale en points à l'aide de Python

Je voudrais convertir certains fichiers shp basés sur des polygones qui ont plusieurs entités polygonales en points pour chaque entité qui représenteraient essentiellement le centriode de chaque entité polygonale. Je sais que dans le monde ArcGIS, je pourrais utiliser l'outil Feature To Point, mais je voudrais le conserver dans un script qui pourrait être exécuté sur des PC sans arcpy, donc je cherche une alternative open source à cela. Quelqu'un est-il au courant d'une bibliothèque que je pourrais utiliser pour cela, ainsi que des instructions sur la façon de l'utiliser pour y parvenir?

Vous pouvez exécuter une ogr2ogrcommande (par exemple à partir d'un shell OSGeo4w). Par exemple sur un fichier de formes de pays:

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

Le nouveau fichier de formes countries_centroid.shpdoit être similaire à l'entrée, mais contenir juste un point par [Multi] Polygone.

@PEL montre également un bon exemple avec ST_PointOnSurface, qui est simple à remplacer dans cette commande.

Quelque chose de similaire peut être fait en Python, si nécessaire, mais cela peut prendre quelques lignes de code de plus:

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

Utilisez simplement Fiona ou GeoPandas (Python 2.7.x et 3.x)

Quelques polygones

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

Transformation en points (centroïdes)

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

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

Résultat

Une autre façon, peut-être plus «bas niveau», serait d'utiliser directement fionaet shapelypour le traitement des E / S et de la géométrie.

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)

Je pense que la façon la plus simple est d'utiliser le format virtuel gdal / ogr. ( http://www.gdal.org/drv_vrt.html ) et dialecte SQL / SQLITE ( http://www.gdal.org/ogr_sql.html et https://www.gaia-gis.it/spatialite-3.0 .0-BETA / spatialite-sql-3.0.0.html )

Mon fichier de formes polygonales s'appelle poly.shp. Ensuite, je crée ce fichier XML appelé vrt.vrt. A l'intérieur de ce fichier (vrt.vrt), voici le contenu à convertir en points

<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>

À ce stade, vous pouvez intégrer ce fichier dans Qgis pour le valider. Bien sûr, le rendu est plus lent que la source brute car chaque entité est castée en tant que point sur chaque requête de rendu.

Ensuite, convertissez ce fichier (vrt.vrt) en quelque chose d'autre en utilisant les utilitaires gdal / ogr à partir d'un shell / script python

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

Vous obtenez un fichier de formes de points nommé point_from_vrt.shp.