#make_onp.py
#Takes the modified raw text trail information generated in example 1 and
#reformates it into a new shape file

# 1 import modules
# import ogr, os, sys
import osgeo.ogr as ogr
import osgeo.osr as osr
import string
import sys

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

# create a new data source and layer
newDataSource = driver.CreateDataSource('onp_wgs84.shp')
if newDataSource is None:
   print 'Could not create file'
   sys.exit(1)
NewLayer = newDataSource.CreateLayer('onp_wgs84', geom_type=ogr.wkbLineString)
# create all fields 
fieldID=ogr.FieldDefn('id', ogr.OFTInteger)
fieldlow=ogr.FieldDefn('lowpoint', ogr.OFTInteger)
fieldhigh=ogr.FieldDefn('highpoint', ogr.OFTInteger)
fieldwood=ogr.FieldDefn('wood_page', ogr.OFTInteger)
fieldmolvar=ogr.FieldDefn('molvar', ogr.OFTInteger)
fieldarea=ogr.FieldDefn('area', ogr.OFTString);fieldarea.SetWidth(32)
fieldTrailName = ogr.FieldDefn('trial_name',ogr.OFTString);fieldTrailName.SetWidth(32)
fieldComment = ogr.FieldDefn('comment', ogr.OFTString);fieldTrailName.SetWidth(128)
fieldLength=ogr.FieldDefn('Length', ogr.OFTReal)
fieldEffective=ogr.FieldDefn('e_length', ogr.OFTReal)

# add the fields to the shapefile
NewLayer.CreateField(fieldID)
NewLayer.CreateField(fieldTrailName)
NewLayer.CreateField(fieldLength)
NewLayer.CreateField(fieldarea)
NewLayer.CreateField(fieldEffective)
NewLayer.CreateField(fieldlow)
NewLayer.CreateField(fieldhigh)
NewLayer.CreateField(fieldComment)
NewLayer.CreateField(fieldmolvar)
NewLayer.CreateField(fieldwood)

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

master_points= open('./ONP/finial.txt','r')
number_of_trails=string.atoi( master_points.readline())
# loop through the lines in the text file
for i in range(number_of_trails):
      (trail_name,trail_length,e1,e2,effective_length,ID_number,wood_id,molvar_id,area,comments)=string.split(master_points.readline(),',')
      trail_length=float(trail_length)
      e1=int(e1)
      e2=int(e2)
      ID_number=int(ID_number)
      wood_id=int(wood_id)
      molvar_id=int(molvar_id)
      effective_length=float(effective_length)
      comments=comments[:-1]
      number_of_location_points=string.atoi(master_points.readline() )   
      
      # create an empty geometry for nodes of the trail
      trail_nodes = ogr.Geometry(ogr.wkbLineString)
      for i in range(number_of_location_points):
            tokens=string.split(master_points.readline(),',')
            x=float(tokens[0])
            y=float(tokens[1])
            z=float(tokens[2])
            trail_nodes.AddPoint(x, y,z)

      #create a new feature and set its geometry and attribute
      feature = ogr.Feature(featureDefn)
      feature.SetGeometry(trail_nodes)
      feature.SetField('id', ID_number) 
      feature.SetField('trial_name', trail_name)
      feature.SetField('Length', trail_length)
      feature.SetField('e_length', effective_length)
      feature.SetField('lowpoint', e1)
      feature.SetField('highpoint', e2) 
      feature.SetField('comment', comments) 
      feature.SetField('wood_page', wood_id)
      feature.SetField('molvar',molvar_id)
      feature.SetField('area', area)

      #add the feature to the shapefile
      NewLayer.CreateFeature(feature)
      # destroy the geometries and feature
      trail_nodes.Destroy()#Do NOT destroy geometries that come from an EXISTING feature
      feature.Destroy() 


master_points.close()
newDataSource.Destroy() #close and flush data


# create the *.prj file
aSpatialRef = osr.SpatialReference()# create the new SpatialReference
aSpatialRef.ImportFromEPSG(32610)  #EPSG:32610 - WGS 84 / UTM zone 10N
aSpatialRef.MorphToESRI() 
file = open('onp_wgs84.prj', 'w')
file.write(aSpatialRef.ExportToWkt())
file.close()