Using intersects with 2 polygons

Community
1

I am trying to use intersects to get all the polygons from one geopackage that intersect with a polygon from another geopackage using this pipeline

input(relation: 'catchments') as exposure
-> select({*, intersects(geom, bookmark('RegionalCouncil', {filter: $region_filter}))})
-> filter("intersects")
-> save(name: 'output', format: 'geopackage')

But I am getting this error

Failed to validate model for execution
  - Failed to validate 'select({*, intersects(geom, bookmark('RegionalCouncil', {filter: 'RegionalCouncil = \'Manawatu-Whanganui Region\''}))})' step for execution
    - Failed to validate expression '{*, intersects(geom, bookmark('RegionalCouncil', {filter: 'RegionalCouncil = \'Manawatu-Whanganui Region\''}))}' against input type {geom=>Polygon[crs=EPSG:2193], name=>Text, terms=>Text, rec2term=>Text, tocs=>Text, eva_sites=>Text, eva_rec1s=>Text, eva_rec2s=>Text, flood_stats_sites=>Text, flood_stats_rec1s=>Text}
      - Problems found with 'intersects' riskscape function
        - Type mismatch for 'rhs' function argument. Expected 'Geom' but found 'WithMetadata[Relation[{geom=>MultiPolygon[crs=EPSG:2193], REGC2021_V=>Text, REGC2021_1=>Text, REGC2021_2=>Text, LAND_AREA_=>Floating, AREA_SQ_KM=>Floating, RegionalCouncil=>Text}][crs=EPSG:2193]]'

How can I get the geometry out of a WithMetadata object?

2

Are you trying to join them up, i.e find the regional council that intersects the catchment? Easiest way is usually to sample the area layer, like you would with a hazard layer. If you run through the wizard it’ll do that for you when adding an area layer.

Cheers,
Nick

3

Like Nick suggested, sampling is probably the simplest approach, but it depends a bit on what you’re trying to do here. If you just want to filter by a region, then you can use sample_one() to find the corresponding region, e.g.

input(relation: 'catchments') as exposure
-> select({ *, sample_one(geom, to_coverage(bookmark('RegionalCouncil', {filter: $region_filter}))) as region })
-> filter(is_not_null(region))
-> save(name: 'output', format: 'geopackage')

It sounds like the RegionalCouncil layer is getting filtered so it only has one polygon, so sample_one() will be the simplest approach. If it contained multiple polygons, then approaches are:

  • Use sample() instead of sample_one(). This returns all the polygons that intersect in a list. In filtering terms, it makes it a little bit awkward to check for a region match, but this approach can be handy in other situations.
  • Use layer_intersections() instead of intersects() in a pipeline similar to what you’ve got. This would cut the geom into pieces, based on where they intersect the RegionalCouncil layer.
4

Hi Nick and Tim, thanks for the help! Tims pipeline worked perfectly. The final command to run it was

riskscape model run groupCatchmentsByRegion --param "region_filter = 'RegionalCouncil = \'Auckland Region\''"