8

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?

  • I'm still having multiple issues with the answer that Gene provided below. The issues are how it's reorders the attributes from their original order to alphabetical which is a problem. Secondly, the shape file get's corrupt possibly due to the file I'm trying to convert having over 250 attributes. – wilbev Jan 9 '17 at 5:14
  • There is a standard tool called 'Polygon Centroids' in QGIS that does exactly this - do you require a script? It would be easy enough to script using PyQGIS I would think. – Dùn Caan Jan 9 '17 at 13:52
  • It needs to be a script and work on PC's that do not have QGIS on them. – wilbev Jan 10 '17 at 5:00
8
+50

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
  • I think you've proposed and simpliest solution with OGR and SQL. I think it's safer to add a parameter to OGR with -nlt Point – PEL Jan 11 '17 at 16:18
  • Unfortunately I can't get this to work. I get an error claiming that ST_Centroid is not function. – wilbev Jan 14 '17 at 1:09
  • 1
    It needs the SQLite dialect option (as shown) and Spatialite built into GDAL, which is not always guaranteed. OSGeo4W has a good version of GDAL that will run this command correctly. – Mike T Jan 14 '17 at 21:03
  • I can get your top script within the ogr2ogr commend to work without issue. However, I need to do this in a stand-alone python script so I'm trying to get your 2nd set of code to work which is where I continue to the ST_Centroid is not a function error. My code is identical to what you have above including the sqlite dialect. – wilbev Jan 15 '17 at 0:04
  • 1
    The error you describe is when GDAL was built without Spatialite support. Some gdal-python packages support this, but not all. Try opening an OSGeo4W shell and running the Python script from that environment. I think the default set of packages related to gdal-bin include this support. – Mike T Jan 15 '17 at 10:20
9

Simply use Fiona or GeoPandas (Python 2.7.x and 3.x)

Some polygons

enter image description here

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

enter image description here

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

enter image description here

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

Result

enter image description here

  • Thanks for the answer gene. I think you may have one typo above where the 'gdf' variable should be 'poly" right? in the points.crs = gdf.crs code. I'm also experiencing a couple of other issues where the .prj file is not getting created, it's showing up empty and the order of the attribute fields are changing their order from the polygon data since they seem to now be alphabetical. It's important they stay in the same order. You know of a way to keep the field attribute in same order? – wilbev Nov 6 '16 at 23:24
  • Thanks, corrected. For the order of the attributes fields, simply change the order of the GeoPandas GeoDataFrame (= Pandas DataFrame) – gene Nov 7 '16 at 7:19
  • Thanks Gene, but I'm not sure I understand where I would change that order in the code here. I'm also running into two other issues with this. First the *.prj file is blank on the new shp file. Secondly, when I attempt to open the shp file in shp reader it gives an error opening the file like it is corrupt. It seems to work without getting corrupt if the shp file only has a single feature but multiples is where it's having trouble it seems. – wilbev Nov 7 '16 at 14:26
  • Sorry, but you need to know Pandas for that and I have no problem with the script with my data (I use the latest versions of GeoPandas, Fiona and Numpy) – gene Nov 7 '16 at 14:31
  • I can send you the shp file so you can see for yourself but this shp file has over 250 columns of data that I imagine are creating problems for it. I tried this on a shp file with much fewer attributes and it doesn't seem to have any trouble. – wilbev Nov 8 '16 at 1:39
5

Another, perhaps more 'low level', way would be to directly use fiona and 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)
  • Thanks Loic. That definitely fixes the sort issue I was having but it's not fixing the issue with so many attributes causing the file to become corrupt. Would you have any other ideas on how to overcome that issue? I guess I may need to delete attributes. I can send you an example file if it would help. – wilbev Jan 10 '17 at 1:43
  • @wilbev Send a download link to your data if you can. Otherwise I don't see what could be wrong. – Loïc Dutrieux Jan 10 '17 at 10:51
  • Loic, I emailed you an example file. Hopefully that gives you a good idea of the issue I'm running into. – wilbev Jan 11 '17 at 5:45
  • @wilbev what do you mean by 'the file becomes corrupt'? Using the file you sent, I could produce the centroids and open the output shapefile in QGIS without problems. The attributes table remains unchanged between both files. – Loïc Dutrieux Jan 11 '17 at 14:06
  • By corrupt, I mean it's essentially an empty dbf file since after i run the script it creates a 1 KB sized dbf file and when you open that it's completely empty. If I run the same exact script, the one you list above, on a file with fewer attributes, it works without issue. I even tried on a second PC and got the same result. I don't get it. – wilbev Jan 12 '17 at 4:43
2

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 .

  • I was able to get this to work but I would like to keep all of this within a python script since I need to convert 100's of files all with different file names. I would like to use @Mike T solution but I get a "No such function: ST_Centroid if I use that and I also tried ST_PointOnSurface which also says there is no such function. Any idea why it would be stating these are no functions of ExecuteSQL()? – wilbev Jan 14 '17 at 18:57
  • I get 'wkbPoints' is not a valid value of the atomic type – highmaintenance Oct 5 '17 at 13:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.