Extracting SAVI

[fyousef]

Hello there,

Thanks for the excellent package of bfastSpatial. I have come across an issue

When I try to pick SAVI as the VI in the processLandsat function, it throws an error. Here is the thing, I've ordered SAVI from ESPA directly (with a L=0.5) and dont have reflectance layers in the tar file, so I'm only using the above function to extract the tiff and do further processing on it. I guess the function is still looking to do some sort of processing on reflectance layers b/c of the L argument! I tried to avoid the error by not providing the L with no luck. Thank you

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", : Cannot create a RasterLayer object from this file. In addition: There were 12 warnings (use warnings() to see them)

Warning messages: 1: In if (x == "" | x == ".") { ... : the condition has length > 1 and only the first element will be used 2: In if (!start %in% c("http", "ftp")) { ... : the condition has length > 1 and only the first element will be used 3: In if (fileext %in% c(".GRD", ".GRI")) { ... : the condition has length > 1 and only the first element will be used 4: In if (!file.exists(x)) { ... : the condition has length > 1 and only the first element will be used 5: In if ((fileext %in% c(".HE5", ".NC", ".NCF", ".NC4", ... : the condition has length > 1 and only the first element will be used 6: In if (fileext == ".GRD") { ... : the condition has length > 1 and only the first element will be used 7: In if (fileext == ".BIG" | fileext == ".BRD") { ... : the condition has length > 1 and only the first element will be used 8: In if (fileext %in% c(".BIN")) { ... : the condition has length > 1 and only the first element will be used 9: In if (fileext == ".DOC") { ... : the condition has length > 1 and only the first element will be used 10: In if (fileext %in% c(".SGRD", ".SDAT")) { ... : the condition has length > 1 and only the first element will be used 11: In if (nchar(filename) == 0) stop("empty file name") : the condition has length > 1 and only the first element will be used 12: In if (!file.exists(x)) { ... : the condition has length > 1 and only the first element will be used

[loicdtx]

@fyousef Did you pass 'savi' in lowercase (as it should be, although it would make sense to allow uppercase too)? Looking at the source code I don't see why this wouldn't work if you have a savi layer in the tar archive.

[fyousef]

Yes, I did pass it in lowercase. I have successfully processed NDVI and EVI in the same fashion.

[loicdtx]

@fyousef I tried with data I've recently dowloaded and I can't reproduce. See example below.

library(bfastSpatial)

archive <- '/media/loic/data/landsat-test/loic.dutrieux@wur.nl-10122016-110454-137/LE71960292015017-SC20161012155108.tar.gz'

# Ensure that archive contains savi layer (9)
untar(archive, list = TRUE)

[1] "LE71960292015017ASN00_sr_cloud_qa.tif"         
[2] "LE71960292015017ASN00_sr_band7.tif"            
[3] "LE71960292015017ASN00_sr_ddv_qa.tif"           
[4] "LE71960292015017ASN00_sr_ndmi.tif"             
[5] "LE71960292015017ASN00_sr_evi.tif"              
[6] "LE71960292015017ASN00_sr_snow_qa.tif"          
[7] "LE71960292015017ASN00_sr_ndvi.tif"             
[8] "LE71960292015017ASN00_sr_band5.tif"            
[9] "LE71960292015017ASN00_sr_savi.tif"             
[10] "LE71960292015017ASN00_sr_cloud_shadow_qa.tif"  
[11] "LE71960292015017ASN00_cfmask.tif"              
[12] "LE71960292015017ASN00_sr_land_water_qa.tif"    
[13] "LE71960292015017ASN00_sr_adjacent_cloud_qa.tif"
[14] "LE71960292015017ASN00_cfmask_conf.tif"         
[15] "LE71960292015017ASN00_sr_band1.tif"            
[16] "LE71960292015017ASN00_sr_band2.tif"            
[17] "LE71960292015017ASN00_sr_atmos_opacity.tif"    
[18] "LE71960292015017ASN00.xml"                     
[19] "LE71960292015017ASN00_sr_band4.tif"            
[20] "LE71960292015017ASN00_sr_band3.tif"            
[21] "LE71960292015017ASN00_sr_fill_qa.tif"     

processLandsat(archive, vi = 'savi', srdir = tmpDir(), outdir = '/home/loic/sandbox/', mask = 'fmask')

r <- raster('~/sandbox/savi.LE71960292015017.grd')
r

class       : RasterLayer 
dimensions  : 134, 133, 17822  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : -1993.67, 1996.33, -2012.403, 2007.597  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=aea +lat_1=39.78 +lat_2=49.78 +lat_0=44.78 +lon_0=4.71 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /home/loic/sandbox/savi.LE71960292015017.grd 
names       : layer 
values      : 638, 3304  (min, max)

sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=fr_FR.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] bfastSpatial_0.6.2 zoo_1.7-13         bitops_1.0-6       rgdal_1.1-10       stringr_1.1.0     
[6] gdalUtils_2.0.1.7  bfast_1.5.7        raster_2.5-8       sp_1.2-3          

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.6       magrittr_1.5      forecast_7.1      munsell_0.4.3     colorspace_1.2-6 
 [6] lattice_0.20-33   quadprog_1.5-5    foreach_1.4.3     plyr_1.8.4        tools_3.3.1      
[11] nnet_7.3-12       grid_3.3.1        timeDate_3012.100 gtable_0.2.0      R.oo_1.20.0      
[16] iterators_1.0.8   tseries_0.10-35   ggplot2_2.1.0     R.utils_2.4.0     codetools_0.2-14 
[21] strucchange_1.5-1 fracdiff_1.4-2    sandwich_2.3-4    stringi_1.1.1     R.methodsS3_1.7.1
[26] scales_0.4.0

[fyousef]

This is strange.

datadir` <- "E:/UCLA/Landsat/"
list <- list.files(datadir, full.names=TRUE, pattern = "*.gz")

list
[1] "E:/UCLA/Landsat/LT50270331993092-SC20161008165002.tar.gz"

untar(list, lis=TRUE)

[1] "LT50270331993092XXX02_sr_msavi.tif"    "LT50270331993092XXX02_sr_savi.tif"     "LT50270331993092XXX02_sr_evi.tif"     
[4] "LT50270331993092XXX02_cfmask_conf.tif" "LT50270331993092XXX02.xml"             "LT50270331993092XXX02_cfmask.tif"     
[7] "LT50270331993092XXX02_sr_ndvi.tif"  

processLandsat(list, vi='savi', srdir = tmpDir(), outdir="E:/UCLA/Landsat/", mask = 'fmask')

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  : 
Cannot create a RasterLayer object from this file.
In addition: There were 12 warnings (use warnings() to see them)

[loicdtx]

I see, savi is matching twice, against savi and msavi. I'll try to push a fix today.

[loicdtx]

@fyousef I updated the package. It should work now, after you re-install. Thanks for the bug report

[fyousef]

Thank you The very best