#' @import ggplot2
#' @import dplyr
#' @importFrom textreadr read_html
#' @importFrom utils download.file unzip
#' @importFrom raster as.data.frame raster merge
#' @export
SRTM.list <- function(resolution, want.plot = TRUE)
{
cat("retrieving directory listing from 'https://dds.cr.usgs.gov/'...")
res <- paste0("SRTM", resolution)
URL <- paste("https://dds.cr.usgs.gov/srtm/version2_1", res, sep = "/")
page <- textreadr::read_html(URL)
dir <- page[-c(1,2)]
sub_dir <- list()
j = 1
for(i in 1:length(dir))
{
if(substr(dir[i], nchar(dir[i]), nchar(dir[i])) == "/")
{
sub_dir[[j]] <- textreadr::read_html(paste(URL, dir[i], sep = "/"))[-c(1,2)]
names(sub_dir)[[j]] <- gsub('.{1}$', '', dir[i])
j <- j+1
}
}
files <- sapply(1:length(sub_dir), function(x) paste(names(sub_dir)[x], sub_dir[[x]], sep = "/"))
cat("Done~")
if(want.plot)
{
data.plot<- NULL
for(i in 1:length(files))
{
loc <- strsplit(files[[i]], split = "/")
group <- loc[[1]][1]
lat <- sapply(loc, function(x) substr(x[[2]], 1, 3))
lat <- ifelse(substr(lat, 1, 1) == "N", as.numeric(substr(lat, 2, 3)), -as.numeric(substr(lat, 2, 3)))
lon <- sapply(loc, function(x) substr(x[[2]], 4, 8))
lon <- ifelse(substr(lon, 1, 1) == "E", as.numeric(substr(lon, 2, 4)), -as.numeric(substr(lon, 2, 4)))
data.plot <- rbind(data.plot, data.frame(lat = lat, lon = lon, group = group))
}
WorldData <- map_data('world')
WorldData %>% filter(WorldData$region != "Antarctica") -> WorldData
WorldData <- fortify(WorldData)
east <- ifelse(resolution == 1, -50, 180)
south <- ifelse(resolution == 1, 10, -60)
p <- ggplot()
p <- p + geom_map(map=WorldData, fill=NA, colour="gray30", size=0.5) +
geom_tile(data=data.plot, aes(lon+0.5,lat+0.5, color = group), alpha=0, size = 0.1) +
coord_fixed() +
scale_x_continuous(limits = c(-180, east), expand = c(0, 0)) +
scale_y_continuous(limits = c(south, 60), expand = c(0, 0)) +
xlab(expression(paste("Longitude [", degree, "]", sep = ""))) +
ylab(expression(paste("Latitude [", degree, "]", sep = ""))) +
theme_bw() +
theme(plot.margin = unit(c(0,0,0,0), "lines"), legend.position = "bottom", legend.box.spacing = unit(c(-0.5,0,0.5,0), "lines"), legend.key.size = unit(0.3, "cm"), text = element_text(size = 7), legend.background = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.text.x = element_text(margin = margin(0.2,0,0.2,0, "lines")) )
print(p)
}
files
}
#listSRTM(1)
#' @export
SRTM.get <- function(resolution, files, directory = "data-raw")
{
res <- paste0("SRTM", resolution)
#if(missing(zone))
#stop("zone must be defined! Use function 'listSRTM' to check for available zones")
if(missing(files))
stop("file names must be defined! Use function 'listSRTM' to check for available files")
if(!file.exists(directory))
dir.create(directory)
for(i in 1:length(files))
{
URL <- paste("https://dds.cr.usgs.gov/srtm/version2_1", res, files[i], sep = "/")
dest <- paste(directory, strsplit(files[i], split = "/")[[1]][2], sep = "/")
utils::download.file(URL, destfile = dest)
utils::unzip(dest, exdir = directory)
unlink(dest)
}
}
#' @export
SRTM.read <- function(files, as.data.frame = FALSE)
{
cat("Warning: the files are big, you may want to read one file at a time\n")
if(length(files)>1)
{
tmp <- list()
for(i in 1:length(files))
{
cat("Reading", files[i], "...\n")
tmp[[i]] <- raster::raster(files[i])
}
elevation <- do.call(raster::merge, tmp)
}else{
elevation <- raster::raster(files[i])
}
#set voids to NA
raster::values(elevation)[raster::values(elevation) < -10994 | raster::values(elevation) > 8848] = NA
if(as.data.frame)
{
elevation <- raster::as.data.frame(elevation, xy = TRUE)
}
elevation
}
#setwd("/Users/DYang/Git/SolarData/data-raw")
#files <- dir()
#raster::image(elevation)
#raster::image(df2)
#zone <- "North_America"
#files <- c("North_America/N16W090.hgt.zip", "North_America/N15W090.hgt.zip", "North_America/N15W091.hgt.zip", "North_America/N15W092.hgt.zip")
#SRTM.get(resolution = 3, files = files)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.