#this example generated by Chris Garrard for a class on
# "Geoprocessing with Python using Open Source GIS"
#Create polygons from vertices in a text file and write the resulting geometries to a shapefile. The text
#file ut_counties.shp has one line for each county in Utah. You need to read each line in, create a
#polygon from the vertices contained in that line, and write the resulting polygon out to a shapefile. Each
#line looks like this:
#county_name:x1 y1,x2 y2,…,xn yn
#where xn = x1 and yn = y1
#1. Import needed modules and set the working directory
#2. Get the shapefile driver
#3. Create an output shapefile of type polygon and get its layer
#4. Add a field to the output shapefile for county name
## a. Create a FieldDefn for a string field
## b. Add a new field to your output shapefile using your FieldDefn from step 4a
#5. Now that you have added the field, get the FeatureDefn for your layer
#6. Open ut_counties.txt for reading
#7. For each line in ut_counties.txt
# a. Create an empty ring geometry
# b. Split the line on colons to get the county name and a string containing the coordinates
# c. Split the coordinate string from step 7b on commas to get a list of coordinates
# d. For each x,y pair in the list you got in step 7c
# i. Split the pair on spaces to get individual x and y values
# ii. Convert the x and y values to floating point values instead of strings
# iii. Use x and y to add a vertex to the ring created in step 7a
# e. Create a polygon geometry and add your ring to it
# f. Create a new feature using the FeatureDefn you got in step 5
# g. Set the geometry and county name for the feature
# h. Add the feature to your layer
# i. Destroy the two geometries that you created and the feature
#8. Close the text file and shapefile
# import modules
import ogr, os, sys
# set the working directory
os.chdir('d:/data/classes/python/data')
# get the shapefile driver
driver = ogr.GetDriverByName(
'ESRI Shapefile')
# create a new data source and layer
fn = 'hw2a.shp'
if os.path.exists(fn):driver.DeleteDataSource(fn)
ds = driver.CreateDataSource(fn)
if ds is None:
print 'Could not create file'
sys.exit(1)
layer = ds.CreateLayer('hw2a', geom_type=ogr.wkbPolygon)
# create a field for the county name
fieldDefn = ogr.FieldDefn('name', ogr.OFTString)
fieldDefn.SetWidth(30)
# add the field to the shapefile
layer.CreateField(fieldDefn)
# get the FeatureDefn for the shapefile
featureDefn = layer.GetLayerDefn()
# open the input text file for reading
file = open('ut_counties.txt', 'r')
# loop through the lines in the text file
for line in file:
# create an empty ring geometry
ring = ogr.Geometry(ogr.wkbLinearRing)
# split the line on colons to get county name and a string with coordinates
tmp = line.split(':')
name = tmp[0]
coords = tmp[1]
# split the coords on commas to get a list of x,y pairs
coordlist = coords.split(',')
# loop through the list of coordinates
for coord in coordlist:
# split the x,y pair on spaces to get the individual x and y values
xy = coord.split()
x = float(xy[0])
y = float(xy[1])
# add the vertex to the ring
ring.AddPoint(x, y)
# now that we've looped through all of the coordinates, create a polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
# create a new feature and set its geometry and attribute
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
feature.SetField('name', name)
# add the feature to the shapefile
layer.CreateFeature(feature)
# destroy the geometries and feature
ring.Destroy()
poly.Destroy()
feature.Destroy()
# now that we've looped through the entire text file, close the files
file.close()
ds.Destroy()