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):