Select multiple rasters from location

Question

I have a lot of rasters (almost 4000) wich are representing individual glaciers in the Alps. I would like to select only the rasters located within the french border.

As I have no information in the name of the file, I would like to know how I could select them without checking them one by one.

Is there any way of loading only some rasters based on there coordinates or delete rasters that don't fit a mask ?

I have looked for such a tool in QGIS and R with no sucess.

Answer 1

You can use pyqgis to check if each raster extent is within a polygon:

import os
rasterfolder = '/home/bera/GIS/DEM50m/grid50m/'

borderlyr = QgsProject.instance().mapLayersByName("New scratch layer")[0]
bordergeom = [f.geometry() for f in borderlyr.getFeatures()][0] #If you have more than one polygon feature you need to change the code.

for item in os.listdir(rasterfolder):
    if os.path.isfile(os.path.join(rasterfolder, item)) and item.endswith('.tif'):
        rl = QgsRasterLayer(os.path.join(rasterfolder, item),item.split('.')[0])
        if QgsGeometry.fromRect(rl.extent()).within(bordergeom):
            #print(item)
            QgsProject.instance().addMapLayer(rl)

Answer 2

I would set up a loop that reads each raster using the raster::raster function. This creates a pointer so, in this regard, is very fast. You can then use the raster::extent function to create a polygon representing the footprint of the given (i) raster. Using as(extent(i), "SpatialPolygons") will create an actual sp class polygon object. Once you have this polygon you can use rgeos::gOverlaps to test the intersection of the raster polygon with your jurisdictional polygon (eg., France). You can collect these results (image name, TRUE/FALSE) in a data.frame.

Here is an example that uses lapply in lieu of the for loop. The if check makes sure that the test polygon is in the same projection space as the first raster. Just ignore the warnings. Since we are creating an extent polygon on the fly, it does not have an assigned proj4string and so throws a warning regarding non-matching projection strings. I added a warning suppression but, it is not really necessary.

library(raster)
library(rgeos)
library(rgdal)

setwd(..)
r.list <- list.files(getwd(), "tif$")
p <- readOGR(getwd(), "france")
r <- raster(r.list[1])
  if(proj4string(p) != proj4string(r)) p <- spTransform(p, proj4string(r))

options(warn=-1) 
i <- lapply(r.list, FUN = function(x) { rgeos::gOverlaps( as(extent(raster(x)), 
            "SpatialPolygons"), p) })

( r.intersect <- data.frame(raster = r.list, intersects = unlist(i)) )

Now you can subset the rasters that do match.

( rmatch <- r.intersect[r.intersect$intersects == TRUE,]$raster )

If all you are after is the string of rasters that match, you can shortcut creating a data.frame and just use the index of TRUE matches and subset your raster list vector object directly.

r.list[which(unlist(lapply(r.list, FUN = function(x) { rgeos::gOverlaps( as(extent(raster(x)), 
        "SpatialPolygons"), p) })) == TRUE)]

Answer 3

Here is a solution in R using the raster and rgdal packages and a shapefile as the boundary to check against. It will add the file name of rasters that intersect the shapefile to the newlst vector.

library(raster)
library(rgdal)

# Load shapefile with boundary, change path to correct
# Source - http://biogeo.ucdavis.edu/data/diva/adm/FRA_adm.zip
shp <- shapefile('/Path/To/FRA_adm0.shp')

# Change path to folder containing rasters
rasdir <- '/Path/To/Rasters'

# List all GeoTIFF files in folder, change extension in pattern if different format
fllst <- list.files(path=rasdir, pattern='[.]tif$', full.names=T)

# New vector for storing file names of intersecting rasters
newlst <- c()

# Loop through files
for (fl in fllst){
    r <- raster(fl)
    # Transform shapefile to match crs of raster
    shp <- spTransform(shp, crs(r))
    # Check if raster intersects shapefile
    # Suppress warnings from function is optional
    if (suppressWarnings(!(is.null(intersect(shp, extent(r))))))
    {
        # If raster intersects, add file name to vector
        newlst <- c(newlst, fl)
    }
}