Finding the distance of stations from a point without crossing land?

[SophiaPiper]

Hi, 

I am working in ArcGIS pro (v3.0). I have a feature class with points that are stations in a bay. I am trying to find the distance between each station and a point-source of pollution in a separate feature class. The distance cannot cross land.

Most of the tools that I have looked at only finds the most cost-effective or optimal path between the points of a source feature (ie between stations), not the distance of each station from a pollution source. 

Ideally, the product should be a field in the stations attribute table  that contains the distance of each station from the pollution source (without crossing land).

I know that this is doable in ArcGIS, I am just not sure how. I would appreciate any help that I can get. 

Thank you!

[HeatherBell]

Hi @SophiaPiper,

I haven't tried this myself but it caught my interest because of what you are actually trying to do RE: finding the distance to the point-source of pollution incidents. Have you tried to use this workflow that is written in an Esri Blog which accounts for barriers in distance calculations? Anyway, following with interest. 

Best wishes,

Heather. 

[JohannesLindner]

If you're using a raster based approach, it's probably enough to create a cost raster where land has a very high cost.

If you have a polygon feature class representing land, you can also use this script:

land_fc = "Land"
source_fc = "Source"
station_fc = "Stations"
distance_to_coast = 50  # to separate the paths from the land polygons
                        # can be 0
                        # should be less than the distance of the stations to land

output_folder = "memory"
output_name = "ShortestPath"



def get_shortest_line(start_point, end_point, barrier_polygons, max_iter=50):
    iter_count = 0
    do_it_again = True
    line = arcpy.Polyline(arcpy.Array([start_point.firstPoint, end_point.firstPoint]))
    while do_it_again and iter_count < max_iter:
        # loop over the line segments
        line_segments = [arcpy.Polyline(arcpy.Array([line[0][i-1], line[0][i]])) for i in range(1, len(line[0]))]
        for i, segment in enumerate(line_segments):
            # create the array to store the updated segment's vertices
            new_vertices = [segment.firstPoint]
            # loop over the barriers
            for barrier_polygon in barrier_polygons:
                # cut barrier polygon with segment
                try:
                    barrier_split = barrier_polygon.cut(segment)
                except:
                    continue
                # get the smallest half and its vertex that is furthest from the line
                # -> shortest way around the barrier
                try:
                    small_split = sorted(barrier_split, key=lambda s: s.area)[0]
                    vertices = [arcpy.PointGeometry(v) for v in small_split[0]]
                    furthest_vertex = sorted(vertices, key=lambda v: v.distanceTo(segment))[-1]
                except:
                    continue                
                new_vertices.append(furthest_vertex.firstPoint)
            # update the segment
            new_vertices.append(segment.lastPoint)
            line_segments[i] = arcpy.Polyline(arcpy.Array(new_vertices))
        # update the line
        old_line = line
        line = line_segments[0]
        for segment in line_segments:
            line = line.union(segment)
        # if the line didn't change, we're finished
        do_it_again = not line.equals(old_line)
        iter_count += 1
    return line


# dissolve and buffer the land fc
arcpy.env.addOutputsToMap = False
land_dissolve = arcpy.management.Dissolve(land_fc, "memory/LandDissolve", multi_part=False)
if distance_to_coast:
    land_dissolve = arcpy.analysis.Buffer(land_dissolve, "memory/Coastline", distance_to_coast)
arcpy.env.addOutputsToMap = True

# extract the land polygons and the source point
land_polygons = [row[0] for row in arcpy.da.SearchCursor(land_dissolve, ["SHAPE@"])]
source = [row[0] for row in arcpy.da.SearchCursor(source_fc, ["SHAPE@"])][0]

# create the output fc
out_fc = arcpy.management.CreateFeatureclass(output_folder, output_name, "POLYLINE")
arcpy.management.AddField(out_fc, "FID", "LONG")
arcpy.management.AddField(out_fc, "Length", "DOUBLE")

# calculate the shortest path for each station
with arcpy.da.InsertCursor(out_fc, ["SHAPE@", "FID", "LENGTH"]) as i_cursor:
    with arcpy.da.SearchCursor(station_fc, ["OID@", "SHAPE@"]) as s_cursor:
        for oid, station in s_cursor:
            shp = get_shortest_line(source, station, land_polygons)
            i_cursor.insertRow([shp, oid, shp.length])

This will create a Polyline fc with the shortest paths between the source and the stations. You can join it to the stations via Station.OBJECTID = ShortestPath.FID

The tuqouise polygons represent land around a semi-fictional bay, the red dot is the source.