topic Re: average of ice edge lines (distance from coast) in Spatial Data Science Questions

average of ice edge lines (distance from coast)

AarielRocha1 — Mon, 25 Jan 2016 21:54:37 GMT

I have developed landfast sea ice extents that can be visualized as polygons, lines or rasters. Now I need to calculate monthly and yearly averages and it is much more difficult to create an average line from a set of polylines than I had imagined.

I first tried to use raster world and the Euclidean Distance tools to determine the distance from shore, but since the ice lines and the shoreline are so complex I needed a way to measure the distance from the same point on land or the averages are meaningless.

so, I developed a set of tangent lines and intersected them with my ice edge lines to get a set of points. I have been able to use spatial join to get the distances from shore and then add field and field calculator to get average distances along each of my tangent lines.

My plan is to then plot the average distances for each tangent line as a point and then connect the points to show the 'average' of my original ice edge lines for each month and year.

I think that I need to use linear referencing tools to get past this latest hurdle, so I created routes from my set of tangent lines. Now I need to create points that are my 'average distance' from shore but I am at a loss. I have minimal experience scripting and would appreciate help using the geoprocessing tools in modelbuilder.

Thanks,

Aariel

Re: average of ice edge lines (distance from coast)

DanPatterson_Retired — Mon, 25 Jan 2016 21:59:18 GMT

Now you did say you had lines radiating out to those lines... just a thought Splitting lines at a specified distance or percentage—Help | ArcGIS for Desktop

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Tue, 26 Jan 2016 18:14:19 GMT

Thanks for the response Dan.

I have already split my tangent lines at the intersections with each of my ice edge lines. That is how I got the lengths (=distance from shore) for each date that I needed to get the average distances.

Now I need to get the endpoints from lines that are the average distance from shore.

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Tue, 26 Jan 2016 19:27:33 GMT

looks like this approach may help.

Creating points at a specific distance from other points along a line

Re: average of ice edge lines (distance from coast)

DarrenWiens2 — Tue, 26 Jan 2016 19:39:11 GMT

I believe you can follow these steps (sorry, can't test right now):

1.) Create Routes from your tangent lines

2.) Intersect the ice extent lines with the tangent lines to get intersection points (specify point as output type)

3.) Locate Features Along Routes (features = your intersection points) to get measures for each point

4.) Summary Statistics to get stats for each tangent line

Re: average of ice edge lines (distance from coast)

XanderBakker — Sun, 12 Dec 2021 04:35:38 GMT

OK, I gave it a shot and this is what I got so far...

I took a random part of the globe:

... OK, I admit I georeferenced the image that was attached.

Digitized the coast line, and 4 ice edges (no to mucho detail) and started coding.

The process consists of the following steps:

get the coast line polyline
create a featureclass to store the routes (ice edges where each vertice has a m-value representing the relative position on the coast line)
loop through the ice edges and fill the m-enabled featureclass with the ice edge routes
create a table with events and fill with YEAR values and M values
execute MakeRouteEventLayer_lr to generate a layer with point events based on the ice edge routes and the event table
create a nested dictionary with the event points and relevant information
get min, max and mean statistics from the nested dictionary
... generate polylines and write the resulting geometries and statistics to the output featureclass

Example area showing the ice edges with the event points created on the lines with the m value (the relative position on the coast line).

... and the resulting MIN, MAX and MEAN lines (comparing the event points on the ice edge lines based on the m-values and distance to the coast line) :

Code used:

import arcpy

