Filter exposure layer by raster extents

Community
1

Hi,

I want to filter a large, national exposure dataset by the extents of regional hazard rasters, so that RiskScape doesn’t check the exposure layer on one side of the country when the hazard is at the other end. What would be the most efficient way to do that?

Thank you!
Luisa

2

Hi Luisa,

It depends somewhat on what your model is doing. RiskScape is generally pretty quick at checking if a given point is within a GeoTIFF’s bounds or not. However, if you’re dealing with large exposure features like linestring or polygon geometry, it becomes a lot more work.

Generally you’ll want to take the bounds of the GeoTIFF and check the exposure falls within those bounds. You can get the bounds with:

bounds(bookmark('region-specific.tif'))

If both features are in the same CRS, you can check for an intersection (true/false) with:

intersects(exposure, bounds(bookmark('region-specific.tif')))

This works best as the ‘on’ condition in a join step. However, if the layers are in different CRSs then you’ll have to reproject the geometry first.

If you have a large CSV full of GeoTIFFs then you can turn the bounds of all the GeoTIFFs into a coverage lookup, and then sample that to find matching regions/GeoTIFFs.

Hope that helps.

Tim

3

Hi Tim,
Do you mean something more or less like this?

project.ini

[project]
description = Clips an input geopackage down to the extent of a raster

[model clipInputGpkg]
framework = pipeline
location = clipInputGpkg.txt

[bookmark inputGpkg]
location = input.gpkg
format = geopackage

[bookmark clipRaster]
location = raster.tif
format = geotiff

clipInputGpkg.txt

input(relation: 'inputGpkg') as exposure
-> select({*, intersects(geom, bounds(bookmark('clipRaster')))})
-> filter("intersects")
-> save(name: 'output', format: 'geopackage')

What would be the advantage of doing it in a join step rather than filter?

4

Yes, that looks about right.

You technically don’t need the bookmarks in this case - you could just use the .gpkg/.tif files directly in the pipeline (or as model parameters).

You wouldn’t need a join step in your example. Luisa’s original question made it sound like she was processing multiple raster files in a single model, in which case it would probably already have a join step where the exposure-layer is joined to multiple hazard-layers (i.e. like a probabilistic model). It’s probably a tiny bit more efficient to do the bounds check in the join step, but there wouldn’t be a huge difference compared to using a filter step.

5

So if the clipRaster bookmark was a csv of raster files how would that change my pipeline? Is there any way I can check against just one of the rasters if I know that they all have the same extent?

6

Actually there’s a limitation there currently where the bounds() function won’t work with a dynamic coverage (i.e. raster data loaded from a CSV file).

It’s a bit fiddly, but as a work-around you could do something like this:

  1. Run this pipeline, which will take your CSV of rasters and turn it into pipeline code (but as a CSV output).
riskscape pipeline eval -o . "
  input('ABunchOfRasters.csv')
  -> select({ 'input(value: { bounds(bookmark(\'' + file + '\')) as bounds, }) -> union' as pipeline }) as pipeline
"
  1. Take the pipeline.csv that got generated and copy-paste the contents into a new-pipeline.txt file. Skip the first line and use find/replace to delete the " characters. If you’re a Linux user, you could just do:
tail -n +2  pipeline.csv | sed s/\"//g > new-pipeline.txt
  1. Add the following to the end of new-pipeline.txt:
union()
-> group(by: bounds, select: { * })
-> save('hazard-extent', format: 'geopackage')
  1. Run the pipeline, e.g. riskscape pipeline eval new-pipeline.txt. It should produce a hazard-extent.gpkg result.

The hazard-extent.gpkg should contain polygons that represent the extent of all your hazard maps. You can use this in your regular pipeline (instead of ‘clipRaster’) to check which exposures intersect a hazard-layer, e.g.

filter(is_not_null(sample_one(exposure, to_coverage(bookmark('hazard-extent.gpkg')))))

Hope that makes sense.

7

hey Tim

As all my rasters have the same extent I am getting the filename of the first one (from the first column of my csv) and using that as the clipping layer

[bookmark clipRaster]
location = changeme

input(relation: 'inputGpkg') as exposure
-> select({*, intersects(geom, bounds(bookmark('clipRaster', {location: $clipRasterLocation})))})
-> filter("intersects")
-> save(name: 'output', format: 'geopackage')

to run

clipRaster=`awk -F"," 'NR==2 { print $1 }' data/Westport.csv`
riskscape model run clipInputGpkg -p"clipRasterLocation='"$clipRaster"'"

Would be cool to have the functionality to create an extent from a csv of rasters though…

Thanks for your help!