topic Re: I have a Polyline Feature class that I wish to maintain in 3D Editing in Public Transit Questions

I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Wed, 10 Jun 2015 01:09:46 GMT

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

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Wed, 10 Jun 2015 01:18:35 GMT

Do you know the "Add Surface Information (3D Analyst)" tool? It will add this information to a polyline based on a surface.

For a polyline it will add:

3D distance of the line along the surface.
Minimum, maximum, and mean elevation and slope of the line along the surface.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Wed, 10 Jun 2015 03:05:17 GMT

Thanks Xander,

But Sadly that only finds Max and Min plus Mean Z it will not interpolate the Z at the beginning and end of the Polyline which is what i need. generally the Max and Min are not always at the start or finish of each line.

I have got a little further with using ArcScene in that if I have a very small surface model (250 hectares in size) it seems to be a lot more stable in that environment. Strictly speaking I dont even wish to render the surface model in 3d anyway . a 2D view to digitize is fine, I just want the extract Z heights as above.

I could break my surface model into hundred of tiles but that also seems a bit of a compromise. especially when trying to associate the Polyline with the correct surface model under it.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Wed, 10 Jun 2015 03:26:20 GMT

Personally I would be more interested in knowing the elevation difference in the entire trajectory then only at the start and end positions, especially in cases like the one below:

However, I you want to extract start and end height, determine the difference between them and divide it by the length to obtain the slope, this could be scripted with some Python code. This is much faster than editing in ArcScene in a 3D environment.

You may also want to consider editing in ArcGIS Pro (if you have access to it) since it has enhanced performance for 3D editing workflows.

Guidelines for editing features with base height properties in 3D—Help | ArcGIS for Desktop

http://desktop.arcgis.com/en/desktop/latest/guide-books/extensions/3d-analyst/about-setting-feature-heights-z-s-for-3d-e…

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Wed, 10 Jun 2015 03:46:51 GMT

I know what you mean, the truth is engineers place the line at strategic locations to minimize the possibility of creating the variance that you describe, and generally these sections represent the start and finish height of each section of road. to create an overall gradient between the two locations, the peice in the middle is just the effort amount to get to the finished requirement, but we first need to know that the gradient is not too steep overall, the balancing (or lack thereof)of cuts and fills is really the nuts and bolts of getting the grade from one point to the next right.

I have Arc Pro, but need to come up to speed with that too As 3D editing is not so obvious as yet.

I have used the Guidelines for editing features with base height properties in 3D—Help | ArcGIS for Desktop

And yet still am up against it.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Wed, 10 Jun 2015 11:32:05 GMT

Sounds logical what you are explaining. Maybe a field calculator script that reads out the start and end point, reads the elevation from a defined source for those points and calculates the slope between start and end point could be a solution. I will look into that and let you know what I find.

For ArcGIS Pro there is a post by David Reeves with lots of useful resources here: Introduction to ArcGIS Pro for GIS Professionals (1.1) Class Resources

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Sun, 12 Dec 2021 07:02:52 GMT

In case you want to use the FieldCalculator you could do this:

Select Python as parser
Activate the option "Show Codeblock"
Paste the code below in the codeblock
change the path to your DEM (see line 3 of the code below)
Specify "getSlope( !SHAPE! )" to assign to your output field (in my case the field " test")

Code block:

import arcpy
# change the path to your DEM
dem = "D:/Xander/GeoNet/Pasture/ned10m34101g3.tiff"

def createExtentForPoint(pnt, cellsize):
    return arcpy.Extent(pnt.X - cellsize, pnt.Y - cellsize, pnt.X + cellsize, pnt.Y + cellsize)

def getRowColExtentFromXY(ext, pnt, cellsize):
    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 getSlope(polyline):
    cellsize = getRasterCellSize(dem)
    pnts = [polyline.firstPoint, polyline.lastPoint]
    elevs = []
    for pnt in pnts:
        ext1 = createExtentForPoint(pnt, cellsize)
        ext2 = adaptExtentUsingRaster(ext1, dem, cellsize)
        np = createNumpyArrayUsingExtent(dem, ext2, cellsize)
        row, col = getRowColExtentFromXY(ext2, pnt, cellsize)
        z = np.item(row,col)
        elevs.append(z)
    zdif = elevs[1] - elevs[0]
    slope = float(zdif / polyline.length)
    return slope

This was the result on my test set:

