R/skywatch.R

In fionahevans/agric: Functions for accessing & handling agricultural data.

#' Query the SkyWatch API for satellite imagery and climate/atmospheric datasets
#' 
#' This is a slight modification of the SkyWatchr function 'querySW' that checks for 
#' for a lack of responses. See the querySW help for details on usage.
#'
#' @param api_key	 string; personal alphanumeric API key. 
#' @param time_period	 string; one or two UTC timestamps in ISO format (yyyy-mm-ddThh:mm:ss.sssss+|-zzzz). 
#' @param longitude_latitude	string or an Spatial object. 
#' @param instrument_satellite	 string; source of the data, either the instrument on-board the satellite or the satellite itself. 
#' @param data_level	 numeric; data processing levels for Earth observation data. 
#' @param max_resolution	 numeric; maximum spatial resolution (in meters). 
#' @param max_cloudcover	 numeric; maximum cloud cover (in percentage). 
#' @param wavelength_band	 string; wavelength bands for imagery (i.e. Landsat-8) and by file type for non-imagery data (e.g. Hierarchical-Data-Format). 
#'
#' @author Fiona Evans
#' 
#' @return string; either "data.frame" (default) or "html". "html" returns a data.frame and prints it as html.
#' @export
swQuery <-function (api_key = NULL, time_period, longitude_latitude, instrument_satellite = NULL, 
                    data_level = NULL, max_resolution = NULL, max_cloudcover = NULL, 
                    wavelength_band = NULL, output = "data.frame") 
{
  if (is.null(api_key)) {
    api_key <- getOption("SkyWatchr.apikey")
    if (is.null(api_key)) {
      message("You need to set an API key via options(SkyWatchr.apikey = 'your_api_key')")
    }
  }
  if (is.null(time_period)) 
    time_period <- Sys.Date()
  if (is(longitude_latitude, "Spatial")) 
    longitude_latitude <- paste0(bbox(longitude_latitude), 
                                 collapse = ",")
  URL <- paste0("https://api.skywatch.co/data", "/time/", time_period, 
                "/location/", longitude_latitude)
  if (!is.null(instrument_satellite)) 
    URL <- paste0(URL, "/source/", instrument_satellite)
  if (!is.null(data_level)) 
    URL <- paste0(URL, "/level/", data_level)
  if (!is.null(max_resolution)) 
    URL <- paste0(URL, "/resolution/", max_resolution)
  if (!is.null(max_cloudcover)) 
    URL <- paste0(URL, "/cloudcover/", max_cloudcover)
  if (!is.null(wavelength_band)) 
    URL <- paste0(URL, "/band/", wavelength_band)
  query <- GET(paste0(URL, "?"), add_headers(Accept = "application/json", 
                                             `x-api-key` = api_key))
  res <- content(query)
  res <- lapply(res, unlist)
  
  res <- as.data.frame(do.call(rbind, lapply(res, function(x) {
    x
  })))
  
  # res <- as.data.frame(do.call(rbind, lapply(res, function(x) {
  #   x[!grepl("area.coord", names(x))]
  # })))
  
  if (nrow(res) > 0) {
  
    res[] <- lapply(res, as.character)
    res$size_kb <- round(as.numeric(res$size)/1000, 1)
    res$cloud_cover <- as.numeric(res$cloud_cover)
    res$resolution <- as.numeric(res$resolution)
    
    # Bounding box
    res1 <- res[, which(colnames(res) == "area.bbox1"):which(colnames(res) == 
                                                               "area.bbox4")]
    area <- apply(res1, 1, function(x) {
      x.lb <- paste(x[1], x[2], collapse = " ")
      x.lt <- paste(x[1], x[4], collapse = " ")
      x.rt <- paste(x[3], x[4], collapse = " ")
      x.rb <- paste(x[3], x[2], collapse = " ")
      return(paste0("POLYGON((", paste(x.lb, x.lt, x.rt, x.rb, 
                                       x.lb, sep = ","), "))"))
    })
    res <- res[, (max(grep("area.bbox", colnames(res))) + 1):ncol(res)]
    res$area <- area
    
    # Data polygon
    res2 <- res[, which(colnames(res) == "area.coordinates1"):which(colnames(res) == 
                                                                      "area.coordinates10")]
    area2 <- apply(res2, 1, function(x) {
      x.1 <- paste(x[1], x[2], collapse = " ")
      x.2 <- paste(x[3], x[4], collapse = " ")
      x.3 <- paste(x[5], x[6], collapse = " ")
      x.4 <- paste(x[7], x[8], collapse = " ")
      x.5 <- paste(x[7], x[8], collapse = " ")
      return(paste0("POLYGON((", paste(x.1, x.2, x.3, x.4, 
                                       x.5, sep = ","), "))"))
    })
    res <- res[, (max(grep("area.coordinates", colnames(res))) + 1):ncol(res)]
    res$polygon <- area2
    
    
    if (output == "html") {
      html.res <- interactiveTable(res)
      for (i in 1:nrow(res)) {
        html.res <- gsub(res[i, "download_path"], paste0("<a href='", 
                                                         res[i, "download_path"], "'>", res[i, "download_path"], 
                                                         "</a>"), html.res, fixed = TRUE)
      }
      print(html.res)
    }
  }
  return(res)
}



