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.
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:
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