Calculate distance between 2 lon lats but avoid go

2019-07-21 07:10发布

I am trying to calculate the closest distance between locations in the ocean and points on land but not going through a coastline. Ultimately, I want to create a distance to land-features map.

This map was created using rdist.earth and is a straight line distance. Therefore it is not always correct because it not taking into account the curvatures of the coastline.

c<-matrix(coast_lonlat[,1], 332, 316, byrow=T)
image(1:316, 1:332, t(c))

min_dist2_feature<-NULL
for(q in 1:nrow(coast_lonlat)){
diff_lonlat <- rdist.earth(matrix(coast_lonlat[q,2:3],1,2),as.matrix(feature[,1:2]), miles = F)
min_dist2_feature<-c(min_dist2_feature, min(diff_lonlat,na.rm=T))
}

distmat <- matrix( min_dist2_feature, 316, 332)
image(1:316, 1:332, distmat)

Land feature data is a two column matrix of xy coordinates, e.g.:

ant_x <- c(85, 90, 95, 100)
ant_y <- c(-68, -68, -68, -68)
feature <- cbind(ant_x, ant_y)

Does anyone have any suggestions? Thanks

1条回答
太酷不给撩
2楼-- · 2019-07-21 07:41

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.

enter image description here

查看更多
登录 后发表回答