#' Convert time from SkyWatch API query to a Date object
#' 
#' Convert time from SkyWatch API query to a Date object.
#'
#' @param res Query response returned by SkyWatchr function 'querySW'.
#' @param index Index of element in x to convert.
#'
#' @author Fiona Evans
#' 
#' @return Date object.
#' @export
swDate <- function(res, index) {
  d <- as.character(res[index, "time"])
  d <- substr(d, 2, 11) 
  as.Date(d, "%Y-%m-%d")
}

#' Download Landsat data using the SkyWatch API and save bands as a raster stack 
#' 
#' Download satellite imagery and climate/atmospheric datasets using the SkyWatch API 
#' based on a query output object obtained from the querySW function. The API returns
#' data with coordinate reference 
#' CRS=" +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0". 
#' This function returns either the entire scene if region is not specified, or part of the scene 
#' cropped to the specified region with the same coordinate reference. The downloaded data 
#' are returned as a raster stack Object.
#'
#' @param res Query response returned by SkyWatchr function 'querySW'.
#' @param datetime Time of imafe data as specified in res.
#' @param region raster object specifying the extent to crop.
#' @param download If T, data are downloaded. Specify download=F if the data have already been downloaded.
#' @param dir Path of directory for downloaded data.
#' @param bands Vector of bands to download.
#' 
#' @return Returns a raster stack.
swDownload <- function (res, subset, dir, api_key) 
{
  x <- res
  
  if (is.null(api_key)) {
    api_key <- getOption("SkyWatchr.apikey")
    if (is.null(api_key)) {
      message("You need to set an API key")
    }
  }
  if (missing(subset)) {
    r <- rep_len(TRUE, nrow(x))
  }
  else {
    e <- substitute(subset)
    r <- eval(e, x, parent.frame())
  }
  x <- x[r, ]
  for (i in 1:nrow(x)) {
    cat(paste0("Downloading ", i, "/", nrow(x), " ...", "\n"))
    urli <- x[i, "download_path"]
    file_name <- paste0(x[i, "source"], "_", strsplit(urli, 
                        "/")[[1]][5], "_", x[i, "band"], "_", 
                        strsplit(urli, "/")[[1]][4], ".", x[i, "type"])
    GET(x[i, "download_path"], add_headers(`x-api-key` = api_key), 
        write_disk(paste0(dir, file_name)))
  }
  cat(paste0("\n", "All requested files have been downloaded to: ",  getwd(), "\n"))
}

