Calculating Minimum Distance Between a Point and the Coast

calculating minimum distance between a point and the coast

gDistance(...) returns the minimum Cartesian (Euclidean) distance between the point and feature set provided as arguments. Since your map is in long/lat coordinates, you get distance in "degrees", e.g.

d = sqrt { (long1 - long2)2 + (lat1 - lat2)2 }

where long and lat are in decimal degrees. As was pointed out, this doesn't mean much because converting to planar distance (say, km) depends on where you are. So we need to transform your data into a CRS which is approximately planar in the region of interest. It turns out that the appropriate CRS for Spain is EPSG-2062. The projection string for EPSG-2062 is:

+proj=lcc +lat_1=40 +lat_0=40 +lon_0=0 +k_0=0.9988085293 +x_0=600000 +y_0=600000 +a=6378298.3 +b=6356657.142669561 +pm=madrid +units=m +no_defs 

which has +units=m (meters). So we need to reproject both the point (MAD) and the borders to EPSG-2062.

library(rgeos)
library(maptools)

epsg.2062 <- "+proj=lcc +lat_1=40 +lat_0=40 +lon_0=0 +k_0=0.9988085293 +x_0=600000 +y_0=600000 +a=6378298.3 +b=6356657.142669561 +pm=madrid +units=m +no_defs"
wgs.84 <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

coast <- readShapeLines("ne_10m_coastline",CRS(wgs.84))
MAD <- readWKT("POINT(-3.716667 40.383333)",p4s=CRS(wgs.84))
gDistance(MAD,coast) # WGS-84 (long/lat) projection, units in "degrees"
# [1] 3.021808

coast.proj <- spTransform(coast,CRS(epsg.2062))
MAD.proj <- spTransform(MAD,CRS(epsg.2062))
gDistance(MAD.proj,coast.proj) #EPSG-2062 projection, units are in meters.
# [1] 305171.2

So the minimum distance is ~305.2km.

Finally, note that your coastline file has all the coastlines of the world, so this is the minimum distance to some coastline, not necessarily the Spanish coast (although in this case it does turn out to be on the northern coast of Spain). If your reference point was very near the border with Portugal, the the nearest coastal point would be to the western coast of Portugal.

calculating minimum distance between a point and the coast in the UK

AFAICT you are not doing anything wrong. The error message is telling you that some of the coordinates in the global coastline shapefile map to Inf. I'm not sure why this is happening exactly, but generally expecting a planar projection in a specific region to work globally is not a good idea (although it did work in the other question...). One workaround is to clip the coastline shapefile to the bounding box for your specific projection, then transform the clipped shapefile. The bounding box for EPSG.27700 can be found here.

# use bounding box for epsg.27700
# found here: http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
bbx <- readWKT("POLYGON((-7.5600 49.9600, 1.7800 49.9600, 1.7800 60.8400, -7.5600 60.8400, -7.5600 49.9600))",
p4s=CRS(wgs.84))
coast <- gIntersection(coast,bbx) # just coastlines within bbx
# now transformation works
coast.proj <- spTransform(coast,CRS(epsg.27700))
MAD.proj <- spTransform(MAD,CRS(epsg.27700))
dist.27700 <- gDistance(MAD.proj,coast.proj) # distance in sq.meters
dist.27700
# [1] 32153.23

So the closest coastline is 32.2 km away. We can also identify where on the coast this is

# identify the closest point
th <- seq(0,2*pi,len=1000)
circle <- cbind(1.00001*dist.wgs84*cos(th)+MAD$x,1.00001*dist.wgs84*sin(th)+MAD$y)
sp.circle <- SpatialLines(list(Lines(list(Line(circle)),ID="1")),proj4string=CRS(wgs.84))
sp.pts <- gIntersection(sp.circle,coast)
sp.pts
# SpatialPoints:
# x y
# 1 -0.2019854 50.83079
# 1 -0.1997009 50.83064
# Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

And finally we can plot the results. You should always do this.

