Continuing a bit on this cut fill part of the problem, I created a script based on the previous code (or actually based on this post Extract Raster Values using Numpy) that will create a point featureclass with the elevation of the projected road, the difference between elevation of the road and the DEM and an indication of Cut or Fill.
When you visualize the result it may look like this:
This is a 3D view of the road (black, draped over the DEM) and the points represent the projected road with the gradual elevation. The points on the road indicate cut (blue) or fill (red) and the amount (elevation difference between projected road and DEM).
The code used for this is:
import arcpy
import os
def main():
arcpy.env.overwriteOutput = True
# inputs
dem = r"D:\Xander\GeoNet\RoadSlope\gdb\test.gdb\DEM03"
fc = r"D:\Xander\GeoNet\AcresMultiSDE\roads.gdb\roads"
fc_out = r"D:\Xander\GeoNet\RoadSlope\gdb\test.gdb\points_v03"
fld_objid = "RoadObjectID"
fld_slope = "Slope"
fld_demZ = "ElevationDEM"
fld_roadZ = "ElevationRoad"
fld_difz = "ElevationDif"
fld_cutfill = "CutOrFill"
# describe input
sr = arcpy.Describe(fc).spatialReference
# create output fc
fc_ws, fc_name = os.path.split(fc_out)
arcpy.CreateFeatureclass_management(out_path=fc_ws, out_name=fc_name,
geometry_type="POINT", has_m="ENABLED",
has_z="ENABLED", spatial_reference=sr)
# add fields
AddField(fc_out, fld_objid, "LONG")
AddField(fc_out, fld_slope, "DOUBLE")
AddField(fc_out, fld_demZ, "DOUBLE")
AddField(fc_out, fld_roadZ, "DOUBLE")
AddField(fc_out, fld_difz, "DOUBLE")
AddField(fc_out, fld_cutfill, "TEXT", 12)
# get polyline and extract a part
cellsize = getRasterCellSize(dem)
# start insert cursor output points
flds_out = ("SHAPE@", fld_objid, fld_slope, fld_demZ, fld_roadZ, fld_difz, fld_cutfill)
with arcpy.da.InsertCursor(fc_out, flds_out) as curs_out:
flds = ("SHAPE@", "OID@")
with arcpy.da.SearchCursor(fc, flds) as curs:
for row in curs:
polyline = row[0]
oid = row[1]
pnts = [polyline.firstPoint, polyline.lastPoint]
# slope of projected road
elevs = []
for pnt in pnts:
# get extent from point and adapt extent to raster
ext1 = createExtentForPoint(pnt, cellsize)
ext2 = adaptExtentUsingRaster(ext1, dem, cellsize)
# create numpy array
np = createNumpyArrayUsingExtent(dem, ext2, cellsize)
r, c = getRowColExtentFromXY(ext2, pnt, cellsize)
z = np.item(r, c)
elevs.append(z)
zdif = elevs[1] - elevs[0]
slope = float(zdif / polyline.length)
# get extent of line
ext1 = getExtentFromPolyline(polyline)
ext2 = adaptExtentUsingRaster(ext1, dem, cellsize)
# create numpy array
np = createNumpyArrayUsingExtent(dem, ext2, cellsize)
# loop through points and extract values from raster
d = 0
pnts = getListOfPointsFromPolyline(polyline, cellsize)
l = polyline.length
for pnt in pnts:
r, c = getRowColExtentFromXY(ext2, pnt, cellsize)
try:
z = np.item(r, c)
frac = float(d) / float(l)
z_pnt = zdif * frac + elevs[0]
dif_road_dem = z_pnt - z
if dif_road_dem > 0:
cutfill = "Fill"
elif dif_road_dem < 0:
cutfill = "Cut"
else:
cutfill = "-"
except:
z = None
frac = float(d) / float(l)
z_pnt = zdif * frac + elevs[0]
dif_road_dem = None
cutfill = "Error"
# write to output
pnt2 = arcpy.Point(pnt.X, pnt.Y, z_pnt, d)
pnt_geo = arcpy.PointGeometry(pnt2, sr, True, True)
row_out = ((pnt_geo, oid, slope, z, z_pnt, dif_road_dem, cutfill, ))
curs_out.insertRow(row_out)
d += cellsize
def AddField(fc, fld_name, fld_type, fld_length=25):
if len(arcpy.ListFields(dataset=fc, wild_card=fld_name)) == 0:
if fld_type == "TEXT":
arcpy.AddField_management(in_table=fc, field_name=fld_name,
field_type=fld_type, field_length=fld_length)
else:
arcpy.AddField_management(fc, fld_name, fld_type)
def createExtentForPoint(pnt, cellsize):
return arcpy.Extent(pnt.X - cellsize, pnt.Y - cellsize, pnt.X + cellsize, pnt.Y + cellsize)
def getRowColExtentFromXY(ext, pnt, cellsize):
# r, c start at 0 from upper left
col = int(((pnt.X - ext.XMin) - ((pnt.X - ext.XMin) % cellsize)) / cellsize)
row = int(((ext.YMax - pnt.Y) - ((ext.YMax - pnt.Y) % cellsize)) / cellsize)
return row, col
def createNumpyArrayUsingExtent(raster, ext, cellsize):
lowerLeft = arcpy.Point(ext.XMin, ext.YMin)
ncols = int(ext.width / cellsize)
nrows = int(ext.height / cellsize)
return arcpy.RasterToNumPyArray(raster, lowerLeft, ncols, nrows)
def adaptExtentUsingRaster(ext, raster, cellsize):
ras_ext = arcpy.Describe(raster).extent
xmin = ext.XMin - ((ext.XMin - ras_ext.XMin) % cellsize)
ymin = ext.YMin - ((ext.YMin - ras_ext.YMin) % cellsize)
xmax = ext.XMax + ((ras_ext.XMax - ext.XMax) % cellsize)
ymax = ext.YMax + ((ras_ext.YMax - ext.YMax) % cellsize)
return arcpy.Extent(xmin, ymin, xmax, ymax)
def getRasterCellSize(raster):
desc = arcpy.Describe(raster)
return (desc.meanCellHeight + desc.meanCellWidth) / 2
def getListOfPointsFromPolyline(line, interval):
pnts = []
d = 0
max_len = line.length
while d < max_len:
pnt = line.positionAlongLine(d, False)
pnts.append(pnt.firstPoint)
d += interval
# pnts.append(line.lastPoint)
return pnts
def getExtentFromPolyline(line):
return line.extent
if __name__ == '__main__':
main()
Usage:
- line 8: reference to DEM
- line 9: reference to 2D roads
- line 10: output point featureclass (will be created)
Please note that code and the previous code will only work correct if the cellsize is round number.