Creating series of lines offset from starting points (ArcGIS, arcpy)

Question

I've got a series of points representing generalized locations of trawling activity over many years. A single XY location can represent many trawling records. I'm attempting to create a series of lines, one for each record, representing an approximation of the trawling footprint. I have a script which takes each record and location, adds a random offset to create a new point, calculates a second point using a random bearing and distance, and creates a trawl line from that pair of points.

However, some of my created trawl lines either start or end on land. What I need to do now is test each point and if it intersects with a land feature class, create a new point and test again. I just have no idea how to do that within script I've written.

Code used in ArcGIS

fc = r"G:\CreateTrawl.gdb\TestPoints"
csvFile = open(r"c:\temp\ArcGISTrawlLines3.csv","w")
writer = csv.writer(csvFile, delimiter = ',', lineterminator = '\n')
writer.writerow(["ID","X1","Y1","X2","Y2"])
fields = ['Long1', 'Lat1', 'OBJECTID']
with arcpy.da.UpdateCursor(fc, fields) as cursor:
 for row in cursor:
 ID2 = row[2]
 bearing = random.randint(0,360) # Random bearing
 dist = random.randint(1000,3500) # Random distance between 1000 and 3500 m
 signy = random.randint(0,1) * 2 - 1 # Random positive or negative for Y offset
 signx = random.randint(0,1) * 2 - 1 # Random positive or negative for X offset
 xoffset = random.randint(100,900) #number of metres offset in x from original point
 yoffset = random.randint(100,900) #number of metres offset in y from original point
 yoffset2 = math.sin(bearing * deg2rad) * dist #calculation of change in y from Y1
 xoffset2 = math.cos(bearing * deg2rad) * dist #calculation of change in x from X1
 X1 = int(row[0] + xoffset * signx) #X of starting point
 Y1 = int(row[1] + yoffset * signy) #Y of starting point
 X2 = int(X1 + xoffset2) #X of ending point
 Y2 = int(Y1 + yoffset2) #Y of ending point
 writer.writerow([ID2,X1,Y1,X2,Y2])
