2017年新年快乐!
It’1月1日,我制定了新年决议,今年将发布更多。 ðŸ™,因此,这是有关如何下载SRTM的简短教程 (S 航天飞机雷达地形任务) 90m resolution 数据 (3 arc seconds) for an entire country using the {raster}
and {rgeos}
R packages.
基础知识(下载单个SRTM图块)
如一 以前的帖子 it’s possible to 下载 the 90m SRTM using the getData()
function as follows:
srtm <- getData('SRTM', lon=16, lat=48)
- 选择数据集: 第一个参数指定数据集。‘SRTM’返回SRTM 90高程数据。
- 指定Lon: 第二个参数指定SRTM磁贴的区域。
- 指定纬度:第二个参数指定SRTM磁贴的纬度。
上面的代码将在维也纳附近的某个地方返回一个原始的SRTM Tile:
进阶下载(整个国家/地区)
指定数据下载的经度和纬度通常不是很方便。特别是当我们想下载整个国家的SRTM数据时。通常,一个国家/地区包含多个图块,并且查找每个图块的坐标非常麻烦。
这是一个简短的脚本,您只需要指定ISO国家/地区代码,就会自动拼接和绘制图块。您可以在以下位置找到代码和所需的所有数据 的github.
制备: 在开始之前,请下载包含SRTM磁贴网格的shapefile。可以找到 这里 或上 的github。没有这个文件,下面的代码赢了’t work.
下载shapefile之后,我们就可以启动脚本了。首先加载必要的库,这些库是:
library(raster) library(rgeos) library(rasterVis)
然后为您的SRTM数据下载指定国家/地区:
#Specify target ISO country code and path to 下载ed shapefile country_name <- "AUT" #Austria shp <- shapefile("srtm/tiles.shp") #Path to 下载ed SRTM Tiles
在完成这些初始步骤之后,您可以从下面运行代码并向后倾斜。该脚本会自动确定下载所需的所有切片,镶嵌图并将它们裁剪为一个文件:
#Get country geometry first country <- getData("GADM",                   country = country_name,                   level=0) #Intersect country geometry and tile grid intersects <- gIntersects(country, shp, byid=T) tiles     <- shp[intersects[,1],] #Download tiles 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',                  lon=lon,                  lat=lat)   srtm_list[[i]] <- tile } #Mosaic tiles srtm_list$fun <- mean srtm_mosaic  <- do.call(mosaic, srtm_list) #Crop tiles to country borders srtm_crop    <- mask(srtm_mosaic, country) #Plot p <- levelplot(srtm_mosaic) p + layer(sp.lines(country,                   lwd=0.8,                   col='darkgray'))
在这里,奥地利的结果用 {rasterVis}:
如果您有任何疑问,请在下面发表评论。
干杯!
马丁
10条留言
您可以在这篇文章中发表评论。
嗨,马丁,
很棒的博客。
我无法在全国范围内下载SRTM,进行镶嵌和裁剪工作。第一次迭代后,下载将停止,因为GeoTiff文件可以’找不到。我在Ubuntu 16.04上工作,我想问题在于getData()函数搜索pattern =“.TIF”我解压缩的图片以小写字母结尾=“.tif”。我试图获取getData()函数的源代码,但不幸的是它没有’工作。您在建议方面有什么工作吗?
非常感谢,
瓦伦丁
瓦伦丁·路易斯 在4年前
嗨,瓦伦丁!这很有趣,…我从未在Linux上尝试过该功能。
尝试在下载前手动设置工作目录。即`setwd(“您的/ custom / dir / for / srtm /”)`.
然后将getData()调用包装在一个try子句中:try(getData(“…”))`。这应该确保您完成所有迭代。
下载所有内容后,致电`rasters <- list.files(full.names=T, pattern=".tif$")` and stack them using `stack <- stack(rasters)`. Good luck!
马丁 在4年前
谢谢,马丁。
尝试条款起作用了,但是堆栈没有起作用’t是因为程度不同。栅格上的lapply()完成了这项工作。
最良好的祝愿,
瓦伦丁
瓦伦丁·路易斯 在4年前
很高兴得知您成功了!干杯马丁
马丁 在4年前
很棒的教程!一世’最近我一直在与各种地形数据集进行搏斗,这种方法似乎非常方便。关于如何合并多个国家的任何提示?还是更好的是,世界上的某个特定区域(例如中美洲和加勒比海)?一世’我试图用多个国家制作地形图’仍然种到海岸线上。
-霍利斯
霍利斯·达恩(Hollis Dahn) 在4年前
我不断收到此错误。请指导我:
穆罕默德·瓦卡斯(Muhammad Waqas) 3年前
嘿穆罕默德,
似乎shapefile的路径错误。按着这些次序:
1.确保您下载并解压缩了shapefile: http://www.citystoragesrq.com/wp-content/uploads/2017/01/srtm.zip
2.在计算机上找到shapefile的位置。例如:“C:/下载/srtm/tiles.shp”并将其复制到脚本中,而不是在此处提供的路径中:country_name shp <- shapefile("srtm/tiles.shp") Hope this works! C
马丁 3年前
我一直在尝试使用sf对象执行此工作流程,但我一直在努力了解如何从sf对象中提取lon / lat以在该函数内工作。
您是否建议在sp的sf对象中对此进行操作?
我的目标是为121个国家/地区的海岸线下载SRTM DEM。我也想在一秒钟内得到这个,但是现在尝试一下。
艾米 1年前
好的帖子,我有一个问题:当第二个国家在第一个大国家之内时,如何合并两个国家?如果可以的话
阿布 12个月前
嗨,马丁,
感谢您的代码和所有方便的信息。
我在挪威尝试了您的代码,但似乎一切正常,但是绘图图没有’似乎没有超过60°N,所以我才到达挪威最南端。有什么方法可以使用R中的相同方法获得高于60°N的数据?
非常感谢你!
约瑟芬
约瑟芬 9个月前
发表回复