Find nearest location coordinates in Land using python

Question

A geocode api returns no location information for coordinates in ocean/sea. For those records, I would like to find the nearest possible coordinates that has a valid location information (that is closest land coordinates) Below is the code for fetching location information by passing coordinates

import requests    
request_url = "https://api.mapbox.com/geocoding/v5/mapbox.places/{0}%2C{1}.json?access_token={2}&types=country&limit=1".format(lng,lat,key)
response = requests.get(request_url)
output = response.json()

I have no clue in finding the nearest location. I'm also new to Python

Sample output:

{'type': 'FeatureCollection',
 'query': [32.12, 54.21],
 'features': [{'id': 'country.10008046970720960',
   'type': 'Feature',
   'place_type': ['country'],
   'relevance': 1,
   'properties': {'short_code': 'ru', 'wikidata': 'Q159'},
   'text': 'Russia',
   'place_name': 'Russia',
   'bbox': [19.608673, 41.185353, 179.9, 81.961618],
   'center': [37.61667, 55.75],
   'geometry': {'type': 'Point', 'coordinates': [37.61667, 55.75]}}],
 'attribution': 'NOTICE: © 2020 Mapbox and its suppliers. All rights reserved. Use of this data is subject to the Mapbox Terms of Service (https://www.mapbox.com/about/maps/). This response and the information it contains may not be retained. POI(s) provided by Foursquare.'}

Output when the coordinates are ocean:

{'type': 'FeatureCollection',
 'query': [0, 0],
 'features': [],
 'attribution': 'NOTICE: © 2020 Mapbox and its suppliers. All rights reserved. Use of this data is subject to the Mapbox Terms of Service (https://www.mapbox.com/about/maps/). This response and the information it contains may not be retained. POI(s) provided by Foursquare.'}

Answer 1

Different package to try out is reverse_geocoder, which will return the nearest city, state, and country. Seems to be better than the reverse_geocode package.

import reverse_geocoder as rg
coordinates = (29,-84.1),(37,-125) #Both located in the ocean
rg.search(coordinates) 

Output:

[OrderedDict([('lat', '29.67106'),
              ('lon', '-83.38764'),
              ('name', 'Steinhatchee'),
              ('admin1', 'Florida'),
              ('admin2', 'Taylor County'),
              ('cc', 'US')]),
 OrderedDict([('lat', '38.71519'),
              ('lon', '-123.45445'),
              ('name', 'Sea Ranch'),
              ('admin1', 'California'),
              ('admin2', 'Sonoma County'),
              ('cc', 'US')])]

Answer 2

Using Haversine Formula to find the nearest point (eg. city) based on latitude and longitude.

def dist_between_two_lat_lon(*args):
    from math import asin, cos, radians, sin, sqrt
    lat1, lat2, long1, long2 = map(radians, args)

    dist_lats = abs(lat2 - lat1) 
    dist_longs = abs(long2 - long1) 
    a = sin(dist_lats/2)**2 + cos(lat1) * cos(lat2) * sin(dist_longs/2)**2
    c = asin(sqrt(a)) * 2
    radius_earth = 6378 # the "Earth radius" R varies from 6356.752 km at the poles to 6378.137 km at the equator.
    return c * radius_earth

def find_closest_lat_lon(data, v):
    try:
        return min(data, key=lambda p: dist_between_two_lat_lon(v['lat'],p['lat'],v['lon'],p['lon']))
    except TypeError:
        print('Not a list or not a number.')
    
# city = {'lat_key': value, 'lon_key': value}  # type:dict()
new_york = {'lat': 40.712776, 'lon': -74.005974}
washington = {'lat': 47.751076,  'lon': -120.740135}
san_francisco = {'lat': 37.774929, 'lon': -122.419418}

city_list = [new_york, washington, san_francisco]

city_to_find = {'lat': 29.760427, 'lon': -95.369804}  # Houston
print(find_closest_lat_lon(city_list, city_to_find))

Which Yields:

{'lat': 47.751076, 'lon': -120.740135}  # Corresponds to Washington

Let's suppose you got four json answers from mapbox and you saved them in a list:

json_answers = list()  # = []

json_answers.append({'type': 'FeatureCollection', 
'query': [32.12, 54.21],
'features': [{'id': 'country.10008046970720960',
'type': 'Feature',
'place_type': ['country'],
'relevance': 1,
'properties': {'short_code': 'ru', 'wikidata': 'Q159'},
'text': 'Russia',
'place_name': 'Russia',
'bbox': [19.608673, 41.185353, 179.9, 81.961618],
'center': [37.61667, 55.75],
'geometry': {'type': 'Point', 'coordinates': [37.61667, 55.75]}}],
'attribution': 'NOTICE: ...'})

