Conversão Centroid do recurso de polígono para pontos usando Python

Gostaria de converter alguns arquivos shp baseados em polígono que possuem vários recursos de polígono em pontos para cada recurso que representariam essencialmente o centríodo de cada recurso de polígono. Eu sei que no mundo ArcGIS eu poderia usar a ferramenta Feature To Point, mas gostaria de manter isso em um script que possa ser executado em PCs que não possuam arcpy, então estou procurando uma alternativa de código aberto. Alguém está ciente de uma biblioteca que eu poderia usar para isso, juntamente com alguma orientação sobre como alavancá-lo para fazer isso?

Você pode executar um ogr2ogrcomando (por exemplo, a partir de um OSGeo4w Shell). Por exemplo, em um shapefile de países:

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

O novo shapefile countries_centroid.shpdeve ser semelhante à entrada, mas conter apenas um ponto por [Multi] Polígono.

O @PEL também mostra um bom exemplo com ST_PointOnSurface, que é simples de substituir neste comando.

Algo semelhante pode ser feito em Python, se necessário, mas pode levar algumas linhas de código a mais:

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

Simplesmente use Fiona ou GeoPandas (Python 2.7.xe 3.x)

Alguns polígonos

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

Transformação em pontos (centróides)

# 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

Outra maneira, talvez mais "de baixo nível", seria usar diretamente fionae shapelypara processamento de E / S e geometria.

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)

Eu acho que a maneira mais fácil é usar o formato virtual gdal / ogr. ( http://www.gdal.org/drv_vrt.html ) e dialeto SQL / SQLITE ( http://www.gdal.org/ogr_sql.html e https://www.gaia-gis.it/spatialite-3.0 .0-BETA / spatialite-sql-3.0.0.html )

Meu shapefile de polígono é nomeado poly.shp. Então eu crio esse XML como um arquivo chamado vrt.vrt. Dentro deste arquivo (vrt.vrt), aqui o conteúdo a ser convertido em pontos

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

No momento, você pode integrar esse arquivo ao Qgis para validar. Com certeza, a renderização é mais lenta que a fonte bruta, pois cada recurso é convertido como ponto em cada consulta de renderização.

Depois, converta esse arquivo (vrt.vrt) em outra coisa usando os utilitários gdal / ogr de um shell / script python

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

Você obtém um shapefile de ponto chamado point_from_vrt.shp.