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
Solved! Go to Solution.
The polyline length in the script is 2D, not 3D. So you can determine the 3D length (assuming gradual increase/decrease of elevation) using the Pythagoras theorem.
@Xander Bakker
Yes because I only really need to know the 3Dlength of height variation, the assumption is that if the line is going from it's start point to it's finish point with a gradual continuous slope along the same line xy vertices to vertices ie as though the surface beneath the road is a continuous angled slope then I would need to apportion the variation in Z from vertices to vertices proportionately based on their xy location. This may be where the other Add surface information calculate functions come back into play, as mentioned in one of your first posts. Basically I am making a new surface from the height of the start point to the height of the finish point and then working out the 3D length of the same line on top of that new surface..
IE cutting it into a hill or filling a valley as appropriate.
If see where I am getting to? I will explain, I will want to eventually compare one surface to the next and determine these cuts and fills as surface volume of difference.
Yep, I understand. To change a copy of your DEM to the new road might be a little more complex. With some small changes to the code, you can create a point featureclass with point using a defined interval on the line and extract the interpolated elevation of the road the the current DEM elevation and symbolize the difference (cut/fill).
Would be interesting...
Hi Xander
Both to symbolize and to quantify cuts and fills into volume in cubic meters would be the ultimate goal. Again sort of goal, because potentially I could then start optimizing the volume variation to try to minimize the volume of soil moved.
As this is a green fields plan that seems to spawn from one need to the next. People are starting to realize the capability of GIS in the planning and civil / survey engineering sense. Allowing early decision models to be built on what if scenarios. In this case during the planning phase of new Harvest roads, and where they are being positioned for the least disturbance to the environment and possibly/probably the maximum return on capital spent. Doing so in a package that is not a CAD or civil road engineering design package seems logical to me too. As it will at minimum fill the gap between engineering design of every forest road section and peg layout in the field without much more than dead reckoning.
Anyway back on Subject I welcome all input into this topic as it can and does spawn ideas that can be taken further in a wide range of 3D modelling scenarios.
Cheers
Mat
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:
Please note that code and the previous code will only work correct if the cellsize is round number.
Very clever, now to give the same points width ie make the road 7 metres total width and to quantify the volume cut or filled from the 3d surface and the new 3d Road surface, I suppose I could get a property of the road direction or buffer these points by a 3.5 metre wide buy half the distance to the next or previous point and make that a rectangle and also place this rectangle on the same gradient. i.e. two from corners are half the height difference to the previous entity, and the two to corners are half the height difference to the next two corners (this may not work too well in a bend sense, maybe I have to buffer with the round function so as each polygon at least sort of meets up.
Realistically I was going to generate a new surface at a set width from the road centreline and be able to use this surface as a slice line to quantify volume above or below the original surface.
I know I dont want much, if too hard I understand, as you lost me in how you got to where we are a while back, all I know is it does work to visualise the result So I thank you for that.
Now I want to quantify it in terms of volume differences between the surfaces. and in so doing built a "pseudo road engineering package" out of GIS. Which is so cool I think Thanks Xander Bakker
Hi Xander Bakker I now have a very complex question to put a spanner in the works of the original Python Script. to apply a 3D value to a polyline. (road)
It appears that the user needs to determine the elevation of landings that cross over or are along parts of the road, the trouble with these landings is that they are a level flat surface at a as yet unknown elevation, but somewhere near the elevation of a neutral cut fill balance of the soil needed to make them flat.
so what this does is determines that the start and finish point on each polyline is at a different elevation than the underlying surface. I have toyed with the idea to make a new surface superimposing these landing flat areas into that surface, that in its own right has an issue because it is very user interactive to determine each elevation of each landing. These landings because they are flat also alter the gradient of the road either side of them as each road section is now shorter .
What I considered doing was to give each landing a static elevation (they are polygons) and superimpose them as a hard breakline in the surface. is this the most logical approach? is there a way to interactively use these landings as a height source without making a surface first.
The other part to this is not every road has these landings on them so I cant just say use them all.
as some (many) roads still need elevations from the underlying surface. being the DEM.
Any guidance would be appreciated, thank you. Mathieu
Hi Hikurangi Forest Farms Ltd ,
Maybe it is better to branch this part to a new discussion and include a link to this discussion. This will attract more people to respond.
Concerning your question, the neutral cut full locations can be determined with the points script that I posted before. It is not completely clear how these landings are determined and the effect they have on the road design. Is there some way you can visualize it or attach a small part of (dummy) data to explain it. Then I can have a look and see if this could and should be included in the original solution or it's better to design something new...
Kind regards, Xander
If you are trying to add Z values to the vertices between the beginning and ending point, then I don't have any suggestions. However, if you have already assigned Z values to your lines based on your surface and want to calculate M values that will provide the 3D length of your line and allow you to determine the 3D length at any vertice or interpolate the 3D length between any vertices of your line you can download the Add-In I wrote from this post.
Thanks @Richard Fairhurst that is another approach to getting the 3D length updated which is great if you shift lines, be good to look at how you did this so I could potentially use the same method to update start and finish point Z values