def main():
    import os

    # input and output featureclasses
    # fc_coast = r'D:\Xander\GeoNet\IceEdges\data.gdb\coast'
    fc_coast = r'D:\Xander\GeoNet\IceEdges\data.gdb\coast_simple'
    fc_ice = r'D:\Xander\GeoNet\IceEdges\data.gdb\ice_edges'

    fld_id = 'YEAR' # should be unique long
    fld_stats = 'Stats'
    fld_dist = 'Distance'
    fld_m_value = 'M_Value'
    fld_loc_err = 'LOC_ERROR'

    fc_routes = r'D:\Xander\GeoNet\IceEdges\data.gdb\routes_v04'
    tbl_events = r'D:\Xander\GeoNet\IceEdges\data.gdb\tbl_events_v04'
    fc_events = r'D:\Xander\GeoNet\IceEdges\data.gdb\events_v04'
    fc_res = r'D:\Xander\GeoNet\IceEdges\data.gdb\result_v04'

    # settings
    ws, fc_name = os.path.split(fc_routes)
    arcpy.env.workspace = ws
    arcpy.env.overwriteOutput = True
    steps = 1000

    # get the coast line (asume first feature in coast fc)
    coast_line = arcpy.da.SearchCursor(fc_coast,('SHAPE@')).next()[0]

    # create output polyline zm featureclass
    sr = arcpy.Describe(fc_coast).spatialReference
    arcpy.CreateFeatureclass_management(ws, fc_name, "POLYLINE", None, "ENABLED", "DISABLED", sr)

    # add field
    arcpy.AddField_management(fc_routes, fld_id, "LONG")

    # loop through ice edges and create dct with polyline m values
    dct_polyline_m = {}
    with arcpy.da.SearchCursor(fc_ice,('SHAPE@', fld_id)) as curs:
        for row in curs:
            year = row[1]
            polyline_ice = row[0]
            lst = []
            for part in polyline_ice:
                for pnt in part:
                    pntg_ice = arcpy.PointGeometry(pnt, sr)
                    pntg_on_line, m_value_coast, dist_to_coast, right_side = coast_line.queryPointAndDistance(pntg_ice, True)
                    pnt_m = arcpy.Point(pnt.X, pnt.Y, None, m_value_coast * 100.0)
                    lst.append(pnt_m)

            polyline_m = arcpy.Polyline(arcpy.Array(lst), sr)
            dct_polyline_m[year] = polyline_m

    # create routes featureclass
    with arcpy.da.InsertCursor(fc_routes, ('SHAPE@', fld_id)) as curs:
        for year, polyline_m in dct_polyline_m.items():
            curs.insertRow((polyline_m, year, ))

    # create table with events
    ws, tbl_name = os.path.split(tbl_events)
    arcpy.CreateTable_management(ws, tbl_name)
    arcpy.AddField_management(tbl_name, fld_id, "LONG")
    arcpy.AddField_management(tbl_name, fld_m_value, "DOUBLE")

    # get list of ID's (years)
    lst_ids = sorted(list(set([r[0] for r in arcpy.da.SearchCursor(fc_ice, (fld_id))])))

    # fill events table
    with arcpy.da.InsertCursor(tbl_events, (fld_id, fld_m_value)) as curs:
        for step in range(steps + 1):
            m = float(step) / float(steps) * 100.0
            for y in lst_ids:
                curs.insertRow((y, m, ))

    # create de route events
    event_props = "{0} POINT {1}".format(fld_id, fld_m_value)
    arcpy.MakeRouteEventLayer_lr(in_routes=fc_routes, route_id_field=fld_id,
                                 in_table=tbl_events,
                                 in_event_properties=event_props,
                                 out_layer="EventPoints", offset_field="",
                                 add_error_field="ERROR_FIELD",
                                 add_angle_field="NO_ANGLE_FIELD",
                                 angle_type="NORMAL",
                                 complement_angle="ANGLE",
                                 offset_direction="LEFT",
                                 point_event_type="POINT")

    # write the events to a fc
    arcpy.CopyFeatures_management("EventPoints", fc_events)

    # generate nested dct with event points and distance to coast (only no error events)
    where = "{0} = 'NO ERROR'".format(arcpy.AddFieldDelimiters(fc_events, fld_loc_err))
    dct_all = {}
    lst_m = []
    with arcpy.da.SearchCursor(fc_events,
                                ('SHAPE@', fld_id, fld_m_value),
                                where_clause=where) as curs:
        for row in curs:
            pntg = row[0]
            y = row[1]
            m = row[2]
            lst_m.append(m)
            d = pntg.distanceTo(coast_line)
            if m in dct_all:
                dct = dct_all
                dct = [pntg, d]
                dct_all = dct
            else:
                dct = {y: [pntg, d]}
                dct_all = dct

    # create list of m values
    lst_m = sorted(list(set(lst_m)))

    # stats
    lst_min, dist_min = getStatPointsFromDct(dct_all, lst_m, lst_ids, "MIN")
    lst_max, dist_max = getStatPointsFromDct(dct_all, lst_m, lst_ids, "MAX")
    lst_mean, dist_mean = getStatPointsFromDct(dct_all, lst_m, lst_ids, "MEAN")

    # create output fc for result
    ws, fc_name = os.path.split(fc_res)
    arcpy.CreateFeatureclass_management(ws, fc_name, "POLYLINE", None, "ENABLED", "DISABLED", sr)

    # add fields to result fc
    arcpy.AddField_management(fc_res, fld_stats, "TEXT", None, None, 25)
    arcpy.AddField_management(fc_res, fld_dist, "DOUBLE")

    # write results to output featureclass
    with arcpy.da.InsertCursor(fc_res, ('SHAPE@', fld_stats, fld_dist)) as curs:
        # MIN:
        polyline = arcpy.Polyline(arcpy.Array(lst_min), sr)
        curs.insertRow((polyline, "MIN", dist_min, ))
        # MAX:
        polyline = arcpy.Polyline(arcpy.Array(lst_max), sr)
        curs.insertRow((polyline, "MAX", dist_max, ))
        # MEAN:
        polyline = arcpy.Polyline(arcpy.Array(lst_mean), sr)
        curs.insertRow((polyline, "MEAN", dist_mean, ))