#' Download Landsat data using the SkyWatch API and save bands as a raster stack 
#' 
#' Download satellite imagery and climate/atmospheric datasets using the SkyWatch API 
#' based on a query output object obtained from the querySW function. The API returns
#' data with coordinate reference 
#' CRS=" +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0". 
#' This function returns either the entire scene if region is not specified, or part of the scene 
#' cropped to the specified region with the same coordinate reference. The downloaded data 
#' are returned as a raster stack Object.
#'
#' @param res Query response returned by SkyWatchr function 'querySW'.
#' @param datetime Time of imafe data as specified in res.
#' @param region raster object specifying the extent to crop.
#' @param download If T, data are downloaded. Specify download=F if the data have already been downloaded.
#' @param dir Path of directory for downloaded data.
#' @param bands Vector of bands to download.
#'
#' @author Fiona Evans
#' 
#' @return Returns a raster stack.
#' @export
getLandsat8UTM <- function(res, datetime, region=NULL, download=T, dir=getwd(), api_key,
                       bands=c("Blue", "Green", "Red", "Near-Infrared",
                               "Short-Wave-Infrared-1",  "Short-Wave-Infrared-2")) {
  
  tmp <- subset(res, time == datetime & band %in% bands)
  
  # Download, each takes ~5 minutes
  if (download) swDownload(tmp, dir=dir, api_key=api_key)
  
  
  # Create rasterStack and data frame
  rasterList <- list()
  
  for (i in 1:length(bands)) {
    
    x <- subset(tmp, band == bands[i])
    
    urli <- x[, "download_path"]
    file_name <- paste0(x[, "source"], "_", 
                        strsplit(urli, "/")[[1]][5], "_", 
                        x[, "band"], "_", 
                        strsplit(urli,"/")[[1]][4], ".", x[, "type"])                     
    
    y <- raster(file_name)
    names(y) <- bands[i]
    
    if (!is.null(region)) y <- crop(y, region)
    
    rasterList[[bands[i]]] <- y
  }
  
  stack(rasterList)
}

#' Download Landsat data using the SkyWatch API and save bands as a raster stack 
#' 
#' Download satellite imagery and climate/atmospheric datasets using the SkyWatch API 
#' based on a query output object obtained from the querySW function.  
#' This function returns the part of the scene 
#' cropped and reprojected to the specified region with the coordinate reference
#' CRS="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0". The downloaded data 
#' are returned as a raster stack Object.
#'
#' @param res Query response returned by SkyWatchr function 'querySW'.
#' @param datetime Time of imafe data as specified in res.
#' @param region raster object specifying the extent to crop.
#' @param download If T, data are downloaded. Specify download=F if the data have already been downloaded.
#' @param dir Path of directory for downloaded data.
#' @param bands Vector of bands to download.
#'
#' @author Fiona Evans
#' 
#' @return Returns a raster stack.
#' @export
getLandsat8LL <- function(res, datetime, region, download=T, dir=getwd(), api_key,
                          bands=c("Blue", "Green", "Red", "Near-Infrared",
                                  "Short-Wave-Infrared-1",  "Short-Wave-Infrared-2")) {
  
  tmp <- subset(res, time == datetime & band %in% bands)
  
  # Download, each takes ~5 minutes
  if (download) swDownload(tmp, dir=dir, api_key=api_key)
  
  
  # Create rasterStack and data frame
  rasterList <- list()
  
  for (i in 1:length(bands)) {
    
    x <- subset(tmp, band == bands[i])
    
    urli <- x[, "download_path"]
    file_name <- paste0(x[, "source"], "_", 
                        strsplit(urli, "/")[[1]][5], "_", 
                        x[, "band"], "_", 
                        strsplit(urli,"/")[[1]][4], ".", x[, "type"])                     
    
    y <- raster(file_name)
    names(y) <- bands[i]
    
    # Reproject region (LL) to UTM and crop y 
    r0 <- projectRaster(region, y)
    rs <- crop(y, extent(r0))
    
    # Project y to LL
    y <- projectRaster(y, region)
    
    rasterList[[bands[i]]] <- y
  }
  
  stack(rasterList)
}