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
Solved! Go to Solution.
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:
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_alldct = [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_alld_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):
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
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.
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
looks like this approach may help.
Creating points at a specific distance from other points along a line
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:
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_alldct = [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_alld_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):
interesting implementation Xander and some great points raised.
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.
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
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.
Thanks so much, I will try this now