def getStatPointsFromDct(dct_all, lst_m, lst_ids, stat):
    lst_polyline = []
    lst_dist = []
    for m in lst_m:
        if m in dct_all:
            dct = dct_all
            d_min = None
            pnt = None
            for y, lst in dct.items():
                pntg = lst[0]
                d = lst[1]
                if d_min is None:
                    d_min = d
                    pnt_min = pntg.firstPoint
                    d_max = d
                    pnt_max = pntg.firstPoint
                    d_sum = d
                    lst_x = [pntg.firstPoint.X]
                    lst_y = [pntg.firstPoint.Y]
                else:
                    d_sum += d
                    lst_x.append(pntg.firstPoint.X)
                    lst_y.append(pntg.firstPoint.Y)
                    if d < d_min:
                        d_min = d
                        pnt_min = pntg.firstPoint
                    if d > d_max:
                        d_max = d
                        pnt_max = pntg.firstPoint
            if stat == "MIN":
                lst_polyline.append(pnt_min)
                lst_dist.append(d_min)
            elif stat == "MAX":
                lst_polyline.append(pnt_max)
                lst_dist.append(d_max)
            elif stat == "MEAN":
                lst_polyline.append(arcpy.Point(sum(lst_x)/len(lst_x), sum(lst_y)/len(lst_y)))
                lst_dist.append(d_sum/len(lst_x))
    return lst_polyline, sum(lst_dist)/len(lst_dist)


if __name__ == '__main__':
    main()

