3

In python I have a shapely polygon (which defines the coastline of the UK), given a point within the UK I would like to find the closest point on the polygon (coastline) to that point. I have been trying to modify the code posted here, but the resulting nearest point I get is incorrect. What should I change?

from shapely.geometry import Point
import geopandas as gpd
from cartopy.io import shapereader as shp_r
from shapely.ops import nearest_points
import pandas as pd
import geoplot as gplt
import matplotlib.pyplot as plt


# Random point in the U.K (close to the south coast, Worthing)
long, lat = 50.907262, -0.420036
p0 = Point(long, lat) 

#10m resolution country data
resolution = '10m'
category = 'cultural'
name = 'admin_0_countries'
shpfilename = shp_r.natural_earth(resolution, category, name)
world_pd = gpd.read_file(shpfilename)

# Filter out United Kingdom
country_name = 'United Kingdom'
country_pd = world_pd[world_pd['ADMIN'] == country_name]


# Extract mainland as single Polygon
uk_main_landmass = country_pd['geometry'].iloc[0].geoms[1]
print (type(uk_main_landmass)) #Out: shapely.geometry.polygon.Polygon

#Plot polygon (UK)
m = gpd.GeoSeries(uk_main_landmass.boundary)
# m.plot()

#Calculate the nearest point on the polygon to the input point (p0)
p1, p0 = nearest_points(uk_main_landmass.boundary, p0)
print (list(p1.coords)) #Result is near Dover rather than on South coast?

nearest_lat, nearest_long   = (p1.coords)[0] #Long and lat are swapped? 

#Create geopandas dataframes for random starting point
df_start = pd.DataFrame(
    {'Latitude': [long],
     'Longitude': [lat]})
start = gpd.GeoDataFrame(
    df_start, geometry=gpd.points_from_xy(df_start.Longitude, df_start.Latitude))

#Create geopandas dataframes for nearest point on polygon (coastline)
df_nearest = pd.DataFrame(
    {'Latitude': [nearest_long],
     'Longitude': [nearest_lat]})
nearest = gpd.GeoDataFrame(
    df_nearest, geometry=gpd.points_from_xy(df_nearest.Longitude, df_nearest.Latitude))

#Shows input point (red) and "closest point" (blue) on coastline
ax = gplt.polyplot(country_pd, figsize=(12, 12))
start.plot(ax = ax, color = 'red')
nearest.plot(ax = ax, color = 'blue')
plt.show()
3

Guess what? ;)

>>> p1, p0 = nearest_points(uk_main_landmass.boundary, p0)
>>> p0.coords.xy
(array('d', [50.907262]), array('d', [-0.420036]))
>>> p1.coords.xy
(array('d', [1.3844507170000497]), array('d', [51.151678778000075]))

Yes, you just inverted your (lat, lon) coordinates when defining p0!

This will definitely work:

(...)
lat, long = 50.907262, -0.420036 # instead of long, lat = 50.907262, -0.420036
p0 = Point(long, lat) 
(...)

Now they seem much closer, don't they?

>>> p1, p0 = nearest_points(uk_main_landmass.boundary, p0)
>>> p0.coords.xy
(array('d', [-0.420036]), array('d', [50.907262]))
>>> p1.coords.xy
(array('d', [-0.420036]), array('d', [50.802069403000075]))

Then fix your nearest point coords:

# No, they are no more swapped after the fix:
# it was: nearest_lat, nearest_long = (p1.coords)[0] 
nearest_long, nearest_lat = (p1.coords)[0] 

You also need to fix both your df_start and df_nearest GeoDataFrame as:

df_start = pd.DataFrame(
    {'Latitude': [lat], # it was long
     'Longitude': [long]}) # it was lat
(...)
df_nearest = pd.DataFrame(
    {'Latitude': [nearest_lat], # it was nearest_long
     'Longitude': [nearest_long]}) # it was nearest_lat

Resulting in:

Results

Notice that even the shape of the country improved.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.