-
Notifications
You must be signed in to change notification settings - Fork 1
/
elev_grab.R
102 lines (88 loc) · 4.8 KB
/
elev_grab.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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#' Elevation Grabber - From NTS Grids
#'
#' This function uses your reference layer to define the NTS grids necessary to pull elevation from the CDEM layers provided by the Government of Canada.
#' @note Requires an internet connection
#'
#' @param reference_grid This is a reference raster to provide a projection and a surface to assign values onto, this should be a grid that registers with the other grids you are using for your project. Can be either the location of the raster or a raster object.
#' @param output_directory Output directory where your elevation and Wind Ninja elevation grids will be stored.
#' @param aoi_e Polygon for the area of interest to be intersected with the NTS grid layer. _(Default: "")_ If Default is used the full extent under the reference_grid will be returned.
#'
#' @details The purpose of this function is to generate a common and rapid elevation layer that is sampled and masked to the reference grid for use within Burn-P3. A second elevation grid is also generated for use in Wind Ninja as that software will fail with NA values in the elevation grid.
#'
#' @importFrom terra rast crop merge writeRaster as.polygons project mask crs
#' @importFrom sf st_as_sf st_read st_crop st_transform
#'
#' @export
#'
#' @references [Open Data Canada](https://open.canada.ca/data/en/dataset/7f245e4d-76c2-4caa-951a-45d1d2051333)
#' @references [Open Data Canada FTP](http://ftp.maps.canada.ca/pub/nrcan_rncan/vector/index/)
#' @references [Wind Ninja](https://www.firelab.org/project/windninja)
#'
#'
#' @examples
#' ## Load example data
#' ref_grid <- rast(system.file("extdata/fuel.tif",package = "BurnP3.HelpR"))
#' temp_dir <- tempdir()
#' test <- elev_grab(reference_grid = ref_grid,
#' output_directory = temp_dir)
#'
#'
#' unlink(temp_dir, recursive = T)
elev_grab <- function(aoi_e = NULL, reference_grid,output_directory){
if ( any(grepl("SpatRaster", class(reference_grid))) ) { grast <- reference_grid }
if ( any(grepl("character", class(reference_grid))) ) { grast <- terra::rast(reference_grid) }
if ( any(!grepl("SpatRaster|character", class(reference_grid))) ) { message("Reference Grid must be the directory of the spatraster or a spatraster object.") }
if( is.null(aoi_e) ) { aoi_e <- sf::st_as_sf(as.polygons(grast, extent=T))}
if ( any(grepl("sf", class(aoi_e))) ) { aoi_e <- aoi_e }
if ( any(grepl("character", class(aoi_e))) ) { aoi_e <- sf:read_sf(aoi_e) }
if ( !any(grepl("sf|character", class(aoi_e))) ) { message("aoi_e must be a simple feature (sf) or a directory to a simple feature.") }
## Download and extract the NTS grid specified by user (will be removed after use), uses the 250k grid as that is what CDEM is based on
nts_temp <- tempfile(fileext = '.zip')
download.file(destfile = nts_temp,url = "http://ftp.maps.canada.ca/pub/nrcan_rncan/vector/index/nts_snrc.zip")
nts_files <- unzip(nts_temp,list = T)
unzip(zipfile = nts_temp,
files = nts_files$Name,
exdir = gsub(".zip","",nts_temp))
nts_grid <- sf::st_read(dsn = gsub(".zip","",nts_temp),
layer = "nts_snrc_250k")
## Determine the NTS grids the data exists across
layers <- sf::st_crop(nts_grid,
sf::st_transform(aoi_e,crs = st_crs(nts_grid)))$NTS_SNRC
## Extract the CDEM tiles needed to generate the grid.
elevation <- lapply(layers,function(i){
loc <- tempfile(fileext = ".zip")
loc <- gsub("\\\\","/",loc)
download.file(url = paste0("http://ftp.geogratis.gc.ca/pub/nrcan_rncan/elevation/cdem_mnec/",
substr(i,1,3),
"/cdem_dem_",
i,
"_tif.zip"),
destfile = loc)
files <- unzip(loc,list = T)
terra::rast(paste0("/vsizip/",loc,"/",files$Name[grep(".tif$",files$Name)]))
}
)
if (length(elevation) > 1) {
mosaic.r <- do.call(merge,elevation)
} else{
mosaic.r <- elevation[[1]]
}
mosaic.r <- terra::crop(terra::project(x = mosaic.r, y = grast),grast)
unlink(c(gsub(".zip","",nts_temp),list.files(tempdir(),pattern = ".zip",full.names = T)),recursive = T)
if (class(aoi_e)[1] == "sf") {
terra::writeRaster(terra::mask(mosaic.r,vect(aoi_e)),
paste0(output_directory,"elevation.tif"),
wopt = list(filetype = "GTiff",
datatype = "INT2S",
gdal = c("COMPRESS=DEFLATE","ZLEVEL=9","PREDICTOR=2")),
NAflag = -9999,
overwrite = T)
}
terra::writeRaster(mosaic.r,
paste0(output_directory,"elevation_wn.tif"),
wopt = list(filetype = "GTiff",
datatype = "INT2S",
gdal = c("COMPRESS=DEFLATE","ZLEVEL=9","PREDICTOR=2")),
NAflag = -9999,
overwrite = T)
}