pythoncoordinates

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


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.


Solution

  • 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.