我在尝试使用 raster::getData() 下载 SRTM 数据时遇到了完全相同的错误。正如samAct 所提到的，问题与不再可访问的 URL 有关。
我的解决方案是编辑 getData 函数，使其首先尝试http://srtm.csi.cgiar.org/SRT-ZIP/SRTM_V41/SRTM_Data_GeoTiff/FILENAME URL。
要使用我编辑的 raster 包版本，只需使用 `devtools::install_github("pokyah/raster")。这将覆盖您当前的版本并使 getData 使用有效的 URL。
devtools::install_github("pokyah/raster")
library(raster)
然后，您可以使用此功能下载您感兴趣区域的 SRTM 数据：
# @param country_code.chr a character specifying the ISO contrycode. Ex : BE for belgium
# @param NAME_1.chr a character specifying the NAME_1 value for lower than country level information
# @param aggregation_factor.num a numeric specifying the aggregation factor to get the desired spatial resolution
# @param EPSG.chr a character specifying the EPSG code of the desired Coordiante Reference System (CRS)
# @param path.chr a character specifying the path where to dowload the SRTM data
build_highRes_terrain_rasters.fun <- function(country_code.chr, NAME_1.chr=NULL, aggregation_factor.num=NULL, EPSG.chr=NULL, path.chr) {
# Path to downloaded SRTM Tiles refs
srtm.tiles.ref <- raster::shapefile("<PATH_TO_DOWNLOADED_TILES_REFS")
# Get country geometry first
extent.sp <- raster::getData('GADM', country=country_code.chr, level=1)
if(!is.null(NAME_1.chr)){
extent.sp <- subset(extent.sp, NAME_1 == NAME_1.chr)
}
# Intersect extent geometry and tile grid
intersects <- rgeos::gIntersects(extent.sp, srtm.tiles.ref, byid=T)
tiles <- srtm.tiles.ref[intersects[,1],]
# Download tiles using getData
# inspired from https://www.gis-blog.com/download-srtm-for-an-entire-country/
srtm_list <- list()
for(i in 1:length(tiles)) {
lon <- extent(tiles[i,])[1] + (extent(tiles[i,])[2] - extent(tiles[i,])[1]) / 2
lat <- extent(tiles[i,])[3] + (extent(tiles[i,])[4] - extent(tiles[i,])[3]) / 2
tile <- getData('SRTM', #data are downloaded from http://www.cgiar-csi.org/. See getData do of pokyah/raster repo on github
lon=lon,
lat=lat,
download = FALSE,
path = path.chr)
srtm_list[[i]] <- tile
}
# Mosaic tiles
srtm_list$fun <- mean
srtm_mosaic.ras <- do.call(raster::mosaic, srtm_list)
# Crop tiles to extent borders
extent.elevation.ras <- raster::crop(srtm_mosaic.ras, extent.sp)
extent.elevation.ras <- raster::mask(extent.elevation.ras, extent.sp)
# transform to desired CRS
if(!is.null(EPSG.chr)){
raster::projectRaster(extent.elevation.ras, crs = toString((dplyr::filter(rgdal::make_EPSG(), code==EPSG.chr))$prj4))
}
# aggregate to lower resolution
# inspired from https://stackoverflow.com/questions/32278825/how-to-change-the-resolution-of-a-raster-layer-in-r
if(!is.null(aggregation_factor.num)){
extent.elevation.ras <- raster::aggregate(extent.elevation.ras, fact=aggregation_factor.num)
}
# compute the slope from the elevation
# inspired from https://rpubs.com/etiennebr/visualraster
extent.slope.ras <- raster::terrain(extent.elevation.ras, opt="slope", unit="degrees")
extent.aspect.ras <- raster::terrain(extent.elevation.ras, opt="aspect", unit="degrees")
extent.roughness.ras <- raster::terrain(extent.elevation.ras, opt="roughness")
# compute the aspect from the elevation
extent.terrain.ras = raster::stack(
extent.elevation.ras,
extent.slope.ras,
extent.aspect.ras,
extent.roughness.ras)
}
希望对您有所帮助！