0

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.

  • Are you trying to identify rasters that are completely within the border or ones where any part is within? – Jacob F Sep 27 '18 at 14:08
  • The one where any part is whitin. I don't want the one entierly outside of the border. – C.Serrano Sep 27 '18 at 14:17
0

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)
    }
}
  • Why are you using spTransform? If there is a concern that the projection strings not matching exactly you can just use proj4string(x) <- proj4string(y). It would also seem a bit simpler to just use something like rgeos::gContains to test the Boolean as this is exactly what the function is intended to do and does not intersect all of the feature geometry and attributes like raster::intersect does (which uses rgeos::gIntersect under the hood). BTW, you can suppress warnings using the global setting: options(warn=-1) but I am not sure why this is necessary. – Jeffrey Evans Sep 27 '18 at 15:16
0

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)]
  • Thanks but this method seems to be working only for the intersected glaciers. I also wanted the one inside the limits. – C.Serrano Oct 2 '18 at 9:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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