Error in runGDAL: gdal_utils warp

[caiomattos]

Hello,

I've been trying to use MODIS to download EVI data but I keep getting the same error. I thought it could be a problem with my code so I tried the runGDAL example in the documentation but the error persists. I tried searching here (and elsewhere) for a solution but couldn't find anything. I appreciate the help! Thanks in advance.

Code is:

runGdal( job="LSTaustria", product="MOD11A1", extent="Austria", begin="2010001", end="2010005",
        SDSstring="101", outProj="+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")

And (entire) error message:

########################
outProj          =  +proj=utm +zone=34 +datum=WGS84 +units=m +no_defs  (if applicable, derived from Raster*/Spatial*/sf* object)
pixelSize        =  asIn  (if applicable, derived from Raster* object)
resamplingType   =  near 
Output Directory =  /tmp/Rtmpl5PXqr/MODIS_ARC/PROCESSED/LSTaustria 
########################
Local structure is up-to-date. Using offline information!
Error in sf::gdal_utils(util = "warp", source = gdalSDS, destination = ofile, : gdal_utils warp: an error occured
Traceback:

1. runGdal(job = "LSTaustria", product = "MOD11A1", extent = "Austria", 
 .     begin = "2010001", end = "2010005", SDSstring = "101", outProj = "+proj=utm +zone=34 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
2. sf::gdal_utils(util = "warp", source = gdalSDS, destination = ofile, 
 .     options = c(params, "-overwrite", "-multi"), quiet = !is.null(opts$quiet) && 
 .         opts$quiet)
3. stop(paste0("gdal_utils ", util, ": an error occured"))

[fdetsch]

Hi,

can you please provide your sessionInfo()? Does MODIS:::checkHdf4Driver() evaluate to TRUE?

[caiomattos]

Hi,

Thanks for the response.

Sure! checkHdf4Driver() does evaluate to TRUE. And the sessionInfo() is below:

R version 3.6.3 (2020-02-29)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /cache/home/cr577/anaconda3/envs/amazon/lib/libopenblasp-r0.3.10.so

locale:
[1] C

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

other attached packages:
[1] MODIS_1.2.2   mapdata_2.3.0 maps_3.3.0    raster_3.3-7  rgdal_1.4-8  
[6] sp_1.4-2     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6       lattice_0.20-41    prettyunits_1.1.1  class_7.3-17      
 [5] ps_1.3.3           assertthat_0.2.1   rprojroot_1.3-2    digest_0.6.25     
 [9] mime_0.9           IRdisplay_0.7.0    R6_2.4.1           repr_1.1.0        
[13] backports_1.1.8    evaluate_0.14      e1071_1.7-3        pillar_1.4.4      
[17] rlang_0.4.6        curl_4.3           uuid_0.1-4         callr_3.4.3       
[21] nloptr_1.2.2.1     desc_1.2.0         devtools_2.3.0     foreign_0.8-71    
[25] shiny_1.5.0        compiler_3.6.3     httpuv_1.5.4       pkgconfig_2.0.3   
[29] base64enc_0.1-3    pkgbuild_1.0.8     rgeos_0.5-3        htmltools_0.5.0   
[33] tidyselect_1.1.0   mapedit_0.6.0      tibble_3.0.1       codetools_0.2-16  
[37] fansi_0.4.1        dplyr_1.0.0        crayon_1.3.4       withr_2.2.0       
[41] later_1.1.0.1      sf_0.9-3           bitops_1.0-6       grid_3.6.3        
[45] lifecycle_0.2.0    jsonlite_1.7.0     xtable_1.8-4       DBI_1.1.0         
[49] magrittr_1.5       units_0.6-7        KernSmooth_2.23-17 cli_2.0.2         
[53] ptw_1.9-15         fs_1.4.1           promises_1.1.1     remotes_2.1.1     
[57] testthat_2.3.2     generics_0.0.2     vctrs_0.3.1        ellipsis_0.3.1    
[61] IRkernel_0.8.15    tools_3.6.3        glue_1.4.1         purrr_0.3.4       
[65] processx_3.4.2     pkgload_1.1.0      parallel_3.6.3     fastmap_1.0.1     
[69] maptools_1.0-1     sessioninfo_1.1.1  classInt_0.4-3     memoise_1.1.0     
[73] pbdZMQ_0.3-3       usethis_1.6.1

[fdetsch]

Can you download an arbitrary .hdf file and then try to extract basic information from it using sf and raster functionality, e.g.

hdfs = getHdf(product="MOD11A1", extent="Austria", begin="2010001", end="2010001")
sds = sf::gdal_subdatasets(unlist(hdfs)[1]) # list subdatasets
sds

rst = raster(sds[[1]]) # read 'LST_Day_1km' sds
rst

Does this work?

[caiomattos]

Yes, that appears to work just fine, no errors.

[fdetsch]

So it seems like there's something going wrong with the raster layer options passed down to sf::gdal_utils(). Can you do me the favor and try the following:

hdfs = getHdf(product="MOD11A1", extent="Austria", begin="2010001", end="2010001")

gdal_sds = sapply(
  lapply(
    unlist(hdfs)
    , sf::gdal_subdatasets
  )
  , "[["
  , 1
)

ofile = tempfile(fileext = ".tif")
sf::gdal_utils(
  util = "warp"
  , source = gdal_sds
  , destination = ofile
  , options = c(
    "-of", "GTiff"
    , "-co", "compress=lzw", "predictor=2"
    , "-r", "near"
    , "-srcnodata", "0"
    , "-dstnodata", "0"
    # , "-t_srs", "+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs"
    , "-overwrite"
    , "-multi"
  )
  , quiet = !is.null(opts$quiet) && opts$quiet
)

rst = raster(ofile)
rst

both with and without the -t_srs part commented?

[itati01]

Hi,

I stumbled upon this issue and tested the code. The second assignment (gdal_sds) gives an indexing error. A second unlist() resolves the issue.

kk<- lapply(unlist(hdfs),sf::gdal_subdatasets) # lapply part ok
gdal_sds = sapply(kk, "[[", 1)  # sapply part error
gdal_sds = sapply(unlist(kk), "[[", 1)  # modified sapply works

gdal_utils works fine.

[fdetsch]

@caiomattos As suggested by @mlampros in #96, can you please

remotes::install_github("MatMatt/MODIS", ref = "develop")

and then try to rerun your code?