Python - GDAL : Polygon to point shapefile

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

As I understand your question, what you ultimately want is the min-X and max-X 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 min-X and max-X.

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 impost ogr
import os


os.chdir(r'C:\Fcounty')
ds = ogr.Open("PolygonFC.shp")
layer = ds.GetLayer(0)
f = lyr.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 = lyr.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.