# plot the results: plot is in CRS(wgs.84)
plot(coast, asp=1)
plot(bbx,add=T)
plot(MAD,add=T, col="red", pch=20)
plot(sp.circle,add=T, col="blue", lty=2)
lines(c(MAD$x,sp.pts$x),c(MAD$y,sp.pts$y),col="red")

Sample Image

R measuring distance from a coastline

Distance to the coastline can be calculated by downloading openstreetmap coastline data. You can then use geosphere::dist2Line to get the distance from your points to the coastline.

I noticed that one of your example points was in France so you may need to expand the coastline data beyond just the UK (can be done by playing with the extents of the bounding box).

library(tidyverse)
library(sf)
library(geosphere)
library(osmdata)

#get initial data frame
d1 <- data_frame(
title = c("base1", "base2", "base3", "base4",
"base5", "base6", "base7"),
lat = c(57.3, 58.8, 47.2, 57.8, 65.4, 56.7, 53.3),
long = c(0.4, 3.4, 3.5, 1.2, 1.5, 2.6, 2.7))

# convert to sf object
d1_sf <- d1 %>% st_as_sf(coords = c('long','lat')) %>%
st_set_crs(4326)

# get bouding box for osm data download (England) and
# download coastline data for this area
osm_box <- getbb (place_name = "England") %>%
opq () %>%
add_osm_feature("natural", "coastline") %>%
osmdata_sf()

# use dist2Line from geosphere - only works for WGS84
#data
dist <- geosphere::dist2Line(p = st_coordinates(d1_sf),
line =
st_coordinates(osm_box$osm_lines)[,1:2])

#combine initial data with distance to coastline
df <- cbind( d1 %>% rename(y=lat,x=long),dist) %>%
mutate(miles=distance/1609)

# title y x distance lon lat miles
#1 base1 57.3 0.4 219066.40 -2.137847 55.91706 136.15065
#2 base2 58.8 3.4 462510.28 -2.137847 55.91706 287.45201
#3 base3 47.2 3.5 351622.34 1.193198 49.96737 218.53470
#4 base4 57.8 1.2 292210.46 -2.137847 55.91706 181.60998
#5 base5 65.4 1.5 1074644.00 -2.143168 55.91830 667.89559
#6 base6 56.7 2.6 287951.93 -1.621963 55.63143 178.96329
#7 base7 53.3 2.7 92480.24 1.651836 52.76027 57.47684

#plot
p <- ggplot() +
geom_sf(data=osm_box$osm_lines) +
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat))

Sample Image

That's just for the distance to the coastline. You also need to know how whether it is inland or at sea. For this, you would need a separate shapefile for the sea: http://openstreetmapdata.com/data/water-polygons and see if each point of your points sits in the sea or not.

#read in osm water polygon data
sea <- read_sf('water_polygons.shp')

#get get water polygons that intersect our points
in_sea <- st_intersects(d1_sf,sea) %>% as.data.frame()

#join back onto original dataset
df %>% mutate(row = row_number()) %>%
#join on in_sea data
left_join(in_sea,by=c('row'='row.id')) %>%
mutate(in_sea = if_else(is.na(col.id),F,T)) %>%
#categorise into 'sea', 'coast' or 'land'
mutate(where = case_when(in_sea == T ~ 'Sea',
in_sea == F & miles <=3 ~ 'Coast',
in_sea == F ~ 'Land'))

# title y x distance lon lat miles row col.id in_sea where
#1 base1 57.3 0.4 219066.40 -2.137847 55.91706 136.15065 1 24193 TRUE Sea
#2 base2 58.8 3.4 462510.28 -2.137847 55.91706 287.45201 2 24194 TRUE Sea
#3 base3 47.2 3.5 351622.34 1.193198 49.96737 218.53470 3 NA FALSE Land
#4 base4 57.8 1.2 292210.46 -2.137847 55.91706 181.60998 4 24193 TRUE Sea
#5 base5 65.4 1.5 1074644.00 -2.143168 55.91830 667.89559 5 25417 TRUE Sea
#6 base6 56.7 2.6 287951.93 -1.621963 55.63143 178.96329 6 24193 TRUE Sea
#7 base7 53.3 2.7 92480.24 1.651836 52.76027 57.47684 7 24143 TRUE Sea

