topic Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation? in Python Questions

Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

PeterWilson — Sat, 11 Dec 2021 23:40:13 GMT

I have converted a feature class (polylines) to a numpy array and have expoded the polylines to vertices using the "explode_to_points" method. I would like to convert the output numpy array to a pandas dataframe that will replicate the line segments that make up the polyline so that I land up with the following columns:

The table above was achieved by using "Split Line At Vertices". The fields "STElev" and "EndElev" are the Start Elevation and End Elevation of the line segments determined using "Calculate Geometry". The field "STElevation" and "EndElevation" is the change in elevation determined by subtracting the elevation from the first "STElev" which in this case is 355.8495.

The output numpy array from converting my feature class (polylines) and exploding the features to vertices is:

The values are:

[(DrainID, X, Y, Z)....] for each vertice

'''
Created on 01 Jul 2017

@author: PeterW
'''
# import modules and site-packages
import numpy as np
import pandas as pd
from pathlib import Path
import arcpy

# set environment settings
arcpy.env.overwriteOutput = True

# check out extensions
arcpy.CheckOutExtension("3D")

def longestflowpath_3d(dtm, lfp_2d):
    """Generate LongestFlowPath 3D
    from LongestFlowPath 2D and DTM"""
    output_gdb = str(Path(lfp_2d).parents[0])
    lfp_3d = "{0}\\{1}".format(output_gdb, "lfp_3d")
    lfp_3d = arcpy.InterpolateShape_3d(dtm, lfp_2d, lfp_3d, vertices_only=True)
    arcpy.FlipLine_edit(lfp_3d)
    return lfp_3d


def area_under_profile(lfp_3d):
    """Determine the Equal Area Slope
    for each Drainage Area's Longest
    Flowpath"""
    arr = arcpy.da.FeatureClassToNumPyArray(lfp_3d, ["DrainID", "SHAPE@X", "SHAPE@Y", "SHAPE@Z"], explode_to_points=True)  # @UndefinedVariable
    print(arr)


def equal_area_slope(dtm, lfp_2d):
    """Determine the Equal Area Slope
    for each Drainage Area's Longest
    Flowpath"""
    lfp_3d = longestflowpath_3d(dtm, lfp_2d)
    area_under_profile(lfp_3d)


if __name__ == "__main__":
    dtm = r"E:\Projects\2016\G113665\geoHydro\DTM2\raw1"
    lfp_2d = r"E:\Projects\2016\G113665\geoHydro\Model02\Results.gdb\LongestFlowPath_2D_Orig"
    equal_area_slope(dtm, lfp_2d)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

My current Python Module

