Total Elevation Gain confusion

[JasonBarnes]

Greetings,

I'm having some confusion with results calculating total elevation gain/accumulation from a polylineZ file. Here's the situation, I received a gps track from a friend of a local bicycle training route. According to his gps track the total elevation gain was about 6,200'. I digitized a track based on the roads the gps track followed, and interpolated it to a 3d shape using 3d analyst and a 10m DEM. I wrote a python script to calculate the total elevation gain using the z values of this new polylineZ. The confusion is that my calculation from they python script is giving me approx 11,460', which is almost double what his gps track said. To test it more, I exported the feature to ascii, brought it into excel and wrote a VBA code which gave me the same exact calculation as the python script.

So, I decided to try something else. I then exported the file to a kml layer, and brought it into GPS visualizer (gpsvisualizer.com), and then looked at the total gain (vertical up) on EveryTrail.com. Their calculation for the KML file was really close to what his initial gps track stated, around 6,300'.

Perhaps something is wrong with my logic in calculating the total elevation gain:

Get z value, if it is larger than previous z value then add the difference to a total gain variable. Do until all z values have been cycled through.

Here is the code:

import arcpy

# Get the input feature class
infc = arcpy.GetParameterAsText(0)

# Name the output location where Total Accumulation result will go
opfl = arcpy.GetParameterAsText(1)

# Name of the file written to the above file location
filename = arcpy.GetParameterAsText(2)
outputfile = open(opfl + "/" + filename +".txt", "a")

# Identify the geometry field
#
desc = arcpy.Describe(infc)
shapefieldname = desc.ShapeFieldName

# Create search cursor
#
rows = arcpy.SearchCursor(infc)

# Enter for loop for each feature/row
#
for row in rows:
    # Create the geometry object
    # in this case Polyline Object
    feat = row.getValue(shapefieldname)

    # Print the current multipoint's ID
    #
    print "Feature %i:" % row.getValue(desc.OIDFieldName)
    partnum = 0

    # Step through each part of the feature
    #
    for part in feat:
        # Print the part number
        #
        print "Part %i:" % partnum

        # Step through each vertex in the feature
        prev_z = 10000
        tot_accum = 0
        for pnt in feat.getPart(partnum):
           
            if pnt.Z > prev_z:
                tot_accum = ((pnt.Z - prev_z) + tot_accum)
                prev_z = pnt.Z
            else:
                tot_accum = tot_accum
                prev_z = pnt.Z
        partnum += 1   

    print "Total Accumulation = %i:" % tot_accum

# Write the Total Accumulation to file and close file
outputfile.write("File Name: " + desc.Name + "\n")
outputfile.write("Total Accumulation = " + repr(tot_accum) + "\n \n")
outputfile.close()

# Add messages to tool output window
arcpy.AddMessage("Total accumulation tool is complete")
arcpy.AddMessage("File Name: " + desc.Name)
arcpy.AddMessage("Total Accumulation = " + repr(tot_accum)

Thanks for any help!

[MathewCoyle]

Is there a particular reason you set the prev_z value to 10000 instead of the z value of the first geometry part?

[JasonBarnes]

Is there a particular reason you set the prev_z value to 10000 instead of the z value of the first geometry part?

Matthew,

Since the If statement needs a value to compare the first z value to I put it at something high so it wouldn't be larger. If I set prev_z at 0 then the first vertex z value would be larger than zero and incorrectly add the difference of the two to the total accumulation. This was the best way I could think of to do it. I'm sure there are other ways.

Jason

[MathewCoyle]

That works well enough I would say. I just wanted to make sure it was high enough that it wouldn't be throwing a wrench in things on the first iteration. I assume since you keep overwriting your tot_accum under your cursor that there is only one feature in the feature class you are running this script on? I can't find any errors in your logic really. If you post the track you are processing it on that might help track down the issue.

[JasonBarnes]

That works well enough I would say. I just wanted to make sure it was high enough that it wouldn't be throwing a wrench in things on the first iteration. I assume since you keep overwriting your tot_accum under your cursor that there is only one feature in the feature class you are running this script on? I can't find any errors in your logic really. If you post the track you are processing it on that might help track down the issue.

Definitely no points above 10000 meters in this shapefile, and yes only one feature in the shapefile. I have a bunch of files I need to run total accumulation on, so I'm wondering if it will happen with all of them. I should perform the same process and see what the gain comparison is between my tool and the website. Glad to hear my logic seems ok, but perhaps something will be found. I did a little playing around with the script to print out the result if the difference between the two z values was greater then 5, just to see if by chance there was wicked vertex in the file. Only 3 values were above 5, but they were just over, so it looks like there's no crazy vertex causing a large output.

I've attached the polylineZ.

Thanks.

[MathewCoyle]

Your tool is calculating the Z values correctly. There may have been a problem with the DEM you used to create the Zpolyline or with the DEM the reference systems you are using calculate on. The only other thing I can imagine being wrong is some meter/feet conversion problem, but I can't find any evidence of that.

This tool calculates it the exact same.
http://resources.arcgis.com/gallery/file/Geoprocessing-Model-and-Script-Tool-Gallery/details?entryID...

[JasonBarnes]

Your tool is calculating the Z values correctly. There may have been a problem with the DEM you used to create the Zpolyline or with the DEM the reference systems you are using calculate on. The only other thing I can imagine being wrong is some meter/feet conversion problem, but I can't find any evidence of that.

This tool calculates it the exact same.
http://resources.arcgis.com/gallery/file/Geoprocessing-Model-and-Script-Tool-Gallery/details?entryID...

Thanks again for all your help Matthew. I do believe the tool is calculating the z values correctly, and to test it I created a short polyline with about 10 vertices and hand calculated the total gain and compared it to the tool's result. Got the exact same thing. Perhaps something is funky with the DEM. I'm using the USGS 10m NED. Is there any other options, and if so anything you would suggest? There shouldn't be a feet/meters conversion issue. The DEM is in meters so of course the z values for the polyline are in meters. I just convert the results to feet as a comparison to the original gps line and the website total gain results.

Thanks for the link to the tool. I had searched around quite a bit to find one and couldn't, how I missed this I'll never know. Either way the programming practice was good. I'll continue to try and figure this out.

[JasonBarnes]

Not that surprising, but after downloading a new set of DEMs from a different source, mosaic them and running the same analysis on my polylineZ I get different results:

Old DEM = 3493.5 meters (11,461')
New DEM = 3110.5 meters (10,205')

I'm a bit closer, however still about 4 thousand feet over what the gps track log and Everytrail.com website posted.

So, either the gps track and Everytrail.com are way off and what I'm getting is much closer, or the DEMs that are available are so inaccurate that they generate very poor results. At this point I'm unable to confidently decide on either one.