Determining if a place (given coordinates) is on land or on ocean

Question

I have some coordinates from the ISS (International Space Station) and I'd like to know whether when the coordinates were recorded the ISS was over land or ocean and I should do this offline, but I'm not sure what approach to use. A part from python standard library, I'm restricted to only using these libraries:

numpy
scipy
tensorflow
pandas
opencv-python
opencv-contrib-python
evdev
matplotlib
logzero
pyephem
scikit-image
scikit-learn
reverse-geocoder

If you know how to do this using some other libraries though, it'd be good anyway.

With this code I get the coordinates and write them to a file:

import logging
import logzero
from logzero import logger
import os
import ephem
import time

dir_path = os.path.dirname(os.path.realpath(__file__))
logzero.logfile(dir_path+"/coordinates.csv")

# Set a custom formatter
formatter = logging.Formatter('%(name)s - %(asctime)-15s - %(levelname)s: %(message)s');
logzero.formatter(formatter)

name = "ISS (ZARYA)"
line1 = "1 25544U 98067A   18282.18499736  .00001222  00000-0  25998-4 0  9992"
line2 = "2 25544  51.6418 170.6260 0003545 261.4423 234.4561 15.53790940136242"
iss = ephem.readtle(name, line1, line2)

iss.compute()

latitude = iss.sublat
longitude = iss.sublong

# Save the data to the file
logger.info("%s,%s", latitude, longitude )

Do you guys have any idea? Thanks in advance.

Answer 1

global-land-mask from Karin Todd is extremely easy to use and efficient:

from global_land_mask import globe
print(globe.is_land(49.22, -2.23))
# → True
print(globe.is_land(49.22, -2.25))
# → False

It is available via pip, and its only dependency is numpy.

Answer 2

mpl_toolkits.basemap may be able to help.

from mpl_toolkits.basemap import Basemap
bm = Basemap()   # default: projection='cyl'
print(bm.is_land(99.0, 13.0))  #True
print(bm.is_land(0.0, 0.0)) # False

Docs: here and relevant method below:

is_land(xpt, ypt) Returns True if the given x,y point (in projection coordinates) is over land, False otherwise. The definition of land is based upon the GSHHS coastline polygons associated with the class instance. Points over lakes inside land regions are not counted as land points.

Note: you may need to be careful with the projection that is being used with the Basemap object.

Answer 3

In the end, I was able to solve my problem only using those libraries. I used this website geoplaner to get a rough outline of the oceans shape (it was really really rough, since I did it by hand, but it worked fine for my purposes, I think online there should be some more accurate polygons but I was lost on how to use them).

