# Part 4 - Applying spatial filters

Previously, Part 3 of this guide series to `arcgis.geometry` module, you have been introduced to two ways of conducting spatial operations such as `union` and `intersection`, with `ArcGIS API for Python` - the OOP style vs calling tools off `arcgis.geometry.functions`. Now in Part 4 let us continue to explore how the spatial filtering can be applied within the `arcgis.geometry.filters` sub-module. This module includes functions to filter query results using a spatial relationship with another geometry. It is used when querying feature layers and imagery layers. The spatial filtering is even more powerful, when integrated with `geoenrichment`.

## Use Case

For example, Jim is a Highway Engineer working for the Department of Transportation, and he is performing some fact checking with interstate highways in California. Let's see how he takes advantage of the `arcgis.geometry.filters` in answering these questions:

• Among all Interstate Highways in the United States, which ones go through San Bernadino County of California, and which ones go through the Los Angeles County.
• Are there any highways that are wholly contained inside San Bernadino County? Are there any ones that only run inside the LA County?
• Check if any of the Interstate Highways are too close to the wilderness protected areas in San Bernadino County.
• In case of wildfires, are the wilderness protected areas in a safe distance from the incidents? and which Interstate Highways are too close to the areas impacted by wildfires?

## Import and Data Preparation

First of all, let us import the necessary libraries, and then create a GIS connection object to the ArcGIS online organization.

``````from arcgis.gis import GIS
from arcgis.geometry import Geometry, Polyline, Point, union
from arcgis.geometry.filters import intersects, contains, overlaps, crosses, touches, within
from arcgis.geometry.filters import envelope_intersects, index_intersects
from arcgis.geoenrichment import Country
from arcgis.features import FeatureLayer``````
``gis = GIS('home')``

### Input Data (Freeway)

``````freeway_item = gis.content.search("USA Freeway System AND type:Feature Layer", outside_org=True)[0]
freeway_item``````
USA Freeway System
This layer presents rural and urban interstate highways.Feature Layer Collection by esri_dm
``````freeway_lyr = freeway_item.layers[0]
freeway_lyr``````
`<FeatureLayer url:"https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/USA_Freeway_System/FeatureServer/1">`

### Input Data (Fire)

``````wildfire_item = gis.content.search("USA Current Wildfires AND type:Feature Layer", outside_org=True)[0]
wildfire_item``````
USA Current Wildfires
This layer shows wildfires that have been updated within the past 7 days in the United States from IRWIN and NIFC information.Feature Layer Collection by esri_livefeeds2
``````wildfire_lyr = wildfire_item.layers[0]
wildfire_lyr``````
`<FeatureLayer url:"https://services9.arcgis.com/RHVPKKiFTONKtxq3/arcgis/rest/services/USA_Wildfires_v1/FeatureServer/0">`

### Input Data (Wilderness Preservation Area)

``````# "Wilderness Areas in the United States"
wild_areas_item = gis.content.search("Wilderness Areas in the United States AND owner:wildernesswebmaster99",
item_type="Feature Layer",
outside_org=True)[0]
wild_areas_item``````
Wilderness Areas in the United States
The National Wilderness Preservation System includes wilderness areas in the United States designated by the Wilderness Act of 1964 and all subsequent wilderness laws. Today, more than 800 wilderness areas exist in 46 states and Puerto Rico encompassing 111 million acres. Through designation as wilderness, public lands that are federally owned by the Bureau of Land Management, Fish and Wildlife Service, Forest Service, or National Park Service or given the highest level of protection available.Feature Layer Collection by wildernesswebmaster99
``````wild_areas_lyr = wild_areas_item.layers[0]
wild_areas_lyr``````
`<FeatureLayer url:"https://services1.arcgis.com/ERdCHt0sNM6dENSD/arcgis/rest/services/Wilderness_Areas_in_the_United_States/FeatureServer/1">`

### State and County Boundaries from Geoenrichment

`GeoEnrichment` provides the ability to get facts about a location or area. Using GeoEnrichment, you can get information about the people, places, and businesses in a specific area or within a certain distance or drive time from a location. It enables you to query and use information from a large collection of data sets including population, income, housing, consumer behavior, and the natural environment.

Next, `arcgis.geoenrichment.Country` is used to derive the geometries of the San Bernadino County, the LA County, and the Riverside County.

``````usa = Country.get('US')
type(usa)``````
`arcgis.geoenrichment.enrichment.Country`

