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.'}

Can you share how `output` looks like? – Ali Beyit Jan 14 '20 at 15:18 — Ali Beyit, Jan 14 '20 at 15:18

score 3 · Answer 1 · edited May 06 '21 at 21:16

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}

Thanks Raphael,this looks good. I guess for this we need to have list of all major cities and its coordinates. But in our business, it can be a small town or anything. When its ocean/sea the features returned in the json response is null. I would like to achieve the nearest coordinates which returns non-zero results — Gang Gang, Jan 14 '20 at 22:32
This function works for all globe, you can search within the precision of mapbox, if it's 1m, than you can find the closest distance within 1m of precision. I called the function find_closest_city() just as an example. If you have `null` in the JSON for the ocean, just don't don't call find_closest_city(). Could you give an example of two points that you want to calculate and the expected result? — Raphael, Jan 15 '20 at 09:39
suppose lat, lng = 50.677168, 0.305130 . This is in English channel. We should get closest coordinates in land which is somewhere 50.734998, 0.239263 Eastbourne, England — Gang Gang, Jan 15 '20 at 09:40
Do you have the coordinates in the water then? I thik you need to do a scan in the globe (or a region you want), save all coords over water and all over land first. I updated my answer with more examples. I guess the `place_type` key will tell you if it's land or other thing. The logic is simple, but it requires more data to work with. — Raphael, Jan 15 '20 at 10:15
Thanks Raphael for the response. I am searching for how to get list of all lat,lng coordinates of coatlines — Gang Gang, Jan 16 '20 at 22:17
You're welcome. If you don't have the coordinates in the water will be difficult to know which part of the land is closer. Suppose you have all coordinates for the land, but sea is always `null`. If your search point is in the Athlantic Ocean, which country is closer? Impossible to know without sea coords. Of course I'm assuming this with the info you provided, maybe you have more info that will help to calculate the coords. Let me know if I can help you further. Good luck. — Raphael, Jan 17 '20 at 08:54

score 1 · Answer 2 · edited May 29 '20 at 13:38

1

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'}]

edited May 29 '20 at 13:38

Jaimil Patel

1,301
6
13

answered May 29 '20 at 12:30

Vinay Kumar Marrapu

11
2

I am getting below error `coordinates = (-37.81, 144.96), (31.76, 35.21) reverse_geocode.search(coordinates) ` UnicodeDecodeError: 'charmap' codec can't decode byte 0x98 in position 429: character maps to – Gang Gang Aug 08 '20 at 16:16
I had to go in here https://pypi.org/project/reverse-geocode/ `pip3 install reverse-geocode` – AnonymousUser Oct 31 '21 at 07:08
To add one do this `coordinates = (-37.81, 144.96),` It needs comma in the end, otherwise it gives error. – AnonymousUser Oct 31 '21 at 07:25

rocksNwaves · Answer 3 · 2021-03-15T20:16:13.397

Here is an unoptimized solution.

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

Run a GeoPy reverse look-up on a point.
If the point is found, return its country name.
If the point is not found, search for the nearest point of land in the world_geometry variable.
Perform a reverse lookup on that closest point.
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.

score 0 · Answer 4 · answered Oct 14 '21 at 01:35

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')])]

Find nearest location coordinates in Land using python

4 Answers4

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

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