Python - GDAL : Polygon to point shapefile

Question

I want to read all the points in a polygon shapefile and use their coordinates to make a new point shapefile, then find the min and max of the X coordinates. So here is the code I wrote for the first part (polygon to point) and it makes the point shapefile but the attribute table is just has the ID field.

from osgeo import gdal
from osgeo import ogr
from osgeo import gdal_array
from osgeo import gdalconst
from osgeo.gdalconst import *
import os
os.chdir('C:\Fcounty')
driver = ogr.GetDriverByName("ESRI Shapefile")
fname = 'PolygonFC.shp'
vector = driver.Open(fname, 0)
layer = vector.GetLayer(0)
ofname = 'PointFC.shp'
if os.path.exists(ofname):
         driver.DeleteDataSource(ofname)

ovector = driver.CreateDataSource(ofname)
olayer = ovector.CreateLayer('points', geom_type=ogr.wkbPoint)
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
olayer.CreateField(fieldDefn)
featureDefn = olayer.GetLayerDefn()
point = ogr.Geometry(ogr.wkbPoint)
k = 0
f = layer.GetNextFeature()
while f:
    geom = f.GetGeometryRef()
    ring = geom.GetGeometryRef(0)
    np = ring.GetPointCount()
    for j in range(np):
        p = ring.GetPoint(j)
        ofeature = ogr.Feature(featureDefn)
        point.AddPoint(p[0], p[1])
        ofeature.SetGeometry(point)
        ofeature.SetField('id', k)
        olayer.CreateFeature(ofeature)
        ofeature.Destroy()
        k += 1
   f.Destroy()
   f = layer.GetNextFeature()
vector.Destroy()
ovector.Destroy()

As I mentioned this part is not giving me any error (I know it does not mean it works). So my questions are:

1. Is this the right code? why there is no other fields beside 'ID' in attribute table?
2. How can I do the second part of it? (the coordinates for lowest and highest X) I want to get it from the point shapefile, not the polygon.

Answer 1

As I understand your question, what you ultimately want is the minX and maxX of each polygon (i.e. its extent in the west-east direction). If so, I think you are over complicating things. You can get this simply by determining the extent of each polygon and there is no need to extract each individual vertex. This should also be considerably quicker than making a mass of csv files and then reading their points to determine the minX and maxX.

I will go with OGR as that is what you are using here (though in your other related post you are using ArcPy):

from osgeo import ogr
import os

os.chdir(r'C:\Fcounty')
ds = ogr.Open("PolygonFC.shp")
layer = ds.GetLayer(0)
f = layer.GetNextFeature()

arrayOfMinsAndMaxs = [] 

while f:
    geom = f.GetGeometryRef()
    bbox = geom.GetEnvelope()  # returns envelope as (minX, minY, maxX, maxY)

    # add the current poly FID, minX and maxX to a list of arrays
    minX = bbox[0]
    maxX = bbox[2]
    arrayOfMinsAndMaxs.append([f.GetFID, minX, maxX])

    f = layer.GetNextFeature()

In this code I iterate over the polygons, get the extent of each one and then save them to an array in the following format: [FID, minX, maxX]. So, you final array will look something like this:

[[0,15,500],[1,23,47],[2,76,128]....[500,2,10],[501,789,903]]

(obviously I've just made up those figures!). If you have some other id field you prefer then you could use that instead of the FID field. Also, I'm just dumping the values into an array, but you could save them to a text file or whatever. The array is just for demonstration purposes and with that in mind I spelled out how to get the minX and maxX.

Answer 2

You create and add only one field to the output layer in your code:

fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
olayer.CreateField(fieldDefn)

You have to do that for every field in the output layer, or you can copy the attributes from the input layer:

inFeatureDefn = layer.GetLayerDefn()
for i in range(inFeatureDefn.GetFieldCount()):
    f = inFeatureDefn.GetFieldDefn(i)
    olayer.CreateField(f)

To get the lowest and highest X coordinate from a layer you can use the GetExtent method in the Layer class, see also the Python API documentation.