``````named_area_sb_county = usa.subgeographies.states['California'].counties['San_Bernardino_County']
display(named_area_sb_county)
named_area_sb_county.geometry.as_arcpy``````
`<NamedArea name:"San Bernardino County" area_id="06071", level="US.Counties", country="United States">`
``````sr_sb_county = named_area_sb_county.geometry["spatialReference"]
sr_sb_county``````
`{'wkid': 4326, 'latestWkid': 4326}`
``````named_area_la_county = usa.search(query='Los Angeles', layers=['US.Counties'])
display(named_area_la_county[0])
named_area_la_county[0].geometry.as_arcpy``````
`<NamedArea name:"Los Angeles County" area_id="06037", level="US.Counties", country="United States">`
``````named_area_riverside_county = usa.search(query='Riverside', layers=['US.Counties'])
display(named_area_riverside_county[0])
named_area_riverside_county[0].geometry.as_arcpy``````
`<NamedArea name:"Riverside County" area_id="06065", level="US.Counties", country="United States">`

## `arcgis.geometry.filters` module

Now all input data are ready, you are all set up to start creating spatial filters! The `arcgis.geometry.filters` module contains functions to filter query results of a feature or imagery layer by a spatial relationship with another geometry. You are going to see how `intersects()`, `contains`, etc. are being used to answer Jim's questions in the following sections. Before we get started, let us visualize how these layers and counties overlay spatially with a help of a map.

``````map1 = gis.map('San Bernardino, California')
map1``````
``````map1.add_layer(freeway_item) # red lines

# draw the counties
map1.draw(named_area_sb_county.geometry) # large gray polygon``````

### `intersects`

To find which parts of interstate highways layer intersect with the county of San Bernardino (or county of Los Angeles for instance), we will use `arcgis.geometry.filters.intersects(geometry, sr=None)` to create a geometry filter object that filters results whose geometry intersects with the specified geometry.

``````# construct a geometry filter using the filter geometry
sb_county_filter = intersects(named_area_sb_county.geometry,
sr=sr_sb_county)
type(sb_county_filter)``````
`dict`
``sb_county_filter.keys()``
`dict_keys(['geometry', 'geometryType', 'spatialRel', 'inSR'])`
``sb_county_filter['spatialRel']``
`'esriSpatialRelIntersects'`
``````# query a feature layer for features that meet filter criteria
df_sb_county_Interstate = freeway_lyr.query(geometry_filter=sb_county_filter,
as_df=True)
df_sb_county_Interstate[["OBJECTID","ROUTE_NUM","DIST_MILES","DIST_KM"]]``````
OBJECTIDROUTE_NUMDIST_MILESDIST_KM
04I21593.43523150.36972
142S6614.1056622.70091
247S58119.34584192.06890
348S6066.98541107.80279
4139I102367.174473809.60565
5639I402388.994333844.72138
6672I151410.539542270.04788
7690S21040.4692565.12908
``df_sb_county_Interstate["ROUTE_NUM"].tolist()``
`['I215', 'S66', 'S58', 'S60', 'I10', 'I40', 'I15', 'S210']`

The results above of an `intersects` filter on the query show 8 Interstate highways go through (or should we say... intersect) San Bernardino County, namely highways `I215, S66, S58, S60, I10, I40, I15, and S210`. If you are interested in finding out more information about any of the resulting highways, such as the specific geometry/shape, you can run the subsequent two cells to create composite filtering criteria (with `where=` and `geometry_filter=`) to get a finer query result. After that, parse the properties to get more details.

``````# query a feature layer for features that meet filter criteria
sb_county_I15 = freeway_lyr.query(where="ROUTE_NUM = 'I15'",
geometry_filter=sb_county_filter)
sb_county_I15``````
`<FeatureSet> 1 features`
``````for feature in sb_county_I15.features:
print(feature.attributes)
geom_dict = feature.geometry
geom_dict["spatialReference"] = sr_sb_county
freeway_obj = Polyline(geom_dict)
display(freeway_obj.as_arcpy)``````
```{'OBJECTID': 672, 'ROUTE_NUM': 'I15', 'CLASS': 'I', 'NUMBER': '15', 'SUFFIX': ' ', 'DIST_MILES': 1410.53954, 'DIST_KM': 2270.04788}
```

Similar to querying for highways intersecting the SB County, you can also query highways that intersect Los Angeles county.

