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

#Reproject the shapefile (unprojected geographic coordinates/NAD83 datum)
#created in extersoze 3 to the UTM projection, Zone 12N, but still in NAD83 by:
#1. Import needed modules and set the working directory
#2. Get the shapefile driver
#3. Create the CoordinateTransformation that is used to reproject each geometry
#            a. Create the input SpatialReference by importing information with EPSG 4269
#            b. Create the output SpatialReference by importing information with EPSG 26912
#            c. Use these two SpatialReferences to create a CoordinateTransformation
#4. Open your output shapefile from Example 3 and get the input layer
#5. Create an output shapefile of type polygon and get its layer
#6. Copy the county name field from the input layer to the output layer
#            a. Get a feature from the input layer
#            b. Get the FieldDefn for the county name field
#            c. Use that FieldDefn to create a new field in the output layer
#7. Now that you have added the field, get the FeatureDefn for your output layer
#8. For each feature in the input layer
#            a. Get the feature’s geometry
#            b. Use the CoordinateTransformation created in step 3 to reproject the geometry
#            c. Create a new feature using the FeatureDefn from step 7
#            d. Set the geometry for the new feature to the reprojected geometry
#            e. Copy the county name from the input feature to the new feature
#            f. Add the feature to the output layer
#            g. Destroy both features (not the geometry because it came from an existing feature)
#            h. Get the next input feature
#9. Close both shapefiles
#10. Create a .prj file for your output shapefile
#            a. Convert your SpatialReference from step 3b to ESRI format
#            b. Open a text file for writing, making sure you give it the correct name
#            c. Write out the WKT for the SpatialReference to the text file
#            d. Close the text file

# import modules
import ogr, osr, os, sys

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

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

# create the input SpatialReference
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(4269)

# create the output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(26912)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

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

# create a new data source and layer
fn = 'hw2b.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('hw2b', geom_type=ogr.wkbPolygon)

# get the FieldDefn for the county name field
feature = inLayer.GetFeature(0)
fieldDefn = feature.GetFieldDefnRef('name')

# add the field to the output shapefile
outLayer.CreateField(fieldDefn)

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

# loop through the input features
inFeature = inLayer.GetNextFeature()
while inFeature:
             # get the input geometry
             geom = inFeature.GetGeometryRef()

             # reproject the geometry
             geom.Transform(coordTrans)

             # create a new feature
             outFeature = ogr.Feature(featureDefn)

            # set the geometry and attribute
             outFeature.SetGeometry(geom)
            outFeature.SetField('name',             inFeature.GetField('name'))

             # add the feature to the shapefile
             outLayer.CreateFeature(outFeature)

             # destroy the features and get the next             input feature
             outFeature.Destroy
             inFeature.Destroy
             inFeature = inLayer.GetNextFeature()

# close the shapefiles
inDS.Destroy()
outDS.Destroy()#and flush all data>

# create the *.prj file
outSpatialRef.MorphToESRI()
file = open('hw2b.prj', 'w')
file.write(outSpatialRef.ExportToWkt())
file.close()