计算 2 lon lats 之间的距离,但避免穿过 R 中的海岸线
Calculate distance between 2 lon lats but avoid going through a coastline in R
我正在尝试计算海洋中的位置与陆地上的点之间的最近距离,但不经过海岸线。最终,我想创建一个到地物距离的地图。
这张地图是使用rdist.earth创建的,是一条直线距离。因此它并不总是正确的,因为它没有考虑到海岸线的曲率。
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)
地物数据是xy坐标的两列矩阵,例如:
ant_x <- c(85, 90, 95, 100)
ant_y <- c(-68, -68, -68, -68)
feature <- cbind(ant_x, ant_y)
有人有什么建议吗?谢谢
尚未完全检查错误,但它可以帮助您入门。我认为您需要从禁止区域设置为 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)
似乎有效,因为左侧海洋的值高于右侧海洋,当您进入罗斯海时也是如此。
我正在尝试计算海洋中的位置与陆地上的点之间的最近距离,但不经过海岸线。最终,我想创建一个到地物距离的地图。
这张地图是使用rdist.earth创建的,是一条直线距离。因此它并不总是正确的,因为它没有考虑到海岸线的曲率。
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)
地物数据是xy坐标的两列矩阵,例如:
ant_x <- c(85, 90, 95, 100)
ant_y <- c(-68, -68, -68, -68)
feature <- cbind(ant_x, ant_y)
有人有什么建议吗?谢谢
尚未完全检查错误,但它可以帮助您入门。我认为您需要从禁止区域设置为 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)
似乎有效,因为左侧海洋的值高于右侧海洋,当您进入罗斯海时也是如此。