``````# construct a geometry filter using the filter geometry
la_county_filter = intersects(named_area_la_county[0].geometry,
sr=sr_sb_county)``````
``````# query a feature layer for features that meet filter criteria
df_la_county_Interstate = freeway_lyr.query(geometry_filter=la_county_filter,
as_df=True)
display(df_la_county_Interstate.shape)
`(24, 8)`
OBJECTIDROUTE_NUMDIST_MILESDIST_KM
02I10519.8304131.91401
110I60527.1368943.67268
211I71020.0254732.22794
322S100.639531.02922
424S1109.0880214.62578

When the same highways feature layer was queried using the LA County filter, the query result had 24 entries (as seen from the shape of the DataFrame object). This means 24 interstate highways intersect LA County.

If interested in finding out more information about route `U101`, you can run the following two cells to create composite filtering criteria with `where` and `geometry_filter`, get a finer search result, and parse the detailed properties.

``````# query a feature layer for features that meet filter criteria
la_county_U101 = freeway_lyr.query(where="ROUTE_NUM = 'U101'",
geometry_filter=la_county_filter)
la_county_U101``````
`<FeatureSet> 1 features`
``````for feature in la_county_U101.features:
print(feature.attributes)
geom_dict = feature.geometry
geom_dict["spatialReference"] = sr_sb_county
freeway_obj = Polyline(geom_dict)
display(freeway_obj.as_arcpy)``````
```{'OBJECTID': 675, 'ROUTE_NUM': 'U101', 'CLASS': 'U', 'NUMBER': '101', 'SUFFIX': ' ', 'DIST_MILES': 722.19665, 'DIST_KM': 1162.26517}
```

## `contains` V.S. `within`

In the previous section, we have sorted out that 8 Interstate Highways, namely `I215, S66, S58, S60, I10, I40, I15, and S210`, intersect with county of San Bernardino, and 24 routes intersect with LA County. Now, let us find out which highways are entirely contained in these counties. You will need to use `arcgis.geometry.filters.contains(geometry, sr=None)` to compose the spatial filter that aims to return a feature if its shape is wholly contained within the search geometry (valid for all shape type combinations).

On a side note, `arcgis.geometry.filters.within(geometry, sr=None)` is the opposite to `contains` in that it returns a feature if its shape wholly contains the search geometry (valid for all shape type combinations).

``````# construct a geometry filter using the filter geometry
sb_county_filter2 = contains(named_area_sb_county.geometry,
sr=sr_sb_county)``````
``````# query a feature layer for features that meet filter criteria
df_sb_county_Interstate2 = freeway_lyr.query(geometry_filter=sb_county_filter2,
as_df=True)
display(df_sb_county_Interstate2.shape)
`(0, 8)`
OBJECTIDROUTE_NUMDIST_MILESDIST_KM

From the cell output above, none of Interstate Highways is contained within the county boundaries of SB County. What about the county of Los Angeles?

``````# construct a geometry filter using the filter geometry
la_county_filter2 = contains(named_area_la_county[0].geometry,
sr=sr_sb_county)``````
``````# query a feature layer for features that meet filter criteria
df_la_county_Interstate2 = freeway_lyr.query(geometry_filter=la_county_filter2,
as_df=True)
display(df_la_county_Interstate2.shape)
`(9, 8)`
OBJECTIDROUTE_NUMDIST_MILESDIST_KM
011I71020.0254732.22794
122S100.639531.02922
224S1109.0880214.62578
329S13412.3029119.79965
434S1706.023929.69458
``````map2 = gis.map('Los Angeles County, California')
map2``````
``````# map2.add_layer(freeway_item)
map2.draw(named_area_la_county[0].geometry)
df_la_county_Interstate2.spatial.plot(map2)``````
`True`

The map above displays the boundary of LA county and the highways that are fully contained within it. The cells below assert this claim by comparing the names of higways that intersect vs is contained withint the county boundaries.

``````num_contained = 0
for route_num in df_la_county_Interstate2["ROUTE_NUM"]:
if route_num not in df_la_county_Interstate["ROUTE_NUM"]:
num_contained+=1
print(num_contained, "out of 9 are contained in LA County, and not intersecting with County Boundaries")``````
```9 out of 9 are contained in LA County, and not intersecting with County Boundaries
```
``````num_intersecting = 0
for route_num in df_la_county_Interstate["ROUTE_NUM"]:
if route_num not in df_la_county_Interstate2["ROUTE_NUM"]:
num_intersecting+=1
print(num_intersecting, "out of 24 are intersecting with County Boundaries, but not contained in LA County")``````
```24 out of 24 are intersecting with County Boundaries, but not contained in LA County
```

