news: esri is coming out with the new version x 10.1 and the included source net and javahere you can see my preview beta http://sit.sistemigis.it/Samples/Elevations/however you download esri version ...
SETUP1) Create a Map Document (mxd) with at least one elevation layer; the elevation layer must be a single band raster layer. Other types of layers are ok to use because the SOE will only expose a resource for the single band raster layers.2) Publish the Map Document as a Map Service and then stop the map service after it has been successfully published.3) Right-click on the map service and select �??Service Properties�??. In the �??Map Service Properties�?? dialog, select the �??Capabilities�?? tab. 4) From the list of capabilities, check on the new �??Get Elevations�?? capability and click the OK button.5) Start the map service.6) Go to the Services Directory admin page and clear the REST cache.7) Navigate to the Services Directory page of the Map Service and you�??ll now see the Elevations SOE listed at the bottom of the page in the �??SUPPORTED EXTENSIONS�?? section. Follow the available links to view the available REST resources and operations for the Elevations SOE.
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()
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)
Signed in members can post, follow updates, and more. New here? Register a free account.
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.