I have a Polyline Feature class that I wish to maintain in 3D Editing
Hi All,
I'm looking for guidance or a miracle..
I have a Polyline feature class that is Z aware when first built but to adjust and add new features to it. It seems I can only create another Feature class with the updated Z aware Polyline features. and generally having to delete and replace data in the database to maintain the static dataset.
I have tried ArcScene to edit in 3d with very limited success, ie it takes about 10 minutes to digitize a line of 15 vertices 250 metres long (display performance is slow). And then often it will not populate the Z attribute of the vertices along the line, I have only been successful once. (I have tried different Elevation Sources from Raster to TIN and a few others TIF etc).
If I restart the same session and reload the ArcScene document the next day and add a new feature it again takes 10 minutes and has lost the z awareness,
The Polyline Layer knows the Elevations are derived from the 3D Surface but it will not populate the Z value.
I would actually prefer to edit the Z aware features in ArcMap as it is linked to another relational database so if I split a polyline in half it's ID also needs to be updated in the relational database with appropriate new records added.
The 3D line itself only really needs it's Start point and Finish point Z values to be known as I am trying to calculate Gradient of the line over the surface.
To me this sounds like basic 3D editing but I am surprised that I cannot do it easily or quickly. as all I am wanting is the Z value interpolated from the surface to be stored with the Polyline itself.
Why change Polyline features you ask? Well it is all about decision processes of determining Gradient and keeping the gradient within physical constraints such as 12% slope between end points. If you shift the polyline up or down the surface in XY, the Z value has to change also, so I need to reflect this in as near as possible to real time during the editing process, whilst maintaining the links to the relational database.
I am just at a loss as to how to achieve this easily. The modifications could encompass 50 -100 alterations in any one session, and one gradient change affects the next polyline gradient in the network of line features.
Thanks
Mat Boonen , Hikurangi . NZ
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.
