Konversi Centroid Fitur Polygon ke Poin menggunakan Python

Saya ingin mengonversi beberapa file shp berbasis poligon yang memiliki beberapa fitur poligon menjadi poin untuk setiap fitur yang pada dasarnya akan mewakili centriod dari setiap fitur poligon. Saya tahu di dunia ArcGIS saya bisa menggunakan alat Feature To Point tetapi saya ingin menyimpannya dalam skrip yang dapat dijalankan di PC yang tidak memiliki arcpy di atasnya, jadi saya mencari alternatif open source untuk itu. Adakah yang tahu perpustakaan yang bisa saya gunakan untuk ini bersama dengan beberapa arahan tentang bagaimana memanfaatkan itu?

Anda dapat menjalankan ogr2ogrperintah (mis. Dari OSGeo4w Shell). Misalnya pada shapefile negara:

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

Shapefile baru countries_centroid.shpharus mirip dengan input, tetapi hanya berisi satu titik per [Multi] Polygon.

@PEL juga menunjukkan contoh yang baik dengan ST_PointOnSurface, yang mudah untuk diganti dalam perintah ini.

Hal serupa dapat dilakukan dengan Python, jika diperlukan, tetapi mungkin perlu beberapa baris kode lagi:

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

Cukup gunakan Fiona atau GeoPandas (Python 2.7.x dan 3.x)

Beberapa poligon

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

Transformasi ke titik (centroid)

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

Hasil

Cara lain, mungkin lebih 'rendah', adalah dengan langsung menggunakan fionadan shapelyuntuk I / O dan pemrosesan geometri.

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)

Saya pikir cara termudah adalah dengan menggunakan gdal / ogr Virtual Format. ( http://www.gdal.org/drv_vrt.html ) dan dialek SQL / SQLITE ( http://www.gdal.org/ogr_sql.html dan https://www.gaia-gis.it/spatialite-3.0 .0-BETA / spatialite-sql-3.0.0.html )

Shapefile poligon saya bernama poly.shp. Kemudian saya membuat file seperti XML ini bernama vrt.vrt. Di dalam file ini (vrt.vrt), di sini konten untuk dikonversi ke poin

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

Saat ini, Anda dapat mengintegrasikan file ini ke Qgis untuk divalidasi. Yang pasti, rendering lebih lambat dari sumber mentah karena setiap fitur ditampilkan sebagai titik pada setiap permintaan rendering.

Setelah itu, konversikan file ini (vrt.vrt) menjadi sesuatu yang lain menggunakan util gdal / ogr dari python shell / script

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

Anda mendapatkan shapefile titik bernama point_from_vrt.shp.