ggplot() +
geom_sf(data=osm_box$osm_lines) +
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat)) +
ggrepel::geom_text_repel(data=df,
aes(x=x,y=y,label=paste0(where,'\n',round(miles,0),'miles')),size=2)

Sample Image

Update 16/08/2018

Since you asked for an approach specifically using a shapefile I have downloaded this one here: openstreetmapdata.com/data/coastlines which I will use to carry out the same approach as above.

clines <- read_sf('lines.shp') #path to shapefile

Next I created a custom bounding box so that we can cut down the size of the shapefile to only include coastlines reasonably close to the points.

# create bounding box surrounding points 
bbox <- st_bbox(d1_sf)

# write a function that takes the bbox around our points
# and expands it by a given amount of metres.
expand_bbox <- function(bbox,metres_x,metres_y){

box_centre <- bbox %>% st_as_sfc() %>%
st_transform(crs = 32630) %>%
st_centroid() %>%
st_transform(crs = 4326) %>%
st_coordinates()


bbox['xmin'] <- bbox['xmin'] - (metres_x / 6370000) * (180 / pi) / cos(bbox['xmin'] * pi/180)
bbox['xmax'] <- bbox['xmax'] + (metres_x / 6370000) * (180 / pi) / cos(bbox['xmax'] * pi/180)
bbox['ymin'] <- bbox['ymin'] - (metres_y / 6370000) * (180 / pi)
bbox['ymax'] <- bbox['ymax'] + (metres_y / 6370000) * (180 / pi)


bbox['xmin'] <- ifelse(bbox['xmin'] < -180, bbox['xmin'] + 360, bbox['xmin'])
bbox['xmax'] <- ifelse(bbox['xmax'] > 180, bbox['xmax'] - 360, bbox['xmax'])
bbox['ymin'] <- ifelse(bbox['ymin'] < -90, (bbox['ymin'] + 180)*-1, bbox['ymin'])
bbox['ymax'] <- ifelse(bbox['ymax'] > 90, (bbox['ymax'] + 180)*-1, bbox['ymax'])
return(bbox)
}

# expand the bounding box around our points by 300 miles in x and 100 #miles in y direction to make nice shaped box.
bbox <- expand_bbox(bbox,metres_x=1609*200, metres_y=1609*200) %>% st_as_sfc

# get only the parts of the coastline that are within our bounding box
clines2 <- st_intersection(clines,bbox)

Now I used the dist2Line function here because it is accurate and it gives you the points on the coastline it is measuring to which makes it good for checking errors. The downside is, it's very slow for our rather large coastline file.

Running this took me 8 minutes:

dist <- geosphere::dist2Line(p = st_coordinates(d1_sf), 
line = as(clines2,'Spatial'))

#combine initial data with distance to coastline
df <- cbind( d1 %>% rename(y=lat,x=long),dist) %>%
mutate(miles=distance/1609)

df

# title y x distance lon lat ID miles
#1 base1 57.3 0.4 131936.70 -1.7711149 57.46995 4585 81.99919
#2 base2 58.8 3.4 98886.42 4.8461433 59.28235 179 61.45831
#3 base3 47.2 3.5 340563.02 0.3641618 49.43811 4199 211.66129
#4 base4 57.8 1.2 180110.10 -1.7670712 57.50691 4584 111.93915
#5 base5 65.4 1.5 369550.43 6.2494627 62.81381 9424 229.67709
#6 base6 56.7 2.6 274230.37 5.8635346 58.42913 24152 170.43528
#7 base7 53.3 2.7 92480.24 1.6518358 52.76027 4639 57.47684

plot:

