1

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

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.

3

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.

  • I got this error when i tried the change u suggested @webrian Traceback (most recent call last): line 20, in <module> f = inFeatureDefn(i) TypeError: 'FeatureDefn' object is not callable – user2841098 Oct 10 '13 at 0:51
  • It's f = inFeatureDefn.GetFieldDefn(i) – webrian Oct 10 '13 at 3:02
  • now my main concern is if I'm using the best way to solve this problem: Using GDAL-Python to read points in a polygon and make a new shape file (point) then find the min and max for X coordinates. @webrian – user2841098 Oct 10 '13 at 3:20
  • Just use the GetExtent method of a Layer, like I described in the answer. – webrian Oct 10 '13 at 4:29

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.