-
Notifications
You must be signed in to change notification settings - Fork 0
/
Download_DEM-file.R
52 lines (39 loc) · 1.41 KB
/
Download_DEM-file.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#Quelle: http://www.gis-blog.com/download-srtm-for-an-entire-country/
library(raster)
library(rgeos)
library(rasterVis)
#------------SETTINGS--------------
#Specify target ISO country code and path to downloaded shapefile
country_name <- "AUT" #Austria
shp <- shapefile("srtm/tiles.shp") #Path to SRTM Tiles (can be found in subfolder srtm)
#------------EXECUTE FROM HERE--------------
#Get country geometry first
country <- getData("GADM",
country = country_name,
level=0)
#Intersect country geometry with 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 <- try(getData('SRTM',
lon=lon,
lat=lat))
srtm_list[[i]] <- tile
}
# AB HIER FUNKTIONIERTS NICHT MEHR (auf Linux)
# #Mosaic
# srtm_list$fun <- mean
# srtm_mosaic <- do.call(mosaic, srtm_list)
#
# #Crop to country borders
# srtm_crop <- mask(srtm_mosaic, country)
#
# #Plot results
# p <- levelplot(srtm_crop)
# p + layer(sp.lines(country,
# lwd=0.8,
# col='darkgray'))