使用R通过WFS获取3D BAG数据时突破服务器5000记录上限求助

问题背景 #

我在开展物种生态位建模工作，需要从代尔夫特理工大学（TU Delft）的3D BAG数据中提取哈勒姆市的建筑高度。手动下载数据切片效率低且易遗漏，还要批量处理多个城市，因此尝试通过WFS服务获取要素。使用扩展1.2倍的边界框请求时，服务器最多返回5000条记录，再加上对WFS语义（命名空间、要素类型、属性）不熟悉，试了多种方法都没突破限制。

已尝试的方法 #

• 分页：教程要求通过设置resultType="hits"获取总要素数才能分页，但无法轻松拿到目标边界框内的总要素数
• 按切片ID过滤：计划利用BAG3D_v2:bag_tiles_3k图层的tile_id属性，先提取匹配边界框的切片ID再逐片获取要素，但连单个切片的CQL过滤条件都没写成功
• 拆分边界框：考虑用R包slippymath把大边界框拆成小切片后逐片请求，但过滤问题仍未解决

基础代码 #

library(httr)

url <- parse_url("https://data.3dbag.nl/api/BAG3D_v2/wfs")
url$query <- list(service = "WFS",
                  version = "2.0.0",
                  request = "GetFeature",
                  typename = "BAG3D_v2:lod22",
                  #cql_filter = "BAG3D_v2:tile_id ='4199'",
                  bbox = "100768.4,482708.5,107923.1,494670.4",
                  startindex = 10000,
                  sortBy = "gid")
request <- build_url(url)
test <- st_read(request)
qtm(test)

解决方案步骤 #

1. 正确获取总要素数（解决分页问题） #

先发送带resultType="hits"的请求，解析响应拿到边界框内的总要素数，再计算分页次数循环请求：

library(xml2)

url_hits <- parse_url("https://data.3dbag.nl/api/BAG3D_v2/wfs")
url_hits$query <- list(service = "WFS",
                       version = "2.0.0",
                       request = "GetFeature",
                       typename = "BAG3D_v2:lod22",
                       bbox = "100768.4,482708.5,107923.1,494670.4",
                       resultType = "hits")
request_hits <- build_url(url_hits)

# 解析响应获取总记录数
hits_response <- GET(request_hits)
total_features <- xml_find_first(read_html(hits_response), "//*[@numberOfFeatures]") %>% 
  xml_attr("numberOfFeatures") %>% 
  as.integer()

# 分页循环获取数据
all_data <- list()
page_size <- 5000
pages <- ceiling(total_features / page_size)

for(page in 0:(pages-1)){
  url_page <- parse_url("https://data.3dbag.nl/api/BAG3D_v2/wfs")
  url_page$query <- list(service = "WFS",
                         version = "2.0.0",
                         request = "GetFeature",
                         typename = "BAG3D_v2:lod22",
                         bbox = "100768.4,482708.5,107923.1,494670.4",
                         startindex = page * page_size,
                         count = page_size,
                         sortBy = "gid")
  request_page <- build_url(url_page)
  page_data <- st_read(request_page)
  all_data[[page+1]] <- page_data
}

final_data <- do.call(rbind, all_data)

2. 正确编写CQL过滤条件（按切片ID获取） #

先获取匹配边界框的切片ID列表，再循环每个ID发送过滤请求：

# 获取边界框内的tile_id列表
url_tiles <- parse_url("https://data.3dbag.nl/api/BAG3D_v2/wfs")
url_tiles$query <- list(service = "WFS",
                        version = "2.0.0",
                        request = "GetFeature",
                        typename = "BAG3D_v2:bag_tiles_3k",
                        bbox = "100768.4,482708.5,107923.1,494670.4",
                        outputFormat = "csv")
request_tiles <- build_url(url_tiles)
tiles_df <- read.csv(request_tiles)
tile_ids <- tiles_df$tile_id

# 循环每个tile_id获取数据
all_tile_data <- list()
for(tile in tile_ids){
  url_tile <- parse_url("https://data.3dbag.nl/api/BAG3D_v2/wfs")
  url_tile$query <- list(service = "WFS",
                         version = "2.0.0",
                         request = "GetFeature",
                         typename = "BAG3D_v2:lod22",
                         cql_filter = paste0("tile_id = '", tile, "'"),
                         outputFormat = "GML3")
  request_tile <- build_url(url_tile)
  tile_data <- st_read(request_tile)
  all_tile_data[[as.character(tile)]] <- tile_data
}

final_tile_data <- do.call(rbind, all_tile_data)

若服务器要求带命名空间，将过滤条件改为paste0("BAG3D_v2:tile_id = '", tile, "'")即可。

3. 拆分边界框的简化处理 #

用slippymath将大边界框拆分为多个子区域，循环每个子区域发送请求：

library(slippymath)

# 转换边界框为sf格式（坐标为EPSG:28992）
bbox_sf <- st_bbox(c(xmin=100768.4, ymin=482708.5, xmax=107923.1, ymax=494670.4), crs=st_crs(28992))
# 拆分为9个小区域（可根据实际调整数量）
grid <- make_grid(bbox_sf, n = c(3,3))

# 循环每个网格获取数据
all_grid_data <- list()
for(g in seq_along(grid)){
  grid_bbox <- st_bbox(grid[g])
  grid_bbox_str <- paste(grid_bbox$xmin, grid_bbox$ymin, grid_bbox$xmax, grid_bbox$ymax, sep=",")
  
  url_grid <- parse_url("https://data.3dbag.nl/api/BAG3D_v2/wfs")
  url_grid$query <- list(service = "WFS",
                         version = "2.0.0",
                         request = "GetFeature",
                         typename = "BAG3D_v2:lod22",
                         bbox = grid_bbox_str,
                         outputFormat = "GML3")
  request_grid <- build_url(url_grid)
  grid_data <- st_read(request_grid)
  all_grid_data[[g]] <- grid_data
}

final_grid_data <- do.call(rbind, all_grid_data)

关键注意事项 #

• 优先使用outputFormat="application/json"（GeoJSON）或csv格式，比GML更易处理
• 可访问WFS的GetCapabilities接口，确认要素类型、属性的完整命名空间和参数规则
• 若遇到请求超时，可在循环中添加Sys.sleep(1)避免频繁请求触发服务器限制

内容的提问来源于stack exchange，提问作者Berry