Intra-version improvements (#242)

## Minor changes - Add NDBI / NDVIre indices (#241) - Add footprint among metadata retrieved using `safe_getMetadata()` (from existing offline SAFE archives) and `s2_list()` (from SciHub online metadata) - Filter offline SAFE archives by footprint instead than by S2 tile ID, so skipping using images which do not contain any data for a specific Area Of Interest - Add method to convert class `safelist` to `sf` (using footprint) ## Fixes - Hide warning using PROJ >= 6 (rspatial/raster#78)
ranghetti · Nov 19, 2019 · 03aa65a · 03aa65a
1 parent fca9aeb
commit 03aa65a
Show file tree

Hide file tree

Showing 17 changed files with 288 additions and 144 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: sen2r
 Type: Package
 Title: Find, Download and Process Sentinel-2 Data
-Version: 1.2.1
+Version: 1.2.1.9001
 Authors@R: c(person("Luigi", "Ranghetti",
                     email = "luigi@ranghetti.info",
                     role = c("aut", "cre"),

diff --git a/NAMESPACE b/NAMESPACE
@@ -4,6 +4,7 @@ S3method(as.character,safelist)
 S3method(as.data.frame,safelist)
 S3method(as.data.table,safelist)
 S3method(print,safelist)
+S3method(st_as_sf,safelist)
 export(abs2rel)
 export(build_example_param_file)
 export(check_gdal)

diff --git a/R/create_indices_db.R b/R/create_indices_db.R
@@ -457,6 +457,26 @@ create_indices_db <- function(xslt_path = NA,
       s2_formula_mathml = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <mrow>\n  <mfrac>\n   <mrow>\n    <mrow>\n     <mrow>\n      <mi mathcolor=\"#443399\">GREEN</mi>\n      <mo>-</mo>\n      <mi mathcolor=\"#443399\">NIR</mi>\n     </mrow>\n    </mrow>\n   </mrow>\n   <mrow>\n    <mrow>\n     <mrow>\n      <mi mathcolor=\"#443399\">GREEN</mi>\n      <mo>+</mo>\n      <mi mathcolor=\"#443399\">NIR</mi>\n     </mrow>\n    </mrow>\n   </mrow>\n  </mfrac>\n </mrow>\n</math>",
       checked = TRUE,
       a = NA, b = NA, x = NA
+    ),
+    "NDVIre" = data.frame(
+      n_index = 321,
+      longname = "Red-edge-based Normalized Difference Vegetation Index",
+      name = "NDVIre",
+      link = "https://doi.org/10.1016/S0034-4257(03)00131-7",
+      s2_formula = "(band_5-band_4)/(band_5+band_4)",
+      s2_formula_mathml = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <mrow>\n  <mfrac>\n   <mrow>\n    <mrow>\n     <mrow>\n      <mi mathcolor=\"#443399\">Rededge1</mi>\n      <mo>-</mo>\n      <mi mathcolor=\"#443399\">RED</mi>\n     </mrow>\n    </mrow>\n   </mrow>\n   <mrow>\n    <mrow>\n     <mrow>\n      <mi mathcolor=\"#443399\">Rededge1</mi>\n      <mo>+</mo>\n      <mi mathcolor=\"#443399\">RED</mi>\n     </mrow>\n    </mrow>\n   </mrow>\n  </mfrac>\n </mrow>\n</math>",
+      checked = TRUE,
+      a = NA, b = NA, x = NA
+    ),
+    "NDBI" = data.frame(
+      n_index = 322,
+      longname = "Normalized Difference Built-up Index",
+      name = "NDBI",
+      link = "https://doi.org/10.1080/01431160304987",
+      s2_formula = "(band_11-band_8)/(band_11+band_8)",
+      s2_formula_mathml = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">\n <mrow>\n  <mfrac>\n   <mrow>\n    <mrow>\n     <mrow>\n      <mi mathcolor=\"#443399\">SWIR1</mi>\n      <mo>-</mo>\n      <mi mathcolor=\"#443399\">NIR</mi>\n     </mrow>\n    </mrow>\n   </mrow>\n   <mrow>\n    <mrow>\n     <mrow>\n      <mi mathcolor=\"#443399\">SWIR1</mi>\n      <mo>+</mo>\n      <mi mathcolor=\"#443399\">NIR</mi>\n     </mrow>\n    </mrow>\n   </mrow>\n  </mfrac>\n </mrow>\n</math>",
+      checked = TRUE,
+      a = NA, b = NA, x = NA
     )
   ), fill=TRUE)
 

diff --git a/R/s2_calcindices.R b/R/s2_calcindices.R
@@ -138,13 +138,16 @@ s2_calcindices <- function(infiles,
     clip <- function(x,min,max) {(x+min+2*max+abs(x-min)-abs(x+min-2*max+abs(x-min)))/4}
 
     out <- raster(x)
-    out <- writeStart(
-      out, out_file,
-      NAflag = NAflag,
-      datatype = datatype,
-      format = ifelse(sel_format=="VRT", "GTiff", sel_format),
-      if(sel_format %in% c("GTiff","VRT")){options = paste0("COMPRESS=",compress)},
-      overwrite = overwrite
+    suppress_warnings(
+      out <- writeStart(
+        out, out_file,
+        NAflag = NAflag,
+        datatype = datatype,
+        format = ifelse(sel_format=="VRT", "GTiff", sel_format),
+        if(sel_format %in% c("GTiff","VRT")){options = paste0("COMPRESS=",compress)},
+        overwrite = overwrite
+      ),
+      "NOT UPDATED FOR PROJ >\\= 6"
     )
 
     # x <- brick(infiles)

diff --git a/R/s2_list.R b/R/s2_list.R
@@ -327,27 +327,29 @@ s2_list <- function(spatial_extent = NULL,
           url <- in_entry[which(grepl("<link href=", in_entry))] %>%
             gsub("^.*<link href=\"([^\"]+)\"/>.*$", "\\1", .)
 
-          id_orbit <- in_entry[which(grepl("relativeorbitnumber", in_entry))] %>%
+          id_orbit <- in_entry[which(grepl("\\\"relativeorbitnumber\\\"", in_entry))] %>%
             gsub("^.*<int name=\"relativeorbitnumber\">([^<]+)</int>.*$", "\\1", .) %>%
             as.numeric() %>% sprintf("%03i", .)
 
-          clouds <- in_entry[which(grepl("cloudcoverpercentage", in_entry))] %>%
+          footprint <- tryCatch(
+            in_entry[which(grepl("\\\"footprint\\\"", in_entry))] %>%
+              gsub("^.*<str name=\"footprint\">([^<]+)</str>.*$", "\\1", .) %>%
+              st_as_sfc(crs = 4326),
+            error = function(e) {st_polygon()}
+          )
+
+          clouds <- in_entry[which(grepl("\\\"cloudcoverpercentage\\\"", in_entry))] %>%
             gsub("^.*<double name=\"cloudcoverpercentage\">([^<]+)</double>.*$", "\\1", .) %>%
             as.numeric()
 
-          proc_level <- in_entry[which(grepl("processinglevel", in_entry))] %>%
+          proc_level <- in_entry[which(grepl("\\\"processinglevel\\\"", in_entry))] %>%
             gsub("^.*<str name=\"processinglevel\">Level\\-([^<]+)</str>.*$", "\\1", .)
 
-          mission <- in_entry[which(grepl("platformserialidentifier", in_entry))] %>%
+          mission <- in_entry[which(grepl("\\\"platformserialidentifier\\\"", in_entry))] %>%
             gsub("^.*<str name=\"platformserialidentifier\">Sentinel\\-([^<]+)</str>.*$", "\\1", .)
 
-          # id_tile <- in_entry[which(grepl("name=\"tileid\"", in_entry))] %>%
-          #   gsub("^.*<str name=\"tileid\">([^<]+)</str>.*$", "\\1", .)
           id_tile <- gsub("^.+_T([0-9]{2}[A-Z]{3})_.+$", "\\1", title)
 
-          # sensing_datetime <- in_entry[which(grepl("name=\"endposition\"", in_entry))] %>%
-          #   gsub("^.*<date name=\"endposition\">([0-9\\-]+)T[0-9\\:\\.]+Z</date>.*$", "\\1", .) %>%
-          #   as.POSIXct()
           sensing_datetime <- gsub(
             "^S2[AB]\\_MSIL[12][AC]\\_([0-9]{8}T[0-9]{6})\\_N[0-9]{4}\\_R[0-9]{3}\\_T[A-Z0-9]{5}\\_[0-9]{8}T[0-9]{6}$",
             "\\1", title
@@ -362,6 +364,9 @@ s2_list <- function(spatial_extent = NULL,
             gsub("^.*<date name=\"ingestiondate\">([0-9\\-]+)T([0-9\\:\\.]+)Z</date>.*$", "\\1 \\2", .) %>%
             as.POSIXct(tz = "UTC")
 
+          uuid <- in_entry[which(grepl("\\\"uuid\\\"", in_entry))] %>%
+            gsub("^.*<str name=\"uuid\">([^<]+)</str>.*$", "\\1", .)
+
           # print(paste0(title, ".SAFE"))
           out_list[[n_entries]] <- data.frame(
             name = paste0(title, ".SAFE"),
@@ -373,6 +378,8 @@ s2_list <- function(spatial_extent = NULL,
             sensing_datetime = sensing_datetime,
             ingestion_datetime = ingestion_datetime,
             clouds = clouds,
+            footprint = st_as_text(footprint),
+            uuid = uuid,
             stringsAsFactors = FALSE
           )
           n_entries <- n_entries + 1

diff --git a/R/s2_mask.R b/R/s2_mask.R
@@ -400,21 +400,31 @@ s2_mask <- function(infiles,
             mask_tmpfiles,
             file.path(sel_tmpdir, basename(tempfile(pattern = "mask_", fileext = ".tif")))
           )
-          raster::calc(inmask[[j]],
-                       function(x){as.integer(!is.na(nn(x)) & !x %in% req_masks[[j]])},
-                       filename = mask_tmpfiles[j],
-                       options = "COMPRESS=LZW",
-                       datatype = "INT1U",
-                       overwrite = TRUE)
+          suppress_warnings(
+            raster::calc(
+              inmask[[j]],
+              function(x){as.integer(!is.na(nn(x)) & !x %in% req_masks[[j]])},
+              filename = mask_tmpfiles[j],
+              options = "COMPRESS=LZW",
+              datatype = "INT1U",
+              overwrite = TRUE
+            ),
+            "NOT UPDATED FOR PROJ >\\= 6"
+          )
           naval_tmpfiles <- c(
             naval_tmpfiles,
             file.path(sel_tmpdir, basename(tempfile(pattern = "naval_", fileext = ".tif")))
           )
-          raster::calc(inmask[[j]],
-                       function(x){as.integer(!is.na(nn(x)))},
-                       filename = naval_tmpfiles[j],
-                       options = "COMPRESS=LZW",
-                       datatype = "INT1U")
+          suppress_warnings(
+            raster::calc(
+              inmask[[j]],
+              function(x){as.integer(!is.na(nn(x)))},
+              filename = naval_tmpfiles[j],
+              options = "COMPRESS=LZW",
+              datatype = "INT1U"
+            ),
+            "NOT UPDATED FOR PROJ >\\= 6"
+          )
         }
         if(length(mask_tmpfiles)==1) {
           outmask <- mask_tmpfiles
@@ -519,18 +529,21 @@ s2_mask <- function(infiles,
               }
               if (any(!file.exists(binmask), overwrite == TRUE)) {
                 # mask NA values
-                raster::mask(
-                  raster(outmask_smooth),
-                  raster(outnaval_res),
-                  filename = binmask,
-                  maskvalue = 0,
-                  updatevalue = sel_naflag,
-                  updateNA = TRUE,
-                  NAflag = 255,
-                  datatype = "INT1U",
-                  format = sel_format,
-                  options = if(sel_format == "GTiff") {paste0("COMPRESS=",compress)},
-                  overwrite = overwrite
+                suppress_warnings(
+                  raster::mask(
+                    raster(outmask_smooth),
+                    raster(outnaval_res),
+                    filename = binmask,
+                    maskvalue = 0,
+                    updatevalue = sel_naflag,
+                    updateNA = TRUE,
+                    NAflag = 255,
+                    datatype = "INT1U",
+                    format = sel_format,
+                    options = if(sel_format == "GTiff") {paste0("COMPRESS=",compress)},
+                    overwrite = overwrite
+                  ),
+                  "NOT UPDATED FOR PROJ >\\= 6"
                 )
               }
             }
@@ -552,14 +565,17 @@ s2_mask <- function(infiles,
               if (grepl("\\.vrt$", out_file)) {
                 out_file <- gsub("\\.vrt$", ".tif", out_file)
               }
-              out <- writeStart(
-                out,
-                out_file,
-                NAflag=na,
-                datatype = datatype,
-                format = ifelse(sel_format=="VRT","GTiff",sel_format),
-                if(sel_format %in% c("GTiff","VRT")){options = c("COMPRESS=LZW")},
-                overwrite = overwrite
+              suppress_warnings(
+                out <- writeStart(
+                  out,
+                  out_file,
+                  NAflag=na,
+                  datatype = datatype,
+                  format = ifelse(sel_format=="VRT","GTiff",sel_format),
+                  if(sel_format %in% c("GTiff","VRT")){options = c("COMPRESS=LZW")},
+                  overwrite = overwrite
+                ),
+                "NOT UPDATED FOR PROJ >\\= 6"
               )
               #4 bytes per cell, nb + 1 bands (brick + mask), * 2 to account for a copy
               bs <- blockSize(out, minblocks = 8)

diff --git a/R/safe_getMetadata.R b/R/safe_getMetadata.R
@@ -241,7 +241,7 @@ safe_isvalid <- function(s2, allow_oldnames = FALSE, check_file = TRUE) {
                  "id_orbit", "orbit_number", "sensing_datetime", "id_baseline") # information GENERALLY retrieved from name
   info_gdal <- c("clouds","direction","orbit_n","preview_url", # information retrieved by reading the file metadata
                  "proc_baseline","gdal_level","gdal_sensing_datetime",
-                 "nodata_value","saturated_value")
+                 "nodata_value","saturated_value","footprint")
 
   # check format attribute
   if (format == "default") {
@@ -717,45 +717,51 @@ safe_isvalid <- function(s2, allow_oldnames = FALSE, check_file = TRUE) {
       # }
 
       # retrieve metadata[[i]] from file content
-      if (exists("s2_gdal") & s2_exists) {
-
+      if (
+        exists("s2_gdal") & 
+        s2_exists &
+        any(info_gdal %in% sel_info)
+      ) {
         # Read metadata[[i]]
+        s2_gdal_metadata <- py_to_r(s2_gdal$GetMetadata())
         if ("clouds" %in% sel_info) {
-          metadata[[i]][["clouds"]] <- py_to_r(s2_gdal$GetMetadata()[["CLOUD_COVERAGE_ASSESSMENT"]])
+          metadata[[i]][["clouds"]] <- s2_gdal_metadata[["CLOUD_COVERAGE_ASSESSMENT"]]
         }
         if ("direction" %in% sel_info) {
-          metadata[[i]][["direction"]] <- py_to_r(s2_gdal$GetMetadata()[["DATATAKE_1_SENSING_ORBIT_DIRECTION"]])
+          metadata[[i]][["direction"]] <- s2_gdal_metadata[["DATATAKE_1_SENSING_ORBIT_DIRECTION"]]
         }
         if ("orbit_n" %in% sel_info) {
-          metadata[[i]][["orbit_n"]] <- py_to_r(s2_gdal$GetMetadata()[["DATATAKE_1_SENSING_ORBIT_NUMBER"]])
+          metadata[[i]][["orbit_n"]] <- s2_gdal_metadata[["DATATAKE_1_SENSING_ORBIT_NUMBER"]]
         }
         if ("preview_url" %in% sel_info) {
-          metadata[[i]][["preview_url"]] <- py_to_r(s2_gdal$GetMetadata()[["PREVIEW_IMAGE_URL"]])
+          metadata[[i]][["preview_url"]] <- s2_gdal_metadata[["PREVIEW_IMAGE_URL"]]
         }
         if ("proc_baseline" %in% sel_info) {
-          metadata[[i]][["proc_baseline"]] <- py_to_r(s2_gdal$GetMetadata()[["PROCESSING_BASELINE"]])
+          metadata[[i]][["proc_baseline"]] <- s2_gdal_metadata[["PROCESSING_BASELINE"]]
         }
         # if ("level" %in% sel_info) {
-        #   metadata[[i]][["level"]] <- py_to_r(s2_gdal$GetMetadata()[["PROCESSING_LEVEL"]])
+        #   metadata[[i]][["level"]] <- s2_gdal_metadata[["PROCESSING_LEVEL"]]
         # }
         if ("sensing_datetime" %in% sel_info) {
           start_time <- as.POSIXct(
-            py_to_r(s2_gdal$GetMetadata()[["PRODUCT_START_TIME"]]), format="%Y-%m-%dT%H:%M:%S", tz="UTC")
+            s2_gdal_metadata[["PRODUCT_START_TIME"]], format="%Y-%m-%dT%H:%M:%S", tz="UTC")
           stop_time <- as.POSIXct(
-            py_to_r(s2_gdal$GetMetadata()[["PRODUCT_STOP_TIME"]]), format="%Y-%m-%dT%H:%M:%S", tz="UTC")
+            s2_gdal_metadata[["PRODUCT_STOP_TIME"]], format="%Y-%m-%dT%H:%M:%S", tz="UTC")
           metadata[[i]][["sensing_datetime"]] <- if (start_time == stop_time) {
             start_time
           } else {
             c(start_time, stop_time)
           }
         }
         if ("nodata_value" %in% sel_info) {
-          metadata[[i]][["nodata_value"]] <- py_to_r(s2_gdal$GetMetadata()[["SPECIAL_VALUE_NODATA"]])
+          metadata[[i]][["nodata_value"]] <- s2_gdal_metadata[["SPECIAL_VALUE_NODATA"]]
         }
         if ("saturated_value" %in% sel_info) {
-          metadata[[i]][["saturated_value"]] <- py_to_r(s2_gdal$GetMetadata()[["SPECIAL_VALUE_SATURATED"]])
+          metadata[[i]][["saturated_value"]] <- s2_gdal_metadata[["SPECIAL_VALUE_SATURATED"]]
+        }
+        if ("footprint" %in% sel_info) {
+          metadata[[i]][["footprint"]] <- s2_gdal_metadata[["FOOTPRINT"]]
         }
-
       }
 
     } # end of s2_exists IF cycle

diff --git a/R/safelist-class.R b/R/safelist-class.R
@@ -36,11 +36,16 @@
 #' ## Convert to other classes
 #' (s2_char <- as.character(list_safe)) # convert to a simple named character
 #' (s2_df <- as.data.frame(list_safe)) # convert to a data.frame
-#' library(data.table); (s2_dt <- as.data.table(list_safe)) # convert to a data.table
+#' library(data.table)
+#' (s2_dt <- as.data.table(list_safe)) # convert to a data.table
+#' library(sf)
+#' (s2_sf <- st_as_sf(list_safe)) # convert to sf
 #' 
 #' ## Convert from other classes
 #' as(s2_char, "safelist") # this causes the loss of hidden attributes
-#' as(s2_df, "safelist") # this maintains attributes as columns
+#' as(s2_df, "safelist") # this (and followings) maintain attributes as columns
+#' as(s2_dt, "safelist")
+#' as(s2_sf, "safelist")
 #' }
 
 setClass("safelist", contains = "character")
@@ -89,12 +94,15 @@ setAs("data.frame", "safelist", function(from) {
   as(to, "safelist")
 })
 
+setAs("sf", "safelist", function(from) {
+  as(as.data.frame(s2_sf), "safelist")
+})
 
 ## Methods FROM safelist
 
 #' @export
 as.data.frame.safelist <- function(x, row.names = NULL, optional = FALSE, ...) {
-  to <- data.frame(name = names(x), url = as.vector(x))
+  to <- data.frame(name = names(x), url = as.vector(x), stringsAsFactors = FALSE)
   autoRN <- (is.null(row.names) || length(row.names) != nrow(to))
   attr(to, "row.names") <- if (autoRN) {seq_len(nrow(to))} else {row.names}
 
@@ -127,37 +135,20 @@ setAs("safelist", "character", function(from) {
   as.character(from)
 })
 
+#' @export
+st_as_sf.safelist <- function(x, ...) {
+  if (!is.null(attr(x, "footprint"))) {
+    sf::st_as_sf(as.data.frame(x), wkt = "footprint", crs = 4326)
+  } else {
+    stop("cannot convert to safelist (missing footprint)")
+  }
+}
+setAs("safelist", "sf", function(from) {
+  st_as_sf(from)
+})
 
-## Print methods
-
-# #' @export
-# print.s2dt = function(x) {
-#   if (!all(c("name","url") %in% names(x))) {
-#     print(data.table(x))
-#   } else {
-#     printed_cols <- c("name", "url")
-#     x_print <- as.data.frame(x)[seq_len(min(nrow(x),5)), printed_cols]
-#     x_print$url <- paste0(substr(x_print$url,1,20),"...")
-#     cat("A data.table with", nrow(x), "SAFE archives.\n")
-#     print.data.frame(x_print)
-#     if (nrow(x) > 5 | ncol(x) > ncol(x_print)) {
-#       cat("...with")
-#       if (nrow(x) > 5) {
-#         cat("", nrow(x)-5, "more rows")
-#       }
-#       if (nrow(x) > 5 & ncol(x) > ncol(x_print)) {
-#         cat(" and")
-#       }
-#       if (ncol(x) > ncol(x_print)) {
-#         cat("", ncol(x)-ncol(x_print), "more columns: ")
-#         cat(paste(names(x)[!names(x) %in% printed_cols], collapse = ", "))
-#       }
-#       cat(".\n")
-#     }
-#   }
-#   invisible(x)
-# }
 
+## Print method
 #' @export
 print.safelist = function(x, ...) {
   x_print <- as.character(x)[seq_len(min(length(x),5))]