0

I would like to first explain what I aim to do. I am used to use the R package lidR and I find very useful the comand readLAScatalog, which creates a catalog, "a representation of a collection of las/laz files. A computer cannot load all the data at once". This is an example: enter image description here It can be used as a sf object that shows the extent of the lidar files and stores the paths to the actual file that each tile represents. I love it because it is really fast and you can select files according to location withput having to actually reading the file.

Is there a way of creating such thing from a list of raster files with terra?

EDIT: I made a loop that reads the raster files, generates a polygon out of the extend of the rasters and merge them:

require(tictoc)
require(terra)

f = list.files("PATH")

tic()
for( i in 1:length(f)){
  r = rast(f[i])
  p = as.polygons(ext(r), crs = crs(r))
  p$path = f[i]
  
  if(i == 1){ r_catalog = p}
  else{ r_catalog = rbind(r_catalog, p) }
}
toc()
plot(r_catalog)

In my case it took ~5 min to loop through 956 tif files.

enter image description here

7
  • 1
    So your desired output is? A data frame of extent and path for all .tif files in a folder? Terra is smart enough to not read all of large GeoTIFFs into memory so you can get that info without filling memory. Have you tried?
    – Spacedman
    Sep 23, 2022 at 7:19
  • I'd like a multipolygon object that just shows the extent of the rasters and one attribute would be the path to the file. I have +900 files and only need about 20 at at time. I would use this multipolygon object to know wich raster files actually intercept with my sample plots Sep 23, 2022 at 7:35
  • Have you tried looping over the files and "reading" them using terra::rast to get the metadata? If the files are large terra will only read the metadata, if the files are small it won't take long to read anyway. Alternatively have you looked at the gdalUtils package which might have functions that only get the metadata...
    – Spacedman
    Sep 23, 2022 at 7:42
  • I have and it works but I was hoping there is a fast way of doing it. I have modified my question with a solution I found. Probably there is a better way. Sep 23, 2022 at 8:48
  • That's pretty much the first way I would have tried (but I'd write it as a function of a path for re-use and testing). You can use gdalUtils::gdalinfo(tif, raw_output=FALSE)$bbox to get a bounding box - that might be faster. If I can create or find a folder full of TIFFs I'll do a comparison.
    – Spacedman
    Sep 23, 2022 at 9:08

1 Answer 1

2

Here's a version of yours written as a function without the for loop, using an lapply to loop over file names and do.call(rbind,...) to combine the list. It might be nanoseconds faster than a for loop, or nanoseconds slower. But I prefer this style as a bit cleaner than a for-loop for constructing things of unknown length:

tiffcat <- function(f){
    do.call(rbind,
            lapply(f, function(r){
                ras = rast(r)
                poly = as.polygons(ext(ras), crs=crs(ras))
                poly$path = r
                poly
            })
            )
}

Sample usage:

files = list.files(path, "*.tif", full=TRUE)
tc = tiffcat(files)

I'm not sure how to make this much faster. The help for rast says:

 When a SpatRaster is created from a file, it does not load the
 cell (pixel) values into memory (RAM). It only reads the
 parameters that describe the geometry of the SpatRaster, such as
 the number of rows and columns and the coordinate reference
 system. The actual values will be read when needed.

Certainly using gdalinfo from gdalUtils seems much much slower, probably due to a GDAL startup cost on every call, which only has to happen once when R calls GDAL from code in a package.

Another possibility is to use gdaltindex, which writes a shapefile of raster extent polygons:

> gdaltindex("index.shp", files)
NULL
> index = st_read("index.shp", quiet=TRUE)
> index
Simple feature collection with 56 features and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 0 ymin: 0 xmax: 7e+05 ymax: 1300000
CRS:           NA
First 10 features:
                               location                       geometry
1  ./1 250 000 Scale Raster/data/HP.tif POLYGON ((4e+05 1300000, 5e...
2  ./1 250 000 Scale Raster/data/HT.tif POLYGON ((3e+05 1200000, 4e...
3  ./1 250 000 Scale Raster/data/HU.tif POLYGON ((4e+05 1200000, 5e...
4  ./1 250 000 Scale Raster/data/HW.tif POLYGON ((1e+05 1100000, 2e...
5  ./1 250 000 Scale Raster/data/HX.tif POLYGON ((2e+05 1100000, 3e...

This is much faster than gdalinfo and possibly faster than reading the metadata via terra::rast since its all GDAL code. If you want to write a persistent index then this gives it to you in the shapefile. gdaltindex will also work in an "update" mode where it skips if the file is already in the index.

1
  • I tried this way as it is in fact a bit faster (15 seconds faster). I guess that there is not much room for improvement and less than 5 minutes for this amount of rasters is not bad. I also prefer this way of coding but I have not got around lapply that much. I should hough... Sep 23, 2022 at 11:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.