Convierta el centroide de entidades poligonales en puntos usando Python

Me gustaría convertir algunos archivos shp basados ​​en polígonos que tienen varias características de polígono en puntos para cada característica que representarían esencialmente el centro de cada característica de polígono. Sé que en el mundo de ArcGIS podría usar la herramienta Feature To Point, pero me gustaría mantener esto en una secuencia de comandos que se pueda ejecutar en PC que no tengan arcpy, así que estoy buscando una alternativa de código abierto para eso. ¿Alguien sabe de una biblioteca que podría usar para esto junto con alguna dirección sobre cómo aprovecharlo para lograr esto?

Puede ejecutar un ogr2ogrcomando (por ejemplo, desde un Shell OSGeo4w). Por ejemplo, en un archivo de forma de países:

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

El nuevo archivo de forma countries_centroid.shpdebe ser similar a la entrada, pero solo debe contener un punto por [Multi] Polígono.

@PEL también muestra un buen ejemplo con ST_PointOnSurface, que es fácil de sustituir en este comando.

Se puede hacer algo similar en Python, si es necesario, pero puede tomar algunas líneas de código más:

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

Simplemente use Fiona o GeoPandas (Python 2.7.xy 3.x)

Algunos polígonos

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

Transformación a puntos (centroides)

# 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')

Resultado

Otra forma, quizás más 'de bajo nivel', sería usar directamente fionay shapelypara E / S y procesamiento de geometría.

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)

Creo que la forma más fácil es usar gdal / ogr Virtual Format. ( http://www.gdal.org/drv_vrt.html ) y dialecto SQL / SQLITE ( http://www.gdal.org/ogr_sql.html y https://www.gaia-gis.it/spatialite-3.0 .0-BETA / spatialite-sql-3.0.0.html )

Mi archivo de forma poligonal se llama poly.shp. Luego creo este XML como archivo llamado vrt.vrt. Dentro de este archivo (vrt.vrt), aquí el contenido para convertir a puntos

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

En este momento, puede integrar este archivo en Qgis para validarlo. Por supuesto, la representación es más lenta que la fuente sin procesar porque cada característica se convierte como punto en cada consulta de representación.

Después, convierta este archivo (vrt.vrt) en otra cosa usando las utilidades gdal / ogr de un shell / script de python

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

Obtiene un archivo de forma de punto llamado point_from_vrt.shp.