How can I convert the current numpy array to a pandas dataframe that will restructure the vertices (DrainID, X, Y, Z) to (DrainID, STElevation (change in elevation), EndElevation (change in elevation), DistM (Length of segment)?

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

JoshuaBixby — Fri, 14 Jul 2017 13:39:08 GMT

Can you post a few sample lines? It would be much easier to demonstrate a solution if there were a couple of examples for us to test against.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

DanPatterson_Retired — Sat, 11 Dec 2021 23:40:17 GMT

Peter

Haven't got time for a full example, but this will get you thinking.

Consider two lines with 4 points each consisting of an ID, X, Y, and Z field as a structured array (numpy )

The final result shown (dz) is the individual lines.

Note... this is the long way, but the logic is less clear in the short way.

a     # a structured array
array([(1, 0.0, 0.0, 0.0), (1, 1.0, 1.0, 1.0), (1, 2.0, 2.0, 2.0),
       (1, 3.0, 3.0, 3.0), (2, 3.0, 3.0, 4.0), (2, 2.0, 2.0, 5.0),
       (2, 1.0, 1.0, 6.0), (2, 0.0, 0.0, 7.0)], 
      dtype=[('ID', '<i8'), ('X', '<f8'), ('Y', '<f8'), ('Z', '<f8')])

ids = np.unique(a['ID'])   # get the unique IDs for the array
ids
array([1, 2], dtype=int64) # no surprise...

# conceptually you want to extract/slice all the points that share a common id
# A list comprehension will be used as an example of this slicing But...

lines = [a[a['ID'] == i] for i in ids]  # .... this is the long way

# now that you have the points for each line, simply take the last and the
# first point (-1 and 0 respectively) and get the 'Z' values

dz = [ln[-1]['Z'] - ln[0]['Z'] for ln in lines]

dz   # no surprise, 4 points per line Z difference each line is 3
[3.0, 3.0]

# similar process for planar distance
# the whole process can be simplified... but that is for later
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

PeterWilson — Fri, 14 Jul 2017 15:03:58 GMT

Hi Joshua

I've attached a copy of the polylines unsplit as well as a copy of the polylines split with the output fields and values that I'm trying to achieve.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

PeterWilson — Fri, 14 Jul 2017 15:08:28 GMT

Hi Dan, Thanks this will get me started. I'll post shortly what I come up with and if I get stuck.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

JoshuaBixby — Sat, 11 Dec 2021 23:40:20 GMT

I think this will get you what you want, or something close:

from itertools import groupby
import pandas

fields = ["STElev", "EndElev", "STElevation", "EndElevation"]
id_field = "DrainID"

split_lines = arcpy.SplitLine_management("LongestFlowPath_unsplit","in_memory/split")

for field in fields:
    arcpy.AddField_management(split_lines, field, "DOUBLE")
    
with arcpy.da.UpdateCursor(split_lines, [id_field, "SHAPE@"] + fields) as cur:
    for k, grp in groupby(cur, key=lambda x: x[0]):
        shape = next(grp)[1]
        fpz, lpz = shape.firstPoint.Z, shape.lastPoint.Z
        kpz = fpz
        cur.updateRow([k, shape, fpz, lpz, fpz - kpz, lpz - kpz])
        for row in grp:
            shape = row[1]
            fpz, lpz = shape.firstPoint.Z, shape.lastPoint.Z
            cur.updateRow([k, shape, fpz, lpz, fpz - kpz, lpz - kpz])
            
with cur:
    df = pandas.DataFrame.from_records(
        cur, exclude=["SHAPE@"], columns=[id_field, "SHAPE@"] + fields
    )

del cur‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Since I don't work with NumPy much, this approach uses an update cursor to populate the table and then just dumps the populated table directly to a pandas data frame.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

DanPatterson_Retired — Sat, 11 Dec 2021 23:40:23 GMT

back for a bit. Now I don't know if you just wanted the 2D planar distance between the first and last point or a 3D distance along the whole length of each polyline. Let me know.

Check these results against a quick scan of measure tool distances. The results are in sequential order...I didn't want to get too fancy for the output yet. I have kept it verbose, since there is a simpler way of calculating the 3D polyline distance if you need it.

pnt_fc = r'F:\_anaconda\polyline_demo\polyline_unsplit\unsplit\split_pnts.shp'
fc = r'F:\_anaconda\polyline_demo\polyline_unsplit\unsplit\LongestFlowPath_unsplit.shp'

flds = ['OID@', 'SHAPE@X', 'SHAPE@Y', 'SHAPE@Z']
SR = arcpy.Describe(fc).spatialreference
a = arcpy.da.FeatureClassToNumPyArray(fc, 
                                   field_names=flds,
                                   spatial_reference=SR,
                                   explode_to_points=True)


dt = [('ID', '<i4'), ('X', '<f8'), ('Y', '<f8'), ('Z', '<f8')]
a.dtype = dt
ids = a['ID']
x_s = a['X']
y_s = a['Y']
z_s = a['Z']
ids = np.unique(ids)
lines = [a[a['ID'] == i] for i in ids]
dz = np.array([ln[-1]['Z'] - ln[0]['Z'] for ln in lines])
dx = np.array([ln[-1]['X'] - ln[0]['X'] for ln in lines])
dy = np.array([ln[-1]['Y'] - ln[0]['Y'] for ln in lines])
dist = np.sqrt(dx**2 + dy**2)
print("2D distances first point to last point \n{}".format(dist))

# ---- result
2D distances first point to last point 
[ 4645.90189847  1249.26103157  1244.66988595  1458.63110731  1757.96791359
  1676.49347515]‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

PeterWilson — Fri, 14 Jul 2017 19:39:19 GMT

Hi Dan

I actually need the distance between each vertices. I'm trying to re-write the Python module to calculate the "Equal Area Slope" using NumPy and Pandas and not ArcPy "SplitLine_Management" based on my previous post: Extract Polyline Geometry Z and M - Equal Area Slope.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

PeterWilson — Fri, 14 Jul 2017 19:47:22 GMT

Hi Joshua

Thanks for the following, sorry I should have been more clear in my post. I'm trying to re-write an existing Python module based on a previous post Extract Polyline Geometry Z and M - Equal Area Slope where I'd like to replace the use of ArcPy "SplitLine_Management" with NumPy arrays and Pandas DataFrames. I wan't to remove the need for and Advanced Licence for the end user.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

JoshuaBixby — Sat, 11 Dec 2021 23:40:25 GMT

Here you go, I did the splitting myself instead of using the geoprocessing tool. NumPy may be faster, I just took this approach since I am more familiar with it.

from itertools import tee
import pandas

unsplit_lines = "LongestFlowPath_unsplit" # Assumes data loaded in ArcMap
fields = ["STElev", "EndElev", "STElevation", "EndElevation"]
id_field = "DrainID"

SR = arcpy.Describe(unsplit_lines).spatialReference
split_lines = arcpy.CreateFeatureclass_management(
    "in_memory","split","POLYLINE", has_m="ENABLED", has_z="ENABLED", spatial_reference=SR
)
arcpy.AddField_management(split_lines, id_field, "LONG")
for field in fields:
    arcpy.AddField_management(split_lines, field, "DOUBLE")
    
with arcpy.da.SearchCursor(unsplit_lines,["SHAPE@","DrainID"]) as scur:
    with arcpy.da.InsertCursor(split_lines,["SHAPE@",id_field] + fields) as icur:
        for shape, drainid in scur:
            for part in shape.getPart():
                n, n1 = tee(part)
                ptz = next(n1).Z
                for pt1 in n1:
                    pt = next(n)
                    segment = arcpy.Polyline(arcpy.Array([pt, pt1]), SR)
                    icur.insertRow([segment, drainid, pt.Z, pt1.Z, pt.Z - ptz, pt1.Z - ptz])

del scur, icur

df_fields = [id_field] + fields
with arcpy.da.SearchCursor(split_lines, df_fields) as cur:
    df = pandas.DataFrame.from_records(cur, columns=df_fields)

del cur‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

DanPatterson_Retired — Sat, 11 Dec 2021 23:40:28 GMT

Peter ...and from the world of numpy... can be greatly simplified if needed. Values can be tied to segment lengths, but I will save that for a blog.

"""
:polyline_demo.py

"""
import numpy as np
import arcpy

ft = {'bool': lambda x: repr(x.astype('int32')),
      'float': '{: 0.2f}'.format}
np.set_printoptions(edgeitems=3, linewidth=80, precision=3,
                    suppress=True, threshold=100, formatter=ft)


def e_leng(a):
    """Length/distance between points in an array using einsum
    : details in one of my blogs... this can be simplified to 4 lines for 2D
    :-----------------------------------------------------------------------
    """
    def cal(diff):
        """ magic happens for 1D, 2D, 3D and 4D arrays... 
        """
        d_arr = np.sqrt(np.einsum('ijk,ijk->ij', diff, diff))
        d_leng = d_arr.flatten()
        length = np.sum(d_leng)
        return length, d_leng
    # ----
    diffs = []
    a = np.atleast_2d(a)
    if a.shape[0] == 1:
        return 0.0
    if a.ndim == 2:
        a = np.reshape(a, (1,) + a.shape)
    if a.ndim == 3:
        diff = a[:, 0:-1] - a[:, 1:]
        length, d_leng = cal(diff)
        diffs.append(d_leng)
    if a.ndim == 4:
        length = 0.0
        for i in range(a.shape[0]):
            diff = a[i][:, 0:-1] - a[i][:, 1:]
            leng, d_leng = cal(diff)
            diffs.append(d_leng)
            length += leng
    return length, diffs


def report(subs):
    """print out the results"""
    for i in range(len(subs)):
        total, segs = e_leng(subs[i])
        frmt = """
        3D distances along polyline {}
        Segment distances
        {}
        Total = sum of segments ? {}
        """
        print(frmt.format(total, segs, total == np.sum(segs)))

# ---- inputs ----
fc = r'C:\GIS\Geometry_projects\polyline_demo\polyline_demo.gdb\polylines'  # change

flds = ['OID@', 'SHAPE@X', 'SHAPE@Y', 'SHAPE@Z']  # these are horibble
SR = arcpy.Describe(fc).spatialreference
a = arcpy.da.FeatureClassToNumPyArray(fc,
                                      field_names=flds,
                                      spatial_reference=SR,
                                      explode_to_points=True)
dt = [('ID', '<i4'), ('X', '<f8'), ('Y', '<f8'), ('Z', '<f8')]  # these are readable
a.dtype = dt                  # simplify the dtype
ids = np.unique(a['ID'])      # get the unique id values
p_lines = [a[a['ID'] == i] for i in ids]                 # collect polylines
subs = [np.c_[p['X'], p['Y'], p['Z']] for p in p_lines]  # stack coordinates

report(subs)  # print the results
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

DanPatterson_Retired — Sun, 20 Aug 2017 03:36:14 GMT

Peter, did you finish this?

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

PeterWilson — Wed, 23 Aug 2017 13:12:16 GMT

Hi Dan

Unfortunately not yet, I started working through your code and got stuck in understanding your reasoning behind reshaping the array from two dimensional to three dimensional lines 31 to 36. I still have quiet a bit to learn about NumPy and haven't needed to work with anything more than a two dimensional array.

Re: Polyline NumPy Array to Pandas DataFrame: StElevation and EndElevation?

DanPatterson_Retired — Wed, 23 Aug 2017 13:24:25 GMT

Peter, there is a caveat in that function which states ... this can be simplified to 4 lines for 2D ....

That is the full code when needing to determine lengths/distances in higher dimensions....

You need not worry about it unless you are working with 3D data (conventional X,Y,Z) or 4D data(3D shifted by a distance over time)

Plus the 'cal' function uses einsum which is a bit of an overkill to determine distance between points, for instance, that can be written much shorter if one is using polygons that one is normally going to encounter. The methods that you will find on the forum will crawl when you are working with objects with millions of vertices (beyond the scope of ArcMap in any event). I just like to use what works for the cases that I use.

Pythagorean theorem is good enough to determine 2D length if that is all you need.