I managed to hack this together at one point, but found it was too slow for a web service. I rewrote the script using gdal for raster access and it was hundreds of times faster. The code for the arcpy script is below, with an example output graph. The commented code was for grabbing points from a feature set (for a toolbox or gp service). This ugly, but it works. disclaimer: This is hacked, old and ugly. sorry.
from polyline import *
import matplotlib.pyplot as plt
def fetch_cell_value(x, y):
xy_string = "{0} {1}".format(x, y)
band = "1"
arcpy.GetCellValue_management(r, xy_string, band)
results = arcpy.GetMessages()
lines = results.split('\n')
value = float(lines[2])
return value
#for arc script
#fs = arcpy.GetParameter(0)
#rows = arcpy.SearchCursor(fs)
# Assume feature set has only one feature and one shape.
#row = rows.next()
#shape = row.getValue("SHAPE")
#arcpoints = shape.getPart(0)
#points = list()
#for p in arcpoints:
# points.append([p.X, p.Y])
raster = "c:/data/someraster.tif"
r = arcpy.Raster(raster)
#test points
points = [[290626.474, 4875792.8594],
[282358.29,4863026.666],
[297836.3843, 4870699.62]]
p = poly_line(points)
n_points = 100
#if p.length() / n_points < cell_size:
# skip = cell_size
#else:
# skip = int(p.length() / n_points)
skip = p.length() / n_points
dists = list()
i = 0
while i * skip < p.length():
dists.append(i * skip)
i += 1
values = list()
for d in dists:
x,y = p.d_line(d)
v = fetch_cell_value(x,y)
values.append(v)
print("Time:{0}".format(time.clock() - start))
plt.plot(dists, values)
plt.xlabel("Distance (m)")
plt.ylabel("Elevation (m)")
plt.title("Elevation Profile")
plt.grid(True)
plt.show()
The poly_line is a little class I wrote to handle densifying lines. It assumes euclidean geometry (planar) so doesn't handle earth curvature. You could probably use some type of densify function for an arc type polyline and use the points in the shape. I didn't try it. poly_line:
import math
##
# Class representing a segment of a polyline
#
class segment:
p1 = [0,0]
p2 = [0,0]
##
# @brief Construct a segment
#
def __init__(self, p1, p2):
self.p1 = p1
self.p2 = p2
##
# @brief Get the length of the segment
#
# @todo Test for /0
#
def length(self):
a = self.p2[0] - self.p1[0]
b = self.p2[1] - self.p1[1]
l = math.hypot(a,b)
return l
##
# @brief Get the slope of the segment
#
# @todo Test for /0
#
def m(self):
a = self.p2[1] - self.p1[1]
b = self.p2[0] - self.p1[0]
m = a / b
return m
##
# @brief Get the y intercept of the segment
#
# @return y intercept
def b(self):
b = self.p2[1] - (self.m() * self.p2[0])
return b
##
# @brief Get the coordinate of a point some distance from point 1
#
# @param d distance from point 1
# @return x and y coordinate
def d_line(self, d):
if d == 0:
return self.p1[0], self.p1[1]
if d > self.length():
return None
m = self.m()
u = d / math.sqrt(m * m + 1)
if self.p2[0] < self.p1[0]:
u = -u
v = (m * d) / math.sqrt(m * m + 1)
if self.p2[1] < self.p1[1]:
v = -v
return (self.p1[0] + u, self.p1[1] + v)
##
# Class representing a polyline composed of segments
#
# @see segment
#
class poly_line:
segments = list()
def __init__(self, points):
for i in range(len(points) - 1):
s = segment(points, points[i+1])
self.segments.append(s)
##
# @brief Get the length of all the segments
#
# @return total length
#
def length(self):
l = 0.0
for s in self.segments:
l += s.length()
return l
##
# @brief Get the coordinate of a point some distance from the polyline
# origin
#
# @param d distance from the origin
# @return x and y coordinate for the point
#
def d_line(self, d):
if d == 0:
return self.segments[0].d_line(0)
l = 0.0
for s in self.segments:
if l >= d:
break
dd = d - l
l += s.length()
seg = s
return seg.d_line(dd)
I just ran it and wrote the attached output, but no guarentees it will work for anyone. I think that the service ESRI wrote used lower level raster access, because of it's speed. I got comparable speed with the gdal version. Sorry if I didn't help, but it is possible. I stopped working on it, but it may help others. I will watch the thread and try to answer questions.kss
