I would like to convert some polygon based shp files that have several polygon features into points for each feature that would essentially represent the centriod of each polygon feature. I know in ArcGIS world I could use Feature To Point tool but I would like to keep this in a script that could be run on PC's that don't have arcpy on them so I'm looking for an open source alternative to that. Is anyone aware of a library I could use for this along with some direction on how to leverage it accomplish this?
You can run an
ogr2ogr command (e.g. from a OSGeo4w Shell). E.g. on a shapefile of countries:
cd path/to/shapefiles ogr2ogr -sql "SELECT ST_Centroid(geometry), * FROM countries" -dialect sqlite countries_centroid.shp countries.shp
The new shapefile
countries_centroid.shp should be similar to the input, but just contain one point per [Multi]Polygon.
@PEL also shows a good example with ST_PointOnSurface, which is simple to substitute in this command.
Something similar can be done in Python, if needed, but it may take a few lines of code more:
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
import geopandas as gpd # GeoDataFrame creation poly = gpd.read_file("geoch_poly.shp") poly.head()
Transformation to points (centroids)
# 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')
Another, perhaps more 'low level', way would be to directly use
shapely for I/O and geometry processing.
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)
I think the easiest way is to use gdal/ogr Virtual Format. (http://www.gdal.org/drv_vrt.html) and SQL /SQLITE dialect (http://www.gdal.org/ogr_sql.html and https://www.gaia-gis.it/spatialite-3.0.0-BETA/spatialite-sql-3.0.0.html)
My polygon shapefile is named poly.shp. Then I create this XML like file named vrt.vrt . Inside this file (vrt.vrt), here the content to convert to 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>
At this time, you can integrate this file into Qgis to validate. For sure, the rendering is slower than the raw source because every feature is cast as point on every rendering query.
After, convert this file (vrt.vrt) into something else using gdal/ogr utils from a python shell/script
os.system("ogr2ogr point_from_vrt.shp vrt.vrt poly")
You get a point shapefile named point_from_vrt.shp .