0

Currently I'm working with a dataset that contains routes around the sea but some of them either cross land or are on land (due to the fidelity of the data being quite low). I have been using the great https://github.com/toddkarin/global-land-mask tool from toddkarin to find which of the coordinates I have are on land so I can discard them (eventually I'll may find a way of moving them to the nearest point at sea).

My current problem is that I need to find a way of finding if a line (given any two coordinates) crosses land (think of an island between two point in the sea).

My area of operation is the entire globe and I am using WGS84 if that changes anything. I have some very basic experience with matplotlib/Basemap but I'm not at all confident with it and I'm struggling to find where to start with this. Do I try to plot each coordinate along the line at a given distance/resolution and then use Todd's tool or is there a more efficient way?

Thanks in advance for any assistance. I've done a lot of digging and reading before posting but haven't found what I think I need.

I need the tool to be in python ideally but if I need to call another language/library/exe that can give me a True/False output that's good too.

SFSDev
  • 23
  • 4

2 Answers2

1

A possible tool available in Python to perform these sorts of operations is Shapely.

If you're able to extract the polygon data of island and other masses, then you could use Shapely to perform an intersection test (see Line vs. Polygon Intersection Coordinates). This will work for checking intersections between points, lines and arbitrary polygons.

The quick and dirty way is as you propose yourself, to discretize the line between the two points and check each of these.

Thomas
  • 36
  • 4
1

Thanks to some help from the answer here I came up with the following which is now working.

How to find all coordinates efficiently between two geo points/locations with certain interval using python

from global_land_mask import globe

def crosses_land(x1,y1,x2,y2):

    # your geo points
    #x1, y1 = 13.26077,100.81099
    #x2, y2 = 13.13237,100.82993
    # the increment step (higher = faster)
    STEP = 0.0003

    if x1 > x2:           # x2 must be the bigger one here
        x1, x2 = x2, x1
        y1, y2 = y2, y1

    for i in range(int((x2-x1)/STEP) + 1):
        try:
            x = x1 + i*STEP
            y = (y1-y2)/(x1-x2) * (x - x1) + y1
        except:
            continue
        is_on_land = globe.is_land(float(x), float(y))
        #if not is_on_land:
            #print("in water")
        
        if is_on_land:
            #print("crosses land")
            
            return True
print(crosses_land(x1,y1,x2,y2))
SFSDev
  • 23
  • 4