csvFile.close()
arcpy.XYToLine_management(in_table, out_featureclass, startx_field="X1", starty_field="Y1", endx_field="X2", endy_field="Y2", line_type="GEODESIC", id_field="ID", spatial_reference="PROJCS[')
print "Processing complete"

The attached image shows a set of points in a test file. The black points in the image below are the generalized locations, red lines are the created trawl footprints from the above code. The 8 locations shown represent 37 unique records. The green is a 1:10 mil land feature class which is accurate enough for this task to determine whether the trawl intersects with land.

Many of the 816,066 records in the main file are far enough from land so this isn't an issue for most points but there are enough adjacent to land to necessitate a solution in the code (~1% within 2500 m).

Trawl lines

I was wondering whether a Select by Location could be used (if the point is not selected, it's therefore not intersecting the land). The coordinate system units are metres and so I'm using integer values as the excess precision isn't necessary.

Edit: As noted in a comment to GBG below, I've got a test using geometry working, returning a 1 if the new point is on land. But now I'm stuck on how to recalculate the new point within the Search Cursor. Amended partial code is below.

fc_geom = arcpy.CopyFeatures_management(fc_g1, arcpy.Geometry()) #Creates polygon geometry for land
with arcpy.da.SearchCursor(fcsel, fields) as cursor:
 for row in cursor:
 ID2 = row[2]
 bearing = random.randint(0,360) # Random bearing
 dist = random.randint(1000,3500) # Random distance between 1000 and 3500 m
 signy = random.randint(0,1) * 2 - 1 # Random positive or negative for Y offset
 signx = random.randint(0,1) * 2 - 1 # Random positive or negative for X offset
 xoffset = random.randint(100,900) #number of metres offset in x from original point
 yoffset = random.randint(100,900) #number of metres offset in y from original point
 X1 = int(row[0] + xoffset * signx) #X of starting point
 Y1 = int(row[1] + yoffset * signy) #Y of starting point
 pnt1 = arcpy.Point(X1, Y1)
 test1 = pnt1.within(fc_geom[0]) # Test to check whether point is on land
 #From here, if test1 = 1, how do I get back up to the top of the SearchCursor to 
 #calculate new signx, signy, xoffset, and yoffset?
 yoffset2 = math.sin(bearing * deg2rad) * dist #calculation of change in y from Y1
 xoffset2 = math.cos(bearing * deg2rad) * dist #calculation of change in x from X1
 X2 = int(X1 + xoffset2) #X of ending point
 Y2 = int(Y1 + yoffset2) #Y of ending point
 pnt1 = arcpy.Point(X1, Y1)
 test1 = pnt1.within(fc_geom[0]) # Test to check whether point is on land
 writer.writerow([ID2,X1,Y1,X2,Y2])
csvFile.close()

• FelixIP: Create buffers around point. Clip them using water (or erase using land polygons). Randomly generate points inside and draw your lines from centre to each point

Answer 1

If you have the shoreline as a polygon feature class you should be able to use geometry objects to test to see if your points are within the polygon, if it is true that the point is within the polygon, delete the point and try again.

Here is the help page on geometry for 10.3([https://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-classes/geometry.htm])

• phyllo: Based on that link and various others on StackExchange I've been trying to get the Geometry to work but I'm just not clear enough on the various ways geometry is created have it work. If I use a single XY coord pair, taken from the code above, I can create geometry but it returns 'False' when it's clearly within the polygon. >>> fcpt_geom[0].within(fc_geom[0]) True >>> pnt1 = arcpy.Point(-755286, 89147) >>> pnt1_geom = arcpy.PointGeometry(pnt1) >>> pnt1_geom.within(fc_geom[0]) False
• phyllo: In the first case above (fcpt_geom), the fc was a shapefile with a single point. pnt1_geom was a coordinate pair using the same coordinates as the shapefile point. The first was true for 'within' while the second was false. (Sorry for not putting the python output in a code block but I couldn't get the formatting to work).
• GBG: If you can't figure out how to get the geometry object working consider using a select by location. You would 1. Generate a random line base on your selected point. 2. Run select by location-does the line intersect or cross the boundary of the shoreline. 3. If count returns one, then delete the file and try agian. 4. If count returns 0 merge the the new line with the feature class of good trawl paths. This method will probably be slower than using geometry objects but that may not matter to you.
• phyllo: I figured out the geometry test but now I can't figure out how, within the cursor, if that test fails (i.e. point is on land) to calculate new bearings and variables. It seems I have to return to the previous lines of code but I don't have enough python familiarity to do that. I added the updated code in the original question above.
• GBG: Put the code that creates the random points, and the code that does the geometry check inside a function that you define. Use an if statement and if the geometry test fails, delete the point. Next - use recursion (have the function call itself if the geometry check fails) to keep working on the problem until you get points not on land.

Answer 2

This is what I would do:

1. Change your offsets so that the lines are centered over the actual points, i.e., the point is about in the middle of the line. This will make them more accurate and less likely to hit land too.
2. Make your lines shorter. Think about how close a fishing boat would come to land and make sure they are shorter than that.

• phyllo: The first won't work for most of the points. The points in the original file are gridded points representing many different fishing records from many years (more than 8000 of the points represent 20 or more records, some of them represent hundreds of unique records). This was an attempt to estimate spatial coverage of all operations over time, hence the offsets and random bearings. The distances represent a range of trawling distances.
• phyllo: Though perhaps I could use this script for all those points further from the coast than the longest possible line length, and then treat the remaining points differently.
• RHB: @phyllo2 I guess I wasn't clear. I'm talking about the offsets that you are creating in a rather random way. This is something you are doing that has nothing to do with the original data. If you made these lines so they passed through your points instead of randomly off aways, they would be much less likely to hit land. Especially if you made them an appropriate length.
• phyllo: I've rejigged the code to use the original points as a centre for all lines that a location represents but it results in a rosette that over-represents the centre. As mentioned, these points represent the trawl events but are not the true locations. The original trawl would have occured anywhere within an area of about 2.5 sqkm. This method may work for now but the Geometry approach would probably be better if I could get the .within function to work. Thanks.

Answer 3

#Something like this incomplete code...
fc_geom = arcpy.CopyFeatures_management(fc_g1, arcpy.Geometry()) 
Def MakeRandomPoint():
 bearing = random.randint(0,360) # Random bearing
 dist = random.randint(1000,3500) # Random distance between 1000 and 3500 m
 signy = random.randint(0,1) * 2 - 1 # Random positive or negative for Y offset
 signx = random.randint(0,1) * 2 - 1 # Random positive or negative for X offset
 xoffset = random.randint(100,900) #number of metres offset in x from original point
 yoffset = random.randint(100,900) #number of metres offset in y from original point
 X1 = int(row[0] + xoffset * signx) #X of starting point
 Y1 = int(row[1] + yoffset * signy) #Y of starting point
 pnt1 = arcpy.Point(X1, Y1)
 test1 = pnt1.within(fc_geom[0]) # Test to check whether point is on land
 if test1:
 delete test1
 MakeRandomPoint()
 return #coordinates of the good point get returned...
with arcpy.da.SearchCursor(fcsel, fields) as cursor:
 for row in cursor:
 ID2 = row[2]
 goodCoords = MakeRandomPoint()
 yoffset2 = math.sin(bearing * deg2rad) * dist #calculation of change in y from Y1
 xoffset2 = math.cos(bearing * deg2rad) * dist #calculation of change in x from X1
 X2 = int(X1 + xoffset2) #X of ending point
 Y2 = int(Y1 + yoffset2) #Y of ending point
 pnt1 = arcpy.Point(X1, Y1)
 test1 = pnt1.within(fc_geom[0]) 
 writer.writerow([ID2,X1,Y1,X2,Y2])
csvFile.close()

• phyllo: Thanks very much! I think I have it working and will update post when it's complete.