### `crosses` V.S. `overlaps`

Filters `crosses` and `overlaps` share some common aspects in their definitions:

• `arcgis.geometry.filters.crosses(geometry, sr=None)` returns a feature if the intersection of the interiors of the two shapes is not empty and has a lower dimension than the maximum dimension of the two shapes. Two lines that share an endpoint in common do not cross. Valid for Line/Line, Line/Area, Multi-point/Area, and Multi-point/Line shape type combinations.
• `arcgis.geometry.filters.overlaps(geometry, sr=None)` returns a feature if the intersection of the two shapes results in an object of the same dimension, but different from both of the shapes. Applies to Area/Area, Line/Line, and Multi-point/Multi-point shape type combinations.

The major difference is that `overlaps` requires two shapes with the same dimension (e.g. Area/Area, Line/Line, and Multi-point/Multi-point shape type combinations), while `crosses` can operate on two shapes whose intersection has a lower dimension than the maximum dimension of the two.

``````# query a feature layer for features that meet filter criteria
df_sb_county_wildfire = wildfire_lyr.query(geometry_filter=sb_county_filter)
df_sb_county_wildfire``````
`<FeatureSet> 9 features`

The query above returns wildfire incidents inside SB County. Next, we can now loop through each fire incident (point feature), create a buffer around them and apply `overlaps` spatial filter to identify which wilderness preservation areas overlap with the fire buffer.

Note: The `wild_areas_lyr` layer pulls live wildfire information. Hence, when you run this notebook, you might see different set of results for the cells below.
``````fire_buffer_list = []
wild_areas_list = []

# loop through all point features - wildfire incidents
for feature in df_sb_county_wildfire.features:
dict_geom = feature.geometry
dict_geom["spatialReference"] = sr_sb_county

# create buffered area around the wildfire centers
f_buffer = Point(dict_geom).buffer(0.05).generalize(0.01)
fire_buffer_list.append(f_buffer)

# creates a spatial filter
f_buffer_filter = overlaps(f_buffer, sr=sr_sb_county)

# query a feature layer for features that meet filter criteria;
# result == wilderness prevention areas that overlap with the wildfire buffered zones
f_buffer_wild_areas = wild_areas_lyr.query(geometry_filter=f_buffer_filter)
df_f_buffer_wild_areas = f_buffer_wild_areas.sdf

# print results and save
if df_f_buffer_wild_areas.shape[0] > 0:
display(df_f_buffer_wild_areas[["OBJECTID","NAME","Land_Type"]])
wild_areas_list.append(f_buffer_wild_areas)``````
OBJECTIDNAMELand_Type
0704Piute Mountains WildernessPrivate Land
11213Piute Mountains WildernessState Government Land

#### Visualizing the filtered results

Next, let us visualize the wildfires, the buffered zone impacted by wildfires, and the impacted wilderness preservation areas on a single MapView.

``````map3 = gis.map('Redlands, CA')
map3.basemap = 'streets'
map3``````
``````# draw the fires as diamonds
map3.draw(df_sb_county_wildfire,
symbol = {"angle":0,"xoffset":0,"yoffset":0,"type":"esriPMS",
"url":"http://static.arcgis.com/images/Symbols/Shapes/RedDiamondLargeB.png",
"contentType":"image/png","width":24,"height":24})

# draw buffers around fires
for f_b in fire_buffer_list:
map3.draw(f_b)

# draw the wild areas affeted
for w_a in wild_areas_list:
map3.draw(w_a,
symbol = {"type" : "esriSFS","style": "esriSFSSolid",
"color": [115,76,0,255],
"outline": {"type": "esriSLS","style": "esriSLSSolid",
"color": [110,110,110,255],"width": 1
}
})``````

In order to decide what secondary roads and local roads might get affected by the wildfire in SB County, we need to first create a FeatureLayer object from the transportation resources provided by census bureau.

``````fl_url = "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Transportation_LargeScale/MapServer"

We can use the class `PythonSwitch` defined below to explore all possible spatial relations within `arcgis.geometry.filters`, and discover the differences of each spatial filter:

``````# Implement Python Switch Case Statement using Class
class PythonSwitch:
def switch(self, op):
default = "Unsupported Operation"
return getattr(self, 'case_' + str(op), lambda: default)()

def case_1(self):
return crosses(f_buffer, sr=sr_sb_county)

