Home > Blockchain >  check which polygons intersect with raster
check which polygons intersect with raster

Time:08-08

I have a raster with following extent

library(terra)    
library(rnaturalearthdata)

ext(my_raster)
SpatExtent : 4.99958333333333, 9.99958333333333, -0.000416666666666667, 4.99958333333333 (xmin, xmax, ymin, ymax)

ne_shp <- ne_countries()

I want to find which countries in ne_shp does my raster intersect with

terra::intersect(ext(my_raster), ext(ne_shp))

But this returns me only the extent that is intersecting. Is there any way I can get the list of countries in the shapefile which interesct with my tile? I could do a loop and check for condition for every single country in the shapefile. For e.g.

tmp_ls <- list()
for(i in 1:nrow(ne_shp)){
  
    poly_id <- ne_shp[i, ]
    
    if(!is.null(terra::intersect(ext(my_raster), ext(poly_id)))){
      
      tmp_ls[[i]] <- poly_id$admin
    }
}

But then I have so many raster tiles for which I have to check this condition.

CodePudding user response:

Try intersecting the extent of your raster with your actual country features instead of with their (global) extent:

ints <- terra::intersect(ne_shp, ext(my_raster))

ints
#>  class       : SpatVector 
#>  geometry    : polygons 
#>  dimensions  : 5, 168  (geometries, attributes)
#>  extent      : 5.434452, 10, 0, 5  (xmin, xmax, ymin, ymax)
#>  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#>  names       :      featurecla scalerank LABELRANK SOVEREIGNT SOV_A3 ADM0_DIF
#>  type        :           <chr>     <int>     <int>      <chr>  <chr>    <int>
#>  values      : Admin-0 country         0         2    Nigeria    NGA        0
#>                Admin-0 country         0         3   Cameroon    CMR        0
#>                Admin-0 country         0         4      Gabon    GAB        0
#>  LEVEL              TYPE   TLC    ADMIN (and 158 more)
#>  <int>             <chr> <chr>    <chr>               
#>      2 Sovereign country     1  Nigeria               
#>      2 Sovereign country     1 Cameroon               
#>      2 Sovereign country     1    Gabon

plot(ints)

CodePudding user response:

You can use terra::relate

Example data

library(rnaturalearth)
library(terra)
v <- ne_countries() |> vect()
r <- rast(ext(4.99958333333333, 9.99958333333333, -0.000416666666666667, 4.99958333333333))

Solution

x <- relate(v, ext(r), "intersects")
data.frame(v)[which(x), "admin"]
#[1] "Cameroon"          "Gabon"             "Equatorial Guinea" "Nigeria" 
  • Related