You can use the numeric format of the field to represent the calculated fraction as percentage.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Wed, 10 Jun 2015 21:36:21 GMT

Great Stuff @Xander Bakker

I must get into the nut and bolts of geometry using Python myself, it will save me plenty of heartache.

Populating the 18000 records the first time took a bit but looks to be plain sailing running the code on small subsets. as will be the case hourly from now on. So for that, awesome effort.

I was also thinking if I should differentiate between up and down slope (-/+ ve) gradients mainly because I just am lazy and don't want to make a symbol set with both negative and positive values that looks the same color range in both directions.

I might also just populate the Z values of the start and finish point plus (maybe) consider the 3d Length of the line between the start and finish point as it may impact the overall gradient a little . I am not sure reading what you have done thinks in a 3d sense or planimetric sense for the distance between start and finish points.

I also may need to do some error trapping if the polyline falls outside the extent of the DEM. (No such thing as a perfect world especially if it is displayed as incomplete data as is the case here for me.).

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Wed, 10 Jun 2015 22:09:29 GMT

Hi, I'm glad it works for you. You are right, there is no error checking at all in the code I provided.

The slope is calculated by dividing the height difference (end elevation - start elevation) by the length (2D) of the polyline. If you would use the 3D length of the polyline, this would affect the slope in a incorrect way.

For the symbology you can also add another column and calculate the absolute slope and use that for symbology to avoid creating two ranges up- and downhill unless you want to use the symbology to distinguish between up- and downhill.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Wed, 10 Jun 2015 22:43:43 GMT

Yes Thanks Again Xander Bakker

I realise the 3D length will not provide the "right" result but the difference in altitude of the first point vs. the last point in each segment will have a true length between them i.e. not in a flat 2D sense length. I was thinking of just applying basic trig as follows: (But only If I can be sure that the Z has not been accounted for already in this way?)

i.e.

A²+ B² = √C

or

(zdif * zdif ) + (polyline.length * polyline.length) = (Squareroot 3dlength)

Re your Absolute Slope, Thats a great approach, I was going to take the easy way, Reselect all Gradient less than zero and then multiply selected values by negative 1

And as for the lines outside of the DEM I have a similar selection process as above... I have a bounding polygon extent of the DEM, so I perform a select by location of Road features that are completely contained by the bounding polygons, then only perform the calculation on those records...

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Thu, 11 Jun 2015 03:18:51 GMT

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.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

RichardFairhurst — Thu, 11 Jun 2015 03:43:18 GMT

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.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Thu, 11 Jun 2015 04:07:50 GMT

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

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Thu, 11 Jun 2015 04:26:14 GMT

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

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

RichardFairhurst — Thu, 11 Jun 2015 05:24:22 GMT

The tool is written in VB.Net and I have not published the code, just the add-in. I will be working out a process to add Z values to vertices of pipelines soon. A GIS tech is plotting the pipelines and adding a starting Z and ending Z from our improvement plans into two double fields. She is not editing the actual Z values.

I will write a tool to add those Z values to the start and end points and interpolate the Z values for any vertices between using a constant rate of ascent or descent based on the 2D lengths between vertices. For pipelines this should be fine since they do not follow any surface model and are only meant to change vertical direction where elevations are given. Generally the pipelines are being broken where significant changes of direction occur, but only if an elevation is indicated. True curves are being used for pipes that turn corners, so for those I usually only have to apply the start and end Z values.

Anyway, I was waiting for a significant number of pipes to be plotted before developing the tool, but I think that point has been reached. I have to develop my tools in .Net, because we use true curves, which python geometry methods do not support and would corrupt. My tools also have to support both versioned and non-versioned SDE editing and non-SDE editing.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Thu, 11 Jun 2015 16:31:32 GMT

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

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Thu, 11 Jun 2015 20:11:02 GMT

Sounds great, keep me informed if you don't mind I am always looking at ways to improve these processes. I hadn't taken the similarity to pipelines as to what I am trying to achieve but their is certainly a lot of "shared fabric" in a 3D sense.

Thanks again

Mat

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Thu, 11 Jun 2015 20:27:19 GMT

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

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

XanderBakker — Sun, 12 Dec 2021 07:02:55 GMT

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.

Re: I have a Polyline Feature class that I wish to maintain in 3D Editing

MathieuBoonen — Wed, 17 Jun 2015 00:11:57 GMT

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