mesh builder

 library(fmesher)

library(sf)

library(ggplot2)

create_adaptive_mesh <- function(sf_object,

                                 inner_edge_ratio = 30,    # 增大内部边长比例（更密集）

                                 outer_edge_ratio = 5,     # 减小外部边长比例（更稀疏）

                                 inner_offset_ratio = 8,   # 减小内部偏移比例（内部区域更大）

                                 outer_offset_ratio = 20,  # 增大外部偏移比例（外部区域相对更小）

                                 cutoff_ratio = 80,        # 增大cutoff比例（避免过小三角形）

                                 crs_transform = NULL) {   # 可选的坐标转换

  # 坐标系转换（如果需要）

  if (!is.null(crs_transform)) {

    sf_object <- st_transform(sf_object, crs_transform)

  }

  # 创建边界

  boundary <- fm_as_segm(sf_object)

  # 计算自适应参数

  coords <- st_coordinates(sf_object)

  lon_range <- diff(range(coords[, 1]))

  lat_range <- diff(range(coords[, 2]))

  mean_range <- mean(c(lon_range, lat_range))

  # 关键修改：调整max.edge参数

  # 第一个参数控制边界内部，第二个参数控制边界外部

  inner_max_edge <- mean_range/inner_edge_ratio  # 内部：更小的边长 = 更密集

  outer_max_edge <- mean_range/outer_edge_ratio  # 外部：更大的边长 = 更稀疏

  # 创建mesh

  mesh <- fm_mesh_2d(

    boundary = boundary,

    max.edge = c(inner_max_edge, outer_max_edge),  # 内部小，外部大

    offset = c(mean_range/inner_offset_ratio, mean_range/outer_offset_ratio),

    cutoff = mean_range/cutoff_ratio

  )

  # 打印网格信息

  cat("网格参数信息:\n")

  cat("- 内部最大边长:", round(inner_max_edge, 4), "\n")

  cat("- 外部最大边长:", round(outer_max_edge, 4), "\n")

  cat("- 密度比例 (外部/内部):", round(outer_max_edge/inner_max_edge, 2), "\n")

  cat("- 顶点总数:", mesh$n, "\n")

  cat("- 三角形总数:", nrow(mesh$graph$tv), "\n")

  # 返回mesh和元信息

  list(

    mesh = mesh,

    boundary = boundary,

    range_info = c(lon_range = lon_range, lat_range = lat_range, mean_range = mean_range),

    vertices_count = mesh$n,

    triangles_count = nrow(mesh$graph$tv),

    edge_ratio = outer_max_edge/inner_max_edge

  )

}

# 可视化函数

plot_mesh_density <- function(mesh_result) {

  mesh <- mesh_result$mesh

  # 转换为数据框用于ggplot

  vertices_df <- data.frame(

    x = mesh$loc[,1],

    y = mesh$loc[,2]

  )

  triangles_df <- data.frame(

    triangle_id = 1:nrow(mesh$graph$tv),

    mesh$graph$tv

  )

  ggplot() +

    geom_point(data = vertices_df, aes(x = x, y = y),

               size = 0.3, alpha = 0.6, color = "blue") +

    theme_minimal() +

    coord_equal() +

    labs(title = "自适应网格分布图",

         subtitle = paste("内部密集，外部稀疏 | 顶点数:", mesh_result$vertices_count),

         x = "经度", y = "纬度") +

    theme(plot.title = element_text(hjust = 0.5),

          plot.subtitle = element_text(hjust = 0.5))

}

# 使用示例

gs <- st_read('data/省.geojson') %>%

  filter(省=='甘肃省') %>%

  st_set_crs(4326)

# 创建mesh - 参数已经按照内密外疏原则设置

mesh_result <- create_adaptive_mesh(

  gs,

  inner_edge_ratio = 35,   # 内部更密集（增大比例）

  outer_edge_ratio = 4,    # 外部更稀疏（减小比例）

  inner_offset_ratio = 6,  # 内部区域更大

  outer_offset_ratio = 15, # 外部区域适中

  cutoff_ratio = 100       # 避免过细小的三角形

)

mesh <- mesh_result$mesh

plot(mesh_result$mesh)