#this example generated by Chris Garrard for a class on
# "Geoprocessing with Python using Open Source GIS"


#It illustrates creating a new shapefile file and copying selected features from the old shapefile to the new one by:
#1. Import needed modules and set the working directory
#2. Get the shapefile driver
#3. Open sites.shp and get the input layer
#4. Create an output shapefile and get its layer
#5. Copy the ‘id’ and ‘cover’ fields from the input layer to the output layer
#            a. Get a feature from the input layer
#            b. Get the FieldDefn’s for the ‘id’ and ‘cover’ fields
#            c. Use that those FieldDefn’s to create two new fields in the output layer
#6. Get the FeatureDefn for the output layer
#7. Loop through the features in the input layer
#            a. Get the ‘cover’ attribute for the current input feature
#            b. If cover is ‘trees’ (remember the difference between = and ==)
#                        i. Create a new output feature using the FeatureDefn you got in step 6
#                        ii. Get the geometry for the input feature and add it to the output feature
#                        iii. Add the cover attribute of the input feature to the output feature
#                        iv. Get the id attribute for the input feature and add it to the output feature
#                        v. Add the output feature to the output layer
#                        vi. Destroy the output feature
#            c. Destroy the input feature and get the next input feature
#8. Destroy both of the DataSources

# 1 import modules
import ogr, os, sys

# set the working directory
os.chdir('d:/data/classes/python/data')

# 2 get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# 3 open the input data source and get the layer
inDS = driver.Open('sites.shp', 0) #shows cover at given points
if inDS is None:
print 'Could not open file'
sys.exit(1)
inLayer = inDS.GetLayer()

# 4 create a new data source and layer
fn = 'hw1b.shp'
if os.path.exists(fn):driver.DeleteDataSource(fn)
outDS = driver.CreateDataSource(fn)
if outDS is None:
   print 'Could not create file'
   sys.exit(1)
outLayer = outDS.CreateLayer('hw1b', geom_type=ogr.wkbPoint)

# 5 get the FieldDefn's for the id and cover fields in the input shapefile
feature = inLayer.GetFeature(0)
idFieldDefn = feature.GetFieldDefnRef('id')
coverFieldDefn = feature.GetFieldDefnRef('cover')

#create new id and cover fields in the output shapefile
outLayer.CreateField(idFieldDefn)
outLayer.CreateField(coverFieldDefn)

# 6 get the FeatureDefn for the output layer
featureDefn = outLayer.GetLayerDefn()

# 7 loop through the input features
inFeature = inLayer.GetNextFeature()
while inFeature:

         # get the cover attribute for the input feature
         cover = inFeature.GetField('cover')

         # check to see if cover == trees
         if cover == 'trees':

                # create a new feature
                outFeature = ogr.Feature(featureDefn)#using featureDefn created in step 6

                # set the geometry
                geom = inFeature.GetGeometryRef()
                outFeature.SetGeometry(geom) #move it to the new feature

                # set the attributes
                id = inFeature.GetField('id')
                outFeature.SetField('id', id) #move it to the new feature
                outFeature.SetField('cover', cover) #move it to the new feature

                # add the feature to the output layer
                outLayer.CreateFeature(outFeature)

               # destroy the output feature
               outFeature.Destroy()

         # destroy the input feature and get a new one
         inFeature.Destroy()
         inFeature = inLayer.GetNextFeature()

# close the data sources
inDS.Destroy()
outDS.Destroy()#flush out the last changes here