#Merge lake features from several diverse (e.g. different projections and area coverage) 
#GIS files  into one consistent shape file 
import os, sys 
import osgeo.ogr as ogr
import osgeo.osr as osr

old_file=[ ('C:\\_hiker_vectors\\all_alpine_lakes\\water\\raw_wen\\waterbody','wen'),
          ('C:\\_hiker_vectors\\all_alpine_lakes\\water\\raw_baker\\waterbody','baker')]

def add_to_shape_file(old_file_name, outSpatialRef,featureDefn,outLayer,this_area='washington',field_ID='name'):
   inDS = ogr.Open(old_file_name, 0)
   if inDS is None:
      print 'Could not open',old_file_name 
      sys.exit(1)
   old_layer = inDS.GetLayer('PAL') #all water features of interest are polygons 
   OldSpatialRef=old_layer.GetSpatialRef()
   coordTrans = osr.CoordinateTransformation(OldSpatialRef, outSpatialRef)# create the CoordinateTransformation
   inFeature = old_layer.GetNextFeature()
   while inFeature:
      outFeature = ogr.Feature(featureDefn)#create a new feature 
      geom = inFeature.GetGeometryRef().Clone() #get lake polygon
      geom.Transform(coordTrans)# reproject the geometry
      outFeature.SetGeometry(geom) #move geometry to the new feature
      outFeature.SetField('area', this_area)
      outFeature.SetField('name',inFeature.GetField(field_ID))#name of the lake
      outLayer.CreateFeature(outFeature) #add new attributes(area & name) and lake polygon to the new shape file
      outFeature.Destroy()# destroy the output feature
      inFeature.Destroy()# destroy the input feature  
      inFeature = old_layer.GetNextFeature() #get the next feature
   # close the old data source
   inDS.Destroy()


#set up the new file !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!**********************
driver = ogr.GetDriverByName('ESRI Shapefile')#get the shapefile driver
outDS = driver.CreateDataSource('alllakes.shp') 
if outDS is None:
   print 'Could not create file'
   sys.exit(1)
outLayer = outDS.CreateLayer('lakes', geom_type=ogr.wkbPolygon)
fieldarea = ogr.FieldDefn('area',ogr.OFTString);fieldarea.SetWidth(16)
fieldname = ogr.FieldDefn('name',ogr.OFTString);fieldname.SetWidth(32)
outLayer.CreateField(fieldname)
outLayer.CreateField(fieldarea)
featureDefn = outLayer.GetLayerDefn() #get the FeatureDefn for the output layer
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!************************

outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(26710) #EPSG Projection 26710(NAD27 / UTM zone 10N) 


for a_file,area in old_file:
   print "adding file: ",a_file,' for area =',area
   add_to_shape_file(a_file,outSpatialRef,featureDefn,outLayer,area)



outDS.Destroy()#flush the last changes here
#add projection information here
outSpatialRef.MorphToESRI() 
file = open('alllakes.prj', 'w')
file.write(outSpatialRef.ExportToWkt())
file.close()