如何用for循环下载GeoTiff并提取地理ID与风险预测数据？

我正在尝试下载并解析GeoTIFF文件，从中提取字母数字数据。这些文件是加州各次级地理区域的热力风险分类预测数据（如“高”“极高”等），我需要提取的是与各预测关联的地理ID（如县名）和预测结果，而非栅格数据（目前暂不需要）。

注：经与NWS沟通确认，数据页面提供的.kml文件为空。

下载.tif文件 #

以下代码可正常运行：

#Load Packages
#install.packages("raster")
#install.packages("terra")
#install.packages("vapour")

library(raster)
library(httr)
library(terra)


#Read In GeoTiff Files

#Base URL
base_url <- "https://www.wpc.ncep.noaa.gov/heatrisk/data/"

#Function to Read in Files
download_heatrisk_geotiff <- function(day_number) {
  url <- paste0(base_url, "?day=", day_number)
  response <- GET(url)
  
  if (status_code(response) == 200) {
    # Parse the TIFF content
    tif_content <- content(response, "text")
    
    # Save KML content to a file
    writeLines(tif_content, paste0(getwd(), "/data/HeatRisk_", day_number, "_Mercator.tif"))
    cat(paste0("Downloaded HeatRisk_", day_number, "_Mercator.tif\n"))
  } else {
    cat(paste0("Error downloading HeatRisk_", day_number, "_Mercator.tif\n"))
  }
}

# Download .tif files for individual days
for (day in 1:7) {
  download_heatrisk_geotiff(day)
}

已下载文件截图 #

不确定所需数据是否在已下载的文件中：

提取数据尝试 #

#Extract Data
geotiff_file <- raster(paste0(getwd(), "/data/HeatRisk_", day_number, "_Mercator.tif"))
geotiff_file2 <- terra::readValues(paste0(getwd(), "/data/HeatRisk_", day_number, "_Mercator.tif"))
geotiff_file3 <- vapour::vapour_read_raster(paste0(getwd(), "/data/HeatRisk_", day_number, "_Mercator.tif"))

已查阅的资料：

• 空间栅格数据标准教程
• R从Web URL直接读取GeoTIFF原始内容
• R语言中读取GeoTIFF并返回特定数据类型

解决方案 #

1. 修正下载代码的核心问题 #

当前下载代码存在致命错误：content(response, "text")会将二进制TIFF文件以文本格式读取，导致文件损坏。正确做法是获取二进制内容并写入：

download_heatrisk_geotiff <- function(day_number) {
  url <- paste0(base_url, "?day=", day_number)
  response <- GET(url)
  
  if (status_code(response) == 200) {
    # 获取二进制TIFF内容
    tif_content <- content(response, "raw")
    # 确保data目录存在
    if (!dir.exists("data")) dir.create("data")
    # 写入二进制文件
    writeBin(tif_content, paste0("data/HeatRisk_", day_number, "_Mercator.tif"))
    cat(paste0("Downloaded HeatRisk_", day_number, "_Mercator.tif\n"))
  } else {
    cat(paste0("Error downloading HeatRisk_", day_number, "_Mercator.tif\n"))
  }
}

2. 提取地理ID与风险分类 #

这些GeoTIFF是分类栅格，存储的是风险编码值，需关联地理边界提取对应区域的风险结果：

步骤1：加载并查看栅格属性

用terra加载修复后的文件，查看编码对应的风险标签：

library(terra)
# 加载某一天的栅格
r <- rast("data/HeatRisk_1_Mercator.tif")
# 查看分类编码与对应标签
cat("分类编码与风险标签:\n")
print(cats(r)[[1]])

步骤2：获取加州地理边界数据

用tigris包获取加州县边界（若需要更细的次级区域，需替换对应矢量数据）：

#install.packages("tigris")
library(tigris)
# 获取加州县边界
ca_counties <- counties(state = "CA", cb = TRUE)
# 转换为与栅格一致的投影
ca_counties <- project(ca_counties, crs(r))

步骤3：提取每个区域的风险值

提取每个县的主导风险类别（出现次数最多的编码），并关联标签：

# 提取每个县的栅格值，取众数作为主导类别
extracted <- extract(r, ca_counties, fun = modal, na.rm = TRUE)
# 合并县名与风险编码
result <- cbind(ca_counties$NAME, extracted)
colnames(result) <- c("County", "Risk_Code")
# 关联风险标签
risk_labels <- cats(r)[[1]]
final_result <- merge(result, risk_labels, by.x = "Risk_Code", by.y = "ID")
# 查看最终结果
print(final_result)

3. 补充说明 #

• 若栅格存在无数据区域（NA），需在extract中通过na.rm = TRUE过滤；
• 若需要次级地理区域（非县）的数据，需获取对应级别的官方矢量边界文件；
• 分类栅格的属性表中已包含编码对应的风险描述（如1=低、2=中等、3=高、4=极高），可直接匹配使用。

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