ggplot() + 
geom_sf(data=clines2) +
geom_sf(data=bbox,fill=NA)+
geom_sf(data=d1_sf) +
geom_segment(data=df,aes(x=x,y=y,xend=lon,yend=lat))

Sample Image

If you don't mind about the slight loss of accuracy (results differ by around 0.3% on your data), and are not fussed about knowing where exactly on the coastline it is measuring to, you can measure the distance to the polygon:

# make data into polygons
clines3 <- st_intersection(clines,bbox) %>%
st_cast('POLYGON')

#use rgeos::gDistance to calculate distance to nearest polygon
#need to change projection (I used UTM30N) to use gDistance
dist2 <- apply(rgeos::gDistance(as(st_transform(d1_sf,32630), 'Spatial'),
as(st_transform(clines3,32630),'Spatial'),
byid=TRUE),2,min)

df2 <- cbind( d1 %>% rename(y=lat,x=long),dist2) %>%
mutate(miles=dist2/1609)

df2

# title y x dist2 miles
#1 base1 57.3 0.4 131917.62 81.98733
#2 base2 58.8 3.4 99049.22 61.55949
#3 base3 47.2 3.5 341015.26 211.94236
#4 base4 57.8 1.2 180101.47 111.93379
#5 base5 65.4 1.5 369950.32 229.92562
#6 base6 56.7 2.6 274750.17 170.75834
#7 base7 53.3 2.7 92580.16 57.53894

By contrast this took just 8 seconds to run!

The rest is as in the previous answer.

Calculate minimum distance from polygon to spatial point

There are two issues with your current approach.
1.) The coordinates assignment isn't quite correct. It should be long/lat, but you're assigning it as lat/long.
2.) Directly setting the CRS in the current way will not actually alter the points in the necessary way. You need to assign an appropriate long/lat CRS first, and then perform a spTransform action.

#Open data on tin mines and define as spatial points
tin.01 <- read.csv("random/Tin Mine Location_01.csv")
coordinates(tin.01) <- c("Longitude","Latitude") ## Should be `9:8` if you wanted to stick with indexes, but using the names here is generally lower risk.
proj4string(tin.01) <- CRS(proj4string(Kinta)) ## Setting the initial projection to the same one the polygons are using. You should change this if your original data source uses some other known long/lat projection.

tin.km <- spTransform(tin.01,CRS(EPSG.3375)) ## Creating a transformed set of points for the distance calculation.

#Find distance between district and mines
gDistance(Kinta.km,tin.km,byid=TRUE) ## '0' distance means the mine is inside the district.

59
1 194384.372
2 223773.999
3 0.000
4 36649.914
5 102944.361
6 0.000
7 0.000
8 6246.066
9 0.000

Calculate distance between 2 lon lats but avoid going through a coastline in R

Not fully errorchecked yet but it may get you started. Rather than coastlines, I think you need to start with a raster whose the no-go areas are set to NA.

library(raster)
library(gdistance)
library(maptools)
library(rgdal)

# a mockup of the original features dataset (no longer available)
# as I recall it, these were just a two-column matrix of xy coordinates
# along the coast of East Antarctica, in degrees of lat/long
ant_x <- c(85, 90, 95, 100)
ant_y <- c(-68, -68, -68, -68)
feature <- cbind(ant_x, ant_y)

# a projection I found for antarctica
antcrs <- crs("+proj=stere +lat_0=-90 +lat_ts=-71 +datum=WGS84")

# set projection for your features
# function 'project' is from the rgdal package
antfeat <- project(feature, crs(antcrs, asText=TRUE))

# make a raster similar to yours
# with all land having "NA" value
# use your own shapefile or raster if you have it
# the wrld_simpl data set is from maptools package
data(wrld_simpl)
world <- wrld_simpl
ant <- world[world$LAT < -60, ]
antshp <- spTransform(ant, antcrs)
ras <- raster(nrow=300, ncol=300)
crs(ras) <- crs(antshp)
extent(ras) <- extent(antshp)
# rasterize will set ocean to NA so we just inverse it
# and set water to "1"
# land is equal to zero because it is "NOT" NA
antmask <- rasterize(antshp, ras)
antras <- is.na(antmask)

