在R中，我创建了一个覆盖世界（不包括南极洲）的栅格，大约60，000个点，并且只保留距离其最近的4个相邻点大约100公里的网格中的点。我进一步减少了陆地上的点的数量，得到了这些点的海拔，然后我减少了海拔为0米的点的数量<2000m and >。最后，我有大约14，000个点。

library(raster) ##load package to create raster of points over entire world
library(sp)  ##load package for st_points, st_transform, st_buffer, and st_intersection functions 
library(rgdal) ##elevatr package dependency
library(elevatr) ##load package for get_elev_point function 

r <- raster(ncol = 360, nrow = 180, ymn = -55, ymx = 90) ## Remove Antarctica
points <- as.data.frame(rasterToPoints(r))
points.sf <- st_as_sf(points, coords = 1:2, crs = 4326) ## Convert points to sf object

## Create a bounding box for your points to set the extent for the grid
bbox <- st_bbox(points.sf) ##cannot use r, because will lose precision of Long/Lat after the decimal

## Define the desired spacing between grid points: ~ 100km
#You can convert this distance to degrees using the approximate conversion factor of 1 degree ≈ 111.32 km at the equator. 
grid_spacing_degrees <- 1

## Create a grid of points within the bounding box
x_values <- seq(from = bbox["xmin"], to = bbox["xmax"], by = grid_spacing_degrees)
y_values <- seq(from = bbox["ymax"], to = bbox["ymin"], by = -grid_spacing_degrees)

grid <- expand.grid(x = x_values, y = y_values)

## Convert the grid to an sf object:
grid.sf <- st_as_sf(grid, coords = 1:2, crs = 4326)

## Keep only the points on land
intersection <- st_intersection(grid.sf, world)
##took 2.5 hours

land.grid.points <- intersection

## Check that it worked
plot(st_geometry(world))
plot(land.grid.points, pch=20, cex = 0.00005, max.plot = 10)

## For points on land, get elevation
elevation.land.grid.points <- get_elev_point(land.grid.points, prj = "EPSG:4326", z = 6, src = "aws") #zoom 6 is ground resolution per zoom in km at a given latitude: 0° = 2446m; 45° = 1729.6m; 60° = 1223m; uses Amazon Web Services; elevation in meters
        
elev <- as.data.frame(elevation.land.grid.points)

elev.constrained <- elev[elev$elevation > 0 & elev$elevation < 2000, ] ## fewer than 30 points have elevation = 0m, so removing them as done here is not harmful to the overall number of points
##confirmed elevations < 0 were off the coast a bit and in water, so correct to remove elevations < 0

字符串
接下来，我计划将光栅值设置为高程，并将连接两个顶点的边的长度定义为它们的大圆距离，这样，我将在图中连接具有大圆距离的顶点的边作为边长。

library(geosphere) ##load package for distVincentyEllipsoid function

# Create vertices with latitude, longitude, and elevation information
vertices <- data.frame(
  Lon = elev.constrained$coords.x1,
  Lat = elev.constrained$coords.x2,
  Elevation = elev.constrained$elevation
)
# Calculate great circle distances between vertices
distances <- distVincentyEllipsoid(vertices[, c("Lon", "Lat")], a = 6378137, b = 6356752.3142, f = 1/298.257223563) #meters

# Add edges to the graph using the full distance matrix
edge_weights <- as.vector(distances)
g <- add_edges(g, edges = combn(vertex_indices, 2), edge.attr = edge_weights)

型
实际上，我在这里停了下来，因为创建对象g使用了太多的内存，我正在寻找一种方法。
但是，一旦我得到g，那么我可以继续将栅格值设置为高程，并限制从一个顶点到8个最近邻居的移动。这样，旅行的约束将同时包含高程和距离。

# Load your existing elevation raster
# use rasterize to create desired raster
r_data <- rasterize(x=elev.constrained[, 3:4], # lon-lat data
                    y=r, # raster object
                    field=elev.constrained[, 1], # vals to fill raster with
                    fun=mean) # aggregate function

# Create a new raster to store edge lengths (great circle distances)
edge_length_raster <- r_data

# Assign edge lengths to raster cells
values(edge_length_raster) <- distances

# You now have an elevation raster with edge lengths based on great circle distances

# plot(edge_length_raster) ##checkpoint
tr1 <- transition(edge_length_raster, transitionFunction = mean, directions = 8) ##uses Moore neighborhood

image(transitionMatrix(tr1))

h8 <- geoCorrection(tr1,  scl=FALSE)

型
最后，我将使用Dijkstra算法找到从一个位置到另一个位置的最短路径。
我的问题是，有时候我需要通过水上穿越大陆，比如从俄罗斯的阿纳迪尔到北美。我该如何编码？我是进入我的基本网格并为中间的水指定一个点，还是以某种方式将俄罗斯的阿纳迪尔与北美白令海峡另一边的一个点“连接”起来？