Reading multiple MODIS LST files (MOD11A2 product) in HDF format with rast() fails, file does not exist

[dimfalk]

Messing around with some MODIS LST data (MOD11A2 acquired via USGS EarthExpolorer to be precise) in order to test tapp(), I spend some time now understandig why my workflow fails and it seems like rast() would not accept more than two HDF-EOS files before something breaks:

library(terra)
#> terra 1.7.71

# MOD11A2, tile h18v03, year 2023 (=46 files)
fnames <- list.files(pattern = "hdf$")
head(fnames, 3)
#> [1] "MOD11A2.A2023001.h18v03.061.2023015004633.hdf" 
#> [2] "MOD11A2.A2023009.h18v03.061.2023020002349.hdf" 
#> [3] "MOD11A2.A2023017.h18v03.061.2023027004006.hdf"

# 1 file = 12 layers, looking ok
r1 <- rast(fnames[1])
r1 * 1
#> class       : SpatRaster 
#> dimensions  : 1200, 1200, 12  (nrow, ncol, nlyr)
#> resolution  : 926.6254, 926.6254  (x, y)
#> extent      : 0, 1111951, 5559753, 6671703  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#> source(s)   : memory
#> varname     : MOD11A2.A2023001.h18v03.061.2023015004633 
#> names       : LST_Day_1km, QC_Day, Day_v~_time, Day_v~_angl, LST_N~t_1km, QC_Night, ... 
#> min values  :      242.40,      2,         9.7,         -65,      247.56,        2, ... 
#> max values  :      284.94,    177,        12.6,          65,      282.54,      177, ... 

# 2 files = 24 layers, still ok
r2 <- rast(fnames[1:2])
r2 * 1
#> class       : SpatRaster 
#> dimensions  : 1200, 1200, 24  (nrow, ncol, nlyr)
#> resolution  : 926.6254, 926.6254  (x, y)
#> extent      : 0, 1111951, 5559753, 6671703  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#> source(s)   : memory
#> names       : LST_Day_1km, QC_Day, Day_v~_time, Day_v~_angl, LST_N~t_1km, QC_Night, ... 
#> min values  :      242.40,      2,         9.7,         -65,      247.56,        2, ... 
#> max values  :      284.94,    177,        12.6,          65,      282.54,      177, ... 

# all of a sudden, file not found?
r3 <- rast(fnames[1:3])
r3 * 1
#> Error: [*] file does not exist: HDF4_EOS:EOS_GRID:"MOD11A2.A2023017.h18v03.061.2023027004006.hdf":MODIS_Grid_8Day_1km_LST:Emis_31

# this worked before, now it does not
r1 * 1
#> Error: [*] file does not exist: HDF4_EOS:EOS_GRID:"MOD11A2.A2023001.h18v03.061.2023015004633.hdf":MODIS_Grid_8Day_1km_LST:LST_Day_1km

Obviously, you would need this data to reproduce the issue, I guess... But I also don't want to violate any guidelines by just uploading this here.

Do you have access to Earth Explorer (any other tile would work too, probably)? Or is there any other way I can hand over these files?

sessionInfo()
#> R version 4.3.3 (2024-02-29 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)

[dimfalk]

Just wanted to check on this - and the issue still exists using R 4.4.2 and terra 1.7.83.

I realize that his is most likely some file- or format-specific issue, so I decided to provide some files for testing here.

After failing at r3 * 1, rast(fnames[1]) - this worked before, see above - additionally provides the following output:

Error: [rast] cannot open this file as a SpatRaster: 
MOD11A2.A2023001.h18v03.061.2023015004633.hdf
Additionally: Warning:
Failed to open HDF-EOS file "MOD11A2.A2023001.h18v03.061.2023015004633.hdf" for swath reading. (GDAL error 4) 

[rhijmans]

Sorry for taking so long to respond to this. When I looked at it the first time I could not figure it out, but now it is clear that this is because the HDF4 library has a very restrictive limitation of a maximum number of 32 connections (and each sub-dataset is one). That is documented here. It can be changed if you compile yourself.

Please note, that NCSA HDF library compiled with several defaults which is defined in hlimits.h file. For example, hlimits.h defines the maximum number of opened files:

#   define MAX_FILE   32

If you need open more HDF4 files simultaneously you should change this value and rebuild HDF4 library (and relink GDAL if using static HDF libraries).

I improved the handling of this situation as illustrated below. Here using the luna package to get your example files (remotes::install_github("rspatial/luna")); you need an account username/password.

uname <- "*****"
pword <- "*****"
library(terra)
f <- luna::getNASA("MOD11A2", "2023-01-01", "2023-01-18", download=TRUE, 
                          aoi=c(10,11,55,56), path=".", username=uname, password=pword)
f
#[1] "./MOD11A2.A2023001.h18v03.061.2023015004633.hdf" "./MOD11A2.A2023009.h18v03.061.2023020002349.hdf"
#[3] "./MOD11A2.A2023017.h18v03.061.2023027004006.hdf"

The below still fails, but now with a better error message:

x <- rast(f)
### each layer is a separate data source
nlyr(x)
#[1] 36
length(sources(x))
#[1] 36

y <- x * 1
#Error: [*] cannot read from HDF4_EOS:EOS_GRID:"MOD11A2.A2023017.h18v03.061.2023027004006.hdf":MODIS_Grid_8Day_1km_LST:Emis_31
#(Only 32 open datasets allowed with HDF4)

Also, the 32 connections that were opened are now closed when this error occurs, so you can continue, instead of having to restart a new R session).

One approach would be to select the sub-dataset of interest such that you can at least use 32 of these.

x <- rast(f, "LST_Day_1km")
x * 1
#class       : SpatRaster 
#dimensions  : 1200, 1200, 3  (nrow, ncol, nlyr)
#resolution  : 926.6254, 926.6254  (x, y)
#extent      : 0, 1111951, 5559753, 6671703  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
#source(s)   : memory
#names       : LST_Day_1km, LST_Day_1km, LST_Day_1km 
#min values  :      242.40,      241.82,      238.54 
#max values  :      284.94,      284.36,      279.02 

[dimfalk]

Thanks for the clarification, Robert - and your insights explaining this behaviour as well as adding a comprehensible error message!

Also, thank you for pointing out this possible workaround 👍