# originally I sent land to "NA"
# but that seemed to make some of your features not visible
# so at 999 land (ie everything that was zero)
# becomes very expensive to cross but not "impossible"
antras[antras==0] <- 999
# each cell antras now has value of zero or 999, nothing else

# create a Transition object from the raster
# this calculation took a bit of time
tr <- transition(antras, function(x) 1/mean(x), 8)
tr = geoCorrection(tr, scl=FALSE)

# distance matrix excluding the land
# just pick a few features to prove it works
sel_feat <- head(antfeat, 3)
A <- accCost(tr, sel_feat)

# now A still shows the expensive travel over land
# so we mask it out for sea travel only
A <- mask(A, antmask, inverse=TRUE)
plot(A)
points(sel_feat)

Seems to be working because the left side ocean has higher values than the right side ocean, and likewise as you go down into the Ross Sea.

Sample Image

Spatial data: calculating the distance of points from the maximum point value and plotting

I'm still not sure I understand your question but I propose the following answer.

Load packages

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(tidyr)

Load data

dat = structure(
list(
date.hour = structure(
c(
1551057840, 1551057840, 1551057840, 1551057840, 1551057840,
1551057840, 1551057840
),
tzone = "UTC",
class = c(
"POSIXct",
"POSIXt"
)
),
id = c(2, 5, 7, 8, 9, 10, 11),
variable = c(
456, 27, 130, 116, 92, 141, 145
),
xy_coord = c(
"6.2 14.8", "8.2 8.9", "4.2 8.9", "2.2 8.9", "8.2 3.5", "6.2 3.5",
"4.2 3.5"
)
),
row.names = c(NA,-7L),
groups = structure(
list(
id = c(2, 5, 7, 8, 9, 10, 11),
date.hour = structure(
c(
1551057840, 1551057840, 1551057840, 1551057840, 1551057840,
1551057840, 1551057840
),
tzone = "UTC",
class = c(
"POSIXct",
"POSIXt"
)
),
.rows = structure(
list(1L, 2L, 3L, 4L, 5L, 6L, 7L),
ptype = integer(0),
class = c(
"vctrs_list_of", "vctrs_vctr", "list"
)
)
),
row.names = c(NA, -7L),
class = c("tbl_df", "tbl", "data.frame"),
.drop = TRUE
),
class = c("grouped_df", "tbl_df", "tbl", "data.frame")
)

Separate the xy_coord column, convert columns to numeric and create an sf object

dat_sf <- st_as_sf(
separate(dat, xy_coord, c("x", "y"), sep = " ", convert = TRUE),
coords = c("x", "y")
)

Find the maximum of variable

which.max(dat_sf[["variable"]])
#> [1] 1

Compute all distances

dat_sf[["distances"]] <- st_distance(dat_sf, dat_sf[1, ])

Plot

plot(variable ~ distances, data = dat_sf)

Sample Image

Created on 2021-11-22 by the reprex package (v2.0.1)

You can also remove the first point (with distance = 0).

Shortest distance between a point and a line (Google Maps API issue?)

As far as I can see, your formula for D is correct, and the linear approximation is justified at such a small scale (deltas about a quarter of a degree; relative errors due to non-linearity should be on the order of 10^-5).

What you see can be due to the fact that the map projection is not conformant (does not preserve angles), so that the angle is not displayed as right. But the point is correct.

Do you know which projection they use ?


Bingo, the angle is right, just a display artifact due to the projection.

Sample Image

Find distance from a point to a polygon

In the case that the polygon indeed describes a circle, you could save the polygon as a center location (x,y) coordinate and the radius of the circle. The distane of a point to the polygon can be computed as the distance from the center of the circle to the wanted point, and then reduce the raidus size. As a bonus, if the resulted distance is negative, your point is inside the circle.



Related Topics



Leave a reply



Submit