At second intent I simplified the cost line (less detail, since the more detail and strange forms it has, the bigger the effect on erroneous distribution of the events over the edge lines:

after the simplification of the coast line (with a higher step setting):

Re: average of ice edge lines (distance from coast)

DanPatterson_Retired — Thu, 28 Jan 2016 01:46:54 GMT

interesting implementation Xander and some great points raised.

Re: average of ice edge lines (distance from coast)

XanderBakker — Thu, 28 Jan 2016 02:28:00 GMT

Thanks Dan. I left the median out, since that would raise more discussion related to how to define the median point based on 4 lines in my case. I would prefer to take the two median points and determine the middle between them, but not sure if that would be correct.

Re: average of ice edge lines (distance from coast)

DanPatterson_Retired — Thu, 28 Jan 2016 03:05:32 GMT

Short answer.... take the "middle" between the middle two lines for an even set of lines... take the middle line if the line set is odd

cavaet... the sea ice people in my department do not necessarily share the views expressed here

Perhaps branch this as a discussion

read below only if you really need to

------------------------------------

medians are tricky, for polygons of predictable geometry it is fairly well agreed upon, but for complex geometries, you are talking about a medial axis (skeleton) but still the controversy would remain

Is the median the middle of the values of X and Y?
Do you account for orientation when you are performing the above?
If you only have the bounding points and assumed edges, is the median space the area delineated by the normal axis lines extending inwards to the other side? Where those lines meet, they may form more than one point of intersection. Do you then go back to the beginning and proceed downward again?
It is "easy" for regular polygons and few well-behaving lines.

In your case, you have 4 lines. the question becomes, where and how do you collapse them into a single medial point. I would argue against the average for the very reason that 3 of the 4 lines may be closely spaced and the 4th a some distance. As for any median, what represent the distribution's middleness based upon sequential position rather than arithmetic average. So following those rules that apply to other properties, you could take the middle point in the ordered list of attributes. Should the distribution be even, it is generally agreed that the average of the middle two is used... if odd, the middle value is used (although, there is some "issues" associated with using that definition). A plausible solution would be to take the middle two lines out of the four and "collapse to center line" aka, create a medial axis/skeleton if you consider the linear distribution as simply representing a long polygon.

In short...no one can fault any descriptor of the median, since it is a nebulous quantity for univariate discrete properties, it becomes "middle-ness" for continuous, and when you throw in more than one dimension you often find yourself flipping back and forth between boundaries/lines to points to boundaries/lines to points etc.

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Sun, 31 Jan 2016 00:04:45 GMT

Thanks so much, I will try this now

Re: average of ice edge lines (distance from coast)

ChadKopplin — Sun, 31 Jan 2016 03:12:16 GMT

I just got this, but if you have the coast line and have the mean, median, or average distances in fields, you can use this distance field to buffer the coast line to created polygons, then convert the polygons to lines, then simply split the lines at the ends where you need to, you can choose to use square ends versus rounded ends. With this your lines are tied to the coast line. Just a thought.

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Mon, 25 Apr 2016 01:51:28 GMT

Xander,

I was called away for another project, so I am just now trying to apply your code. Could you explain more how files were set up prior to running the code? I have each of my ice lines in separate feature classes within a geodatabase, but am getting errors that runtime cannot open the ice edge files after line 40 (part of # loop through ice edges and create dct with polyline m values). it seems that I don't have the files in the proper format.

Thanks for your help

Re: average of ice edge lines (distance from coast)

XanderBakker — Sun, 12 Dec 2021 04:35:41 GMT

Hi Aariel Rocha , maybe the easiest way to explain it, is attaching the data so you can have a look yourself. With respect to the data the following is important:

# fc_coast = r'D:\Xander\GeoNet\IceEdges\data.gdb\coast'
fc_coast = r'D:\Xander\GeoNet\IceEdges\data.gdb\coast_simple'
fc_ice = r'D:\Xander\GeoNet\IceEdges\data.gdb\ice_edges'
fld_id = 'YEAR' # should be unique long

The fc "coast" and "coast_simple" are featureclasses that contains a single line describing the coast. The simplified coast line works best.

The ice_edges fc contains all the ice_edges and has a field called "YEAR" that has unique values for each ice edge.

I also included some of the intermediate results and the final result "result_v04" based on the simple coast line.

fc_routes = r'D:\Xander\GeoNet\IceEdges\data.gdb\routes_v04'
tbl_events = r'D:\Xander\GeoNet\IceEdges\data.gdb\tbl_events_v04'
fc_events = r'D:\Xander\GeoNet\IceEdges\data.gdb\events_v04'
fc_res = r'D:\Xander\GeoNet\IceEdges\data.gdb\result_v04

If you have any questions, please let me know.

Kind regards, Xander

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Fri, 29 Apr 2016 23:20:32 GMT

Thank you kindly for your support on this problem. I have been able to apply the script properly, however when I modified the ice lines something strange happened. It seems that when the routes are created there is a sort of 'doubling back' in a straight line that throws of the results. Since the problem began when I modified length of the ice lines by clipping off the ends, I think it may have to do with the step setting?

Re: average of ice edge lines (distance from coast)

XanderBakker — Sun, 01 May 2016 00:57:08 GMT

Hi Aariel Rocha , could you attach a sample of the modified ice line?

There could be several reasons for this behavior. Are the features still single part? If not this could be a reason since maybe the parts have different directions (some pointing north some south).

Another reason may be that there may be several features that have the same "YEAR"value. It is supposed that the field contains unique values (one feature per year). If this is not the case, you may need to dissolve the lines by YEAR and use the unsplit option and check if the features are created correctly (single parts).

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Sun, 01 May 2016 17:39:13 GMT

Xander,

Each of my ice lines is single part, and all year values are unique. However, looking closely at my lines, I see that the offending ones have breaks in them (and thus may behave as a multipart, I suppose). I will fix those and try again.

Thank you very much for your quick response. Your assistance has been instrumental in my progress and is much appreciated.

Aariel

Re: average of ice edge lines (distance from coast)

AarielRocha1 — Sun, 01 May 2016 23:11:31 GMT

That was it! I had clipped my ice edges to an outline of my study site and discovered some topology issues between the coastline and the ice-free areas of the coast. When I preformed the clip, some of the ice edges became fragmented and that threw off the results.

Now I have some nice summary statistics to work with. Thanks for your help (and extra brownie points for working weekends!)

Re: average of ice edge lines (distance from coast)

XanderBakker — Mon, 02 May 2016 01:18:39 GMT

Hi Aariel Rocha , I'm glad you found the reason for the problem. Indeed if you select a single feature and you see several "parts" light up, it is a multi part. Having multi parts will have a negative effect on the script. It does not handle them well, I guess.

I am glad to hear that the script is helping you to analyze the ice lines. If you are ready with the analysis and you are willing (and you are allowed) to share some results, it would be very nice to post them back at this thread.

If there is anything else, just let us know.

Re: average of ice edge lines (distance from coast)

curtvprice — Sun, 12 Dec 2021 04:35:44 GMT

I though this was an interesting problem, and, waffle-head that I am, went at it with an entirely raster approach.

This conveniently avoids issues of topology!
There's a lot of variability in these ice extent lines so I took the approach of averaging distances for equal length reaches of coastline, instead of thousands of tangent lines. You could simplify the coastline just for the allocation processing to straighten it enough to allow more bands if you wanted to split the coastline more.

I think this is a neat example of a problem that can be expressed and solved very simply using map algebra. That said, I am in awe of Xander's code up-thread. Xander probably provided a more accurate solution than this raster approach, pretty as it is.

from arcpy.sa import *
env.extent = "coast_simple"
env.cellSize = 10000
# create a raster of coastline cells
arcpy.PolylineToRaster_conversion("coast_simple", "OBJECTID", "ccell")
# Create coastline cell zones
# find the southernmost cell (this approach would vary with coastline config)
arcpy.RasterToPoint_conversion("ccell", "ccellpt")
arcpy.AddXY_management("ccellpt")
where = "ABS(POINT_Y - {}) < .001".format(arcpy.Describe("ccellpt").extent.YMin)
arcpy.MakeFeatureLayer_management("ccellpt", "cmin", where)
# Create 20 equal-length reaches along coast line
cc = ExtractByMask(CostDistance("cmin", Con(IsNull("ccell"), 1e6, 1)), "ccell")
cc1 = Slice(cc, 20, "EQUAL_AREA")
# separate ocean into areas corresponding to each coastline reach
calloc = EucAllocation(cc1)
# calculate distance from (any) coastline
cdist = EucDistance(cc1)
# get statistics for ice lines along coast for each year
for year in [2008, 2009, 2010, 2011]:
    arcpy.SelectLayerByAttribute_management("ice_edges", "", "YEAR = {}".format(year))
    ZonalStatisticsAsTable(ExtractByMask(calloc, "ice_edges"), "VALUE", cdist, "stats_{}".format(year), "DATA")

Re: average of ice edge lines (distance from coast)

XanderBakker — Tue, 03 May 2016 11:18:06 GMT

Very cool!