Thursday 4 August 2016

arcgis desktop - Is there a tool in GIS to find the surface of a profile section in a DEM?


For example I have a DEM of river bed over which we made several profiles – eg sketch: http://prntscr.com/cqwile Now in cross section, it will look like this http://prntscr.com/cqwx59 .



Is there a way in GIS to find the surface of the green area by adding the DEM and a profile line from which it should be measured?


Or maybe set a Plane Height from which to start.



Answer



No such tool, I decided to make one.


In order to test it, I created 4 copies of the same polyline, and used DEM available and interpolate shape tool to convert them into 3D lines.


INPUT


SCRIPT:


import arcpy
from math import hypot
inLines = r"C:\FELIX_DATA\SCRARCH\line3d.shp"


g=arcpy.Geometry()
with arcpy.da.UpdateCursor(inLines,("Shape@","WL","Channels","SectArea")) as cursor:
for shp,wl,channels,sArea in cursor:
part=shp.getPart(0); array=arcpy.Array();L=0
for i,p in enumerate(part):
if i==0:p1=p
else:
p2=p
L+=hypot(p1.X-p2.X,p1.Y-p2.Y)

p1=p
array.add(arcpy.Point(L,p.Z))
section=arcpy.Polyline(array)
LeftPoint=arcpy.Point(0,wl)
RightPoint=arcpy.Point(L,wl)
horizontal=arcpy.Polyline(arcpy.Array([LeftPoint,RightPoint]))
G=arcpy.FeatureToPolygon_management([section,horizontal],g)
channels=0;sArea=0
for pgon in G:
extent=pgon.extent

if 0.5*(extent.YMax+extent.YMin)>=wl: continue
channels+=1;sArea+=pgon.area
arcpy.AddMessage(channels)
cursor.updateRow((shp,wl,channels,sArea))

OUTPUT:


OUTPUT


Note:


I used in memory geometries as inputs to mighty polygon constructing tool. I think this is very GIS solution to the task. Script will work on v.1.0 and later, not on 9.3.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...