# I changed only the 'coordinates' value for this example
json_answers.append({'type': 'FeatureCollection', 
'query': [32.12, 54.21],
'features': [{'id': 'country.10008046970720960',
'type': 'Feature',
'place_type': ['country'],
'relevance': 1,
'properties': {'short_code': 'ru', 'wikidata': 'Q159'},
'text': 'Russia',
'place_name': 'Russia',
'bbox': [19.608673, 41.185353, 179.9, 81.961618],
'center': [37.61667, 55.75],
'geometry': {'type': 'Point', 'coordinates': [38.21667, 56.15]}}],
'attribution': 'NOTICE: ...'})

# I changed only the 'coordinates' value for this example
json_answers.append({'type': 'FeatureCollection', 
'query': [32.12, 54.21],
'features': [{'id': 'country.10008046970720960',
'type': 'Feature',
'place_type': ['country'],
'relevance': 1,
'properties': {'short_code': 'ru', 'wikidata': 'Q159'},
'text': 'Russia',
'place_name': 'Russia',
'bbox': [19.608673, 41.185353, 179.9, 81.961618],
'center': [37.61667, 55.75],
'geometry': {'type': 'Point', 'coordinates': [33.21667, 51.15]}}],
'attribution': 'NOTICE: ...'})

# The last answer is "null"
json_answers.append({'type': 'FeatureCollection',
'query': [0, 0],
'features': [],
'attribution': 'NOTICE: ...'})

coord_list = []
for answer in json_answers:
    if answer['features']:  # check if ['features'] is not empty
        # I'm not sure if it's [lat, lon] or [lon, lat] (you can verify it on mapbox)
        print(f"Coordinates in [lat, lon]: {answer['features'][0]['geometry']['coordinates']}")
        lat = answer['features'][0]['geometry']['coordinates'][0]
        lon = answer['features'][0]['geometry']['coordinates'][1]

        temp_dict = {'lat': lat, 'lon': lon}
        coord_list.append(temp_dict)

print(f"coord_list = {coord_list}")

point_to_find = {'lat': 37.41667, 'lon': 55.05}  # Houston
print(f"point_to_find = {point_to_find}")
print(f"find_closest_lat_lon = {find_closest_lat_lon(coord_list, point_to_find)}")

Which yields:

{'lat': 47.751076, 'lon': -120.740135}
Coordinates in [lat, lon]: [37.61667, 55.75]
Coordinates in [lat, lon]: [38.21667, 56.15]
Coordinates in [lat, lon]: [33.21667, 51.15]

coord_list = [{'lat': 37.61667, 'lon': 55.75}, {'lat': 38.21667, 'lon': 56.15}, {'lat': 33.21667, 'lon': 51.15}]

point_to_find = {'lat': 37.41667, 'lon': 55.05}

find_closest_lat_lon = {'lat': 38.21667, 'lon': 56.15}

Answer 3

Here is an unoptimized solution.

What's going on under the hood of the function:

1. Run a GeoPy reverse look-up on a point.
2. If the point is found, return its country name.
3. If the point is not found, search for the nearest point of land in the world_geometry variable.
4. Perform a reverse lookup on that closest point.
5. Return that point's country name (if it exists) or the locality name (if no country name).

from geopy.geocoders import Nominatim
from shapely.ops import nearest_points

def country_lookup(query, geocoder, land_geometry):
    
    try:
        loc = geocoder.reverse((query.y, query.x))
        return loc.raw['address']['country']
    except (KeyError, AttributeError):
        _, p2 = nearest_points(query, land_geometry)
        loc = geocoder.reverse((p2.y, p2.x)).raw['address']
        if 'country' in loc.keys():
            return loc['country']
        else:
            return loc['locality']
        
# get world (or any land) geometry, instantiate geolocator service
world = gp.read_file(gp.datasets.get_path('naturalearth_lowres'))
world_geometry = world.geometry.unary_union
geolocator = Nominatim(user_agent="GIW")

# Create a column of country names from points in a GDF's geometry.
gdf['country'] = gdf.geometry.apply(country_lookup, args=(geolocator, world_geometry))

The accuracy of the results depends on the accuracy of the land geometry you provide. For example, geopandas's world geometry is pretty good. I was able to find names for all countries except for some of the smallest of the islands in the Bahamas. Those that it could not find were labelled "Bermuda Triangle" by the function, which is good enough for me.

Answer 4

Use reverse_geocode library in python to get nearest city with country.

Example:

import reverse_geocode

coordinates = (-37.81, 144.96), (31.76, 35.21)

reverse_geocode.search(coordinates)

Result:

[{'city': 'Melbourne', 'code': 'AU', 'country': 'Australia'}, {'city': 'Jerusalem', 'code': 'IL', 'country': 'Israel'}]