The goal is to take a shape file that contain trail locations (i.e. line features) for the Olympic National Park and create the corresponding profile for each trail. The shape file is relatively straightforward in that it contains only one layer of line information. Specifically the file contains 872 features where each feature includes 7 attributes (of which the sixth (or 5 if you start counting at 0)field containing the trail name will turn out to be a critically important organizing parameter) and associated set of points (i.e. geometry) with define a small segment of a trail (in rare instances a feature contains all the x,y nodes for a given trail). These trail segments cannot be assumed to follow one another in the ‘table’ or be in any particular order (e.g. the middle trail segment may occur first in the table). Further, it is not guaranteed that the line segments have the same direction. In addition, there is no elevation data (i.e. all z values are set to zero) associated with the points/nodes that define a line. Thus, the first step in creating a profile is to group all the segments associated with a given trail together and add elevation to each x,y node. The following short program does this.
from osgeo import osr
from osgeo import ogr
import string
from math import sqrt
import raster as RS #A 'local' abstraction of GDAL optimized for hyperspectral Image processing
def dump(text,this_object):
print text+'******************************************'+text
help(this_object)
print text+'******************************************'+text
Driver = ogr.GetDriverByName('ESRI Shapefile')#Returns: a new Geo::OGR::Driver object.
dump("driver",Driver)#see dump
DataSource = Driver.Open('nad27_trails.shp')#doc
dump("DataSource",DataSource)#see dump
line_layer = DataSource.GetLayerByIndex(index=0)#Returns:a new Geo::OGR::Layer object
dump("line_layer ",line_layer)#see dump
#OGR::Layer Class Reference
number_of_features=line_layer.GetFeatureCount()
#a feature contains a segment of the line which defines part of a given trail.
#These segments,which in total define the location of a trail, do not necessarily:
# follow directly one after another in the shape file
# 'point' in the same direction OR
# are in order
dem=RS.open_raster_file('D:/hiker/basic/very_basic/5m_dem_park.tif')
Segments_of_a_trail={}
print "number of features=",number_of_features
count=1
for featureID in range(number_of_features):
#OGRFeature Class Reference
feature = line_layer.GetNextFeature()
if(count==1):dump("feature ",feature)#see dump
Geometry =feature.GetGeometryRef().Clone()#Fetch geometry from container
#Geo::OGR::Geometry Class Reference
if(count==1):dump("Geometry",Geometry)#see dump
number_of_points_this_feature=Geometry.GetPointCount()#doc
this_line_segment=[]
for k in range(number_of_points_this_feature):
(x,y,no_elevation_data_for_z)=Geometry.GetPoint(k)#doc
this_line_segment.append( ( x,y,dem.Get_pixel_value_by_GeoCoordinates(x,y)[0] ) )
trail_name=feature.GetFieldAsString(5)#see table
#The spelling of trail names are known to be consistant in this file and can be used to
#group corresponsing trail segments
if ( not (trail_name in Segments_of_a_trail.keys()) ):Segments_of_a_trail[trail_name]=[] #start a list
Segments_of_a_trail[trail_name].append(this_line_segment)
count+=1
feature.Destroy()
#now that I have all the points from the shape file write them out
raw_text= open('ONP_trails.txt','w')
total_number_of_trails=len(Segments_of_a_trail.keys())
raw_text.write(str(total_number_of_trails)+'\n')
for trail_name in Segments_of_a_trail.keys():
raw_text.write(trail_name+'\n')
number_of_segments=len(Segments_of_a_trail[trail_name])
raw_text.write(str( number_of_segments )+'\n')
for this_list_of_points in Segments_of_a_trail[trail_name]:
number_of_points=len(this_list_of_points)
raw_text.write(str(number_of_points)+'\n')
for a_point in this_list_of_points:
x=round(a_point[0],2)
y=round(a_point[1],2)
z=round(a_point[2],2)
raw_text.write(str(x)+','+str(y)+','+str(z)+'\n')
raw_text.close()