I did so for each of the oceans and got this (note that the coordinates I used didn't cover the oceans entirely, for example I avoided the Southern Ocean):

    atlanticOcean = [(-24.6,68.5), (25.3,69.8), (5.7,61.4), (4.6,52.2), (-6.3,48.4),
            (-9.45,43.5), (-9.63,37.6), (-6.3,35.5), (-10.5,31.1), (-10.5,28.4),
            (-16.1,24.5), (-17.2,14.7), (-8.2,4.1), (6.3,3.6), (9.9,3.4),
            (9,-1.7), (13.8,-12.6), (11.7,-16.5), (14.5,-22.3), (16.1,-28.67),
            (18.9,-34.5), (18.9,-55.7), (-66,-55.7), (-68.5,-50.4), (-58.6,-39.3), (-48.1,-28.2),
            (-48.1,-25.7), (-41.6,-22.7), (-38.7,-17.4), (-39.5,-13.7), (-36.9,-12.5),
            (-34.9,-10.4), (-35.0,-5.5), (-50,-0.1), (-53,5.5), (-57.2,6.1),
            (-62.8,10.9), (-67.8,10.9), (-74.2,10.8), (-76.9,8.5), (-81.6,9.4),
            (-82.7,14), (-87.4,16.1), (-86.3,21.6), (-90.2,21.7), (-91.2,19.2),
            (-95.7,18.8), (-97.1,25.5), (-91.0,28.9), (-84,29.7), (-82.9,27.3),
            (-80.9,24.9), (-79.3,26.7), (-81.1,31.3), (-75.4,35.2), (-73.8,40.3),
            (-69.6,41.4), (-65.1,43.5), (-60,45.8), (-52.2,47.1), (-54.9,52.9),
            (-44.5,60.1), (-38.8,65.1)]

indianOcean =  [(21.40,-34.15), (27.37,-33.71), (40.03,-15.61), (39.68,-3.50), (51.80,10.16), 
                (58.84,22.26), (65.69,25.18), (71.32,19.83), (77.47,6.86), (80.24,12.53),
                (80.90,15.85), (89.05,22.12), (91.38,22.08), (94.54,17.74), (94.02,16.02),
                (97.00,16.82), (98.19,8.33), (100.78,3.18), (94.98,6.29), (105.0,-6.52),
                (118.16,-9.26), (123.52,-11.25), (129.93,-11.08), (128.62,-14.51), (125.89,-3.57),
                 (118.51,-20.37), (113.06,-22.18), (115.26,-34.44), (123.52,-34.88), (130.99,-32.09),
                (137.23,-36.59), (137.50,-66.47), (102.26,-65.79), (85.65,-66.22), (75.01,-69.50),
                (69.04,-67.67), (54.18,-65.76), (37.48,-68.65)]

Now, the Pacific Ocean was more complicated, because it stretches over two sides of a map and you can have two successive points with longitude of -179 and 179 which lead to this polygon to be badly represented in a xy plane. What I did was to divide it in two and so I got this:

pacificEast = [(149.9,-37.8),(153.9,-28.5),(143.2,-11.5),(152.1,-0.9),(127.9,5.7),
                (122.9,23.8),(123.4,31),(128.9,33.7),(129.8,29.4),(141.6,35),
                (142.8,41),(148,43.3),(144.6,45.5),(146.2,49.3),(144.9,54.2),
                (136.8,55.2),(143.1,59.1),(153.7,59.2),(159.4,61.6),(160.3,60.5),
                (161.4,60.3),(155.4,57),(156.6,50.3),(160.8,52.8),(164.1,55.8),
                (163.8,58.1),(167.3,60.1),(170.7,59.8), (179.9,-77.1),
                (166.4,-77.1), (173.8,-71.8), (142.9,-66.8), (146.9,-44.8)]

pacificWest = [(-179.9,62.2),(-179.7,64.7),
                (-177.3,65.3),(-173.6,63.4),(-166,62.2),(-165.8,60.9),(-168.4,60.4),
                (-166.6,58.9),(-158.5,57.8),(-153.1,57),(-144.8,59.9),(-136.1,56.9),
                (-131.7,51.9),(-125.2,48.4),(-124.5,44.6),(-124.4,40.7),(-117.6,32.7),
                (-110.7,23.2),(-105.8,19.7),(-96.1,15.3),(-87.9,12.4),(-83.7,7.3),
                (-78.7,6.1),(-80.2,0.9),(-82.2,-0.6),(-81.2,-6.3),(-76.7,-14.4),
                (-70.4,-18.9),(-73.7,-36.7),(-76,-46.2),(-75.1,-53),(-73.4,-55.1),
                (-66.6,-56.3),(-64.6,-55),(-59.6,-63.4),(-68.4,-65.7),(-75.8,-72.2),
                (-98.6,-71.8),(-126.8,-73.2),(-146.8,-75.7),(-162.6,-78.4),(-179.9,-77.1)]

As I understand, using matplotlib you can create the polygons from the verteces (the lists of coordinates) with path and then you can use the contains_point() function to check whether the point is in either of the polygons (therefore it's in 'ocean') or not (in 'land'):

    p1 = path.Path(atlanticOcean)
    p2 = path.Path(indianOcean)
    p3 = path.Path(pacificEast)
    p4 = path.Path(pacificWest)

    target = [(lon, lat)]

    result1 = p1.contains_points(target)
    result2 = p2.contains_points(target)
    result3 = p3.contains_points(target)
    result4 = p4.contains_points(target)

    # if target is in one of the polygons, it is in ocean
    if result1==True or result2==True or result3==True or result4==True: 
        print("In Ocean")     
    else:
        print("Land")

the lon and lat variables for me where the ones of the ISS which I calculated with the program in my question.

Answer 4

I stumbled accross this problem when triying to classify windmills as on- and offshore. All the given solutions here were too inaccurate for my case, because windmills positioned near the shoreline or in lakes like the IJsselmeer in the netherlands tended to be missclassified. So I came up with a own solution, based on overpass: It searches for all the different types of waterbodies in a given country in OSM and matches the given coordinate to them (or not). It is much slower than the other solutions proposed here and you have to know the the country in which you want to classfiy stuff, but on the plus side it is much more accurate. The areaID is the OSM-ID plus 3600000000, for some reason only then can it be used as searching area in overpass

from geopy.geocoders import Nominatim
countries = {}


def getAreaID(_country_name, asArea=True):    
    # Geocoding request via Nominatim
    geolocator = Nominatim(user_agent="country_compare")    
        
    try:
        geo_results = geolocator.geocode(_country_name, exactly_one=False, limit=3)
    except:
        raise Exception("got an error from Nominatim while getting AreaID of current Country.")

    # Searching for relation in result set
    for r in geo_results:
        if r.raw.get("osm_type") == "relation":
            country = r
            break
    
    # Calculating area id
    area_id = int(country.raw.get("osm_id"))
    if asArea:
        area_id += 3600000000
    
    return area_id

def getCountriesWaterBodies(_countryName):
    if _countryName in countries:
        return countries[_countryName]
    else:
        areaID = getAreaID(_countryName, True)
        eezAreaId = getEEZAreaID(_countryName, True)
        print("load water bodies of {0}".format(_countryName))

        if eezAreaId is not None:
            request = '''            
                [out:json][timeout:10000];  
                area({0})->.country;
                area({1})->.eez;    // exclusive economic zone  
                (            
                   way["natural"="water"](area.country);
                   way["natural"="water"](area.eez);
                   relation["natural"="water"](area.country);  
                   relation["natural"="water"](area.eez);  
                   way["boundary"="maritime"](area.country);
                   way["boundary"="maritime"](area.eez);
                   relation["boundary"="maritime"](area.country);
                   relation["boundary"="maritime"](area.eez);
                   relation["place"="sea"](area.eez);
                   relation["place"="sea"](area.country);
                );        
               out geom;
            '''.format(areaID, eezAreaId)
        else:
            request = '''            
                [out:json][timeout:10000];  
                area({0})->.country;                
                (            
                   way["natural"="water"](area.country);
                   relation["natural"="water"](area.country);                       
                   way["boundary"="maritime"](area.country);
                   relation["boundary"="maritime"](area.country);                   
                   relation["place"="sea"](area.country);
                );        
               out geom;
            '''.format(areaID)

        overpass_url = "http://overpass-api.de/api/interpreter"
        response = requests.get(overpass_url, params={'data': request})
        if response.status_code == 200:
            # convert result to GeoDataFrame
            data = response.json()
            geojson = osm2geojson.json2geojson(data)
            geojsonDF = gpd.GeoDataFrame.from_features(geojson)
            countries[_countryName] = geojsonDF
            return geojsonDF
        else:
            return None


def isCoordinateInWater(_lat, _lon, _country):
    countryGDF = getCountriesWaterBodies(_country)
    if countryGDF is None:
        raise Exception("The Countries water bodies could not be loaded")

    # store the coordinate in question into a dataframe
    df = pd.DataFrame({'y': [_lat], 'x': [_lon]})
    pointGDF = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df['x'], df['y']))

    # search for water bodies at this point
    pointInPolys = sjoin(pointGDF, countryGDF, how='left')
    intersectingWaterBodies = pointInPolys[~pointInPolys['index_right'].isnull()]
    isInWater = intersectingWaterBodies.shape[0] > 0
    return isInWater
    
def getEEZAreaID(_countryName, asAreaID=True):
    id = None

    if _countryName == "netherlands":
        id = 3893531
    elif _countryName == "belgium":
        id = 7645548
    elif _countryName == "italy":
        id = 4769816
    elif _countryName == "germany":
        id = 3649463
    elif _countryName == "poland":
        id = 9942670
    elif _countryName == "norway":
        id = 13955915
    elif _countryName == "france":
        id = 2500903

    if id is not None and asAreaID:
        id += 3600000000
    return id