def case_2(self):
return contains(f_buffer, sr=sr_sb_county)

def case_3(self):
return intersects(f_buffer, sr=sr_sb_county)

def case_4(self):
return overlaps(f_buffer, sr=sr_sb_county)

def case_5(self):

def case_6(self):
return within(f_buffer, sr=sr_sb_county)

def case_7(self):
return envelope_intersects(f_buffer, sr=sr_sb_county)

def case_8(self):
return index_intersects(f_buffer, sr=sr_sb_county)``````
``````# For experimental purposes, just target at one of the fire buffered areas
op_list = ['crosses','contains','intersects','overlaps','touches','within','envelope_intersects','index_intersects']
for f_buffer in fire_buffer_list[1:2]:
# switch case to loop thru all spatial relations
for i in range(1,8):
# creates a spatial filter
s = PythonSwitch()
f_buffer_filter = s.switch(i)

# query a feature layer for features that meet filter criteria

# print results and save
else:
print(f"Op {i} {op_list[i]}: 0")``````
```Op 1 contains: 2
```
BASENAMEMTFCCNAMEOBJECTID
0330S1200State Rte 330148784
118S1200State Rte 18149536
```Op 2 intersects: 2
```
BASENAMEMTFCCNAMEOBJECTID
0City CreekS1200City Creek Rd105267
1HilltopS1200Hilltop Blvd149111
```Op 3 overlaps: 7
```
BASENAMEMTFCCNAMEOBJECTID
0City CreekS1200City Creek Rd105267
1330S1200State Rte 330148784
2HilltopS1200Hilltop Blvd149111
318S1200State Rte 18149536
4HilltopS1200Hilltop Con152083
5HilltopS1200Hilltop Blvd152199
6HilltopS1200Hilltop Blvd152200
```Op 4 touches: 0
Op 5 within: 0
Op 6 envelope_intersects: 0
Op 7 index_intersects: 11
```
BASENAMEMTFCCNAMEOBJECTID
0City CreekS1200City Creek Rd105267
1330S1200State Rte 330148784
2HilltopS1200Hilltop Blvd149111
318S1200State Rte 18149536
4VeS1200Ve Trl149657
5173S1200State Rte 173151713
6HilltopS1200Hilltop Con152083
7HilltopS1200Hilltop Blvd152199
8HilltopS1200Hilltop Blvd152200
938S1200State Rte 38222160
10Rim of the WorldS1200Rim of the World Hwy223822

After carefully observing the results, we can tell in this example, `intersects` works the same way as `envelope_intersects`, and `crosses` & `contains` each represent a mutually exclusive subset of what's derived from `intersects`.

``````impacted_sec_roads_list = []

for f_buffer in fire_buffer_list:
# creates a spatial filter
s = PythonSwitch()
f_buffer_filter = s.switch(3)
##f_buffer_filter = intersects(f_buffer, sr=sr_sb_county)

# query a feature layer for features that meet filter criteria

# print results and save

# query a feature layer for features that meet filter criteria

# print results and save

Now that we have obtained `impacted_sec_roads_list` and `impacted_loc_roads_list` that store the secondary and local roads impacted by the wildfire, it is straight forward to just map them in the Map Widget as well.

``````for road in impacted_sec_roads_list:
symbol = {
"type": "simple-line",
"color": "darkblue",
"width": "8px",
"style": "short-dot"
})``````
``````for road in impacted_loc_roads_list:
symbol = {
"type": "simple-line",
"color": "lightblue",
"width": "4px",
"style": "short-dot"
})``````

## Conclusion

This part 4 of the geometry module guide demonstrated how to construct spatial filter using `arcgis.goemetry.filters` sub-module. These filters are can be applied when querying feature and imagery layers.

The choice to make between various spatial relations for spatial queries depends on what is known prior to executing the query and the type of results needed. If the goal of the query is to find features that satisfy a spatial relationship with a single geometry (or a collection of geometries), use a spatial filter.

What's not mentioned here but still important for improving efficiencies and precisons of queries, is that users should consider `spatial caching` and `tolerance` when executing spatial queries. If multiple spatial queries are going to be executed within a common extent, the use of spatial caching can significantly increase performance by reducing round trips to the data source. A feature class x,y tolerance can have an effect on the results of a spatial query and should be considered when executing spatial queries, particularly with feature classes that have unusually large x,y tolerances. For more information, please check out Resouces on ArcGIS engine.