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.

Input
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
Input
gis = GIS('home')

Input Data (Freeway)

Input
freeway_item = gis.content.search("USA Freeway System AND type:Feature Layer", outside_org=True)[0]
freeway_item
Output
USA Freeway System
This layer presents rural and urban interstate highways.Feature Layer Collection by esri_dm
Last Modified: January 07, 2020
2 comments, 2,353,899 views
Input
freeway_lyr = freeway_item.layers[0]
freeway_lyr
Output
<FeatureLayer url:"https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/USA_Freeway_System/FeatureServer/1">

Input Data (Fire)

Input
wildfire_item = gis.content.search("USA Current Wildfires AND type:Feature Layer", outside_org=True)[0]
wildfire_item
Output
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
Last Modified: July 06, 2020
6 comments, 1,682,915 views
Input
wildfire_lyr = wildfire_item.layers[0]
wildfire_lyr
Output
<FeatureLayer url:"https://services9.arcgis.com/RHVPKKiFTONKtxq3/arcgis/rest/services/USA_Wildfires_v1/FeatureServer/0">

Input Data (Wilderness Preservation Area)

Input
# "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
Output
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
Last Modified: April 10, 2020
0 comments, 229,830 views
Input
wild_areas_lyr = wild_areas_item.layers[0]
wild_areas_lyr
Output
<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.

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

To learn more about GeoEnrichment, refer to the guide on this module here.

Input
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">
Output
Input
sr_sb_county = named_area_sb_county.geometry["spatialReference"]
sr_sb_county
Output
{'wkid': 4326, 'latestWkid': 4326}
Input
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">
Output
Input
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">
Output

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.

Input
map1 = gis.map('San Bernardino, California')
map1
Input
map1.add_layer(freeway_item) # red lines
map1.add_layer(wildfire_item) # dots
map1.add_layer(wild_areas_item) # white polygons

# 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.

Input
# 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)
Output
dict
Input
sb_county_filter.keys()
Output
dict_keys(['geometry', 'geometryType', 'spatialRel', 'inSR'])
Input
sb_county_filter['spatialRel']
Output
'esriSpatialRelIntersects'
Input
# 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"]]
Output
OBJECTID ROUTE_NUM DIST_MILES DIST_KM
0 4 I215 93.43523 150.36972
1 42 S66 14.10566 22.70091
2 47 S58 119.34584 192.06890
3 48 S60 66.98541 107.80279
4 139 I10 2367.17447 3809.60565
5 639 I40 2388.99433 3844.72138
6 672 I15 1410.53954 2270.04788
7 690 S210 40.46925 65.12908
Input
df_sb_county_Interstate["ROUTE_NUM"].tolist()
Output
['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.

Input
# 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
Output
<FeatureSet> 1 features
Input
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.

Input
# construct a geometry filter using the filter geometry
la_county_filter = intersects(named_area_la_county[0].geometry, 
                              sr=sr_sb_county)
Input
# 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)
df_la_county_Interstate[["OBJECTID","ROUTE_NUM","DIST_MILES","DIST_KM"]].head()
(24, 8)
Output
OBJECTID ROUTE_NUM DIST_MILES DIST_KM
0 2 I105 19.83041 31.91401
1 10 I605 27.13689 43.67268
2 11 I710 20.02547 32.22794
3 22 S10 0.63953 1.02922
4 24 S110 9.08802 14.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.

Input
# 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
Output
<FeatureSet> 1 features
Input
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).

Input
# construct a geometry filter using the filter geometry
sb_county_filter2 = contains(named_area_sb_county.geometry, 
                             sr=sr_sb_county)
Input
# 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)
df_sb_county_Interstate2[["OBJECTID","ROUTE_NUM","DIST_MILES","DIST_KM"]].head()
(0, 8)
Output
OBJECTID ROUTE_NUM DIST_MILES DIST_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?

Input
# construct a geometry filter using the filter geometry
la_county_filter2 = contains(named_area_la_county[0].geometry, 
                             sr=sr_sb_county)
Input
# 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)
df_la_county_Interstate2[["OBJECTID","ROUTE_NUM","DIST_MILES","DIST_KM"]].head()
(9, 8)
Output
OBJECTID ROUTE_NUM DIST_MILES DIST_KM
0 11 I710 20.02547 32.22794
1 22 S10 0.63953 1.02922
2 24 S110 9.08802 14.62578
3 29 S134 12.30291 19.79965
4 34 S170 6.02392 9.69458
Input
map2 = gis.map('Los Angeles County, California')
map2
Input
# map2.add_layer(freeway_item)
map2.draw(named_area_la_county[0].geometry)
df_la_county_Interstate2.spatial.plot(map2)
Output
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.

Input
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
Input
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.

Input
# 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
Output
<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.
Input
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)
OBJECTID NAME Land_Type
0 704 Piute Mountains Wilderness Private Land
1 1213 Piute Mountains Wilderness State 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.

Input
map3 = gis.map('Redlands, CA')
map3.basemap = 'streets'
map3
Input
# 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
                          }
                      })

Consider the impacted roads

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.

Input
fl_url = "https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/Transportation_LargeScale/MapServer"
secondary_roads_layer = FeatureLayer(fl_url+"/1")
local_roads_layer = FeatureLayer(fl_url+"/2")

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:

Input
# 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):
        return touches(f_buffer, sr=sr_sb_county)
 
    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)
Input
# 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
        f_buffer_sec_roads = secondary_roads_layer.query(geometry_filter=f_buffer_filter)
        df_f_buffer_sec_roads = f_buffer_sec_roads.sdf

        # print results and save
        if df_f_buffer_sec_roads.shape[0] > 0:
            print(f"Op {i} {op_list[i]}: {df_f_buffer_sec_roads.shape[0]}")
            display(df_f_buffer_sec_roads[["BASENAME","MTFCC","NAME","OBJECTID"]])
        else:
            print(f"Op {i} {op_list[i]}: 0")
Op 1 contains: 2
BASENAME MTFCC NAME OBJECTID
0 330 S1200 State Rte 330 148784
1 18 S1200 State Rte 18 149536
Op 2 intersects: 2
BASENAME MTFCC NAME OBJECTID
0 City Creek S1200 City Creek Rd 105267
1 Hilltop S1200 Hilltop Blvd 149111
Op 3 overlaps: 7
BASENAME MTFCC NAME OBJECTID
0 City Creek S1200 City Creek Rd 105267
1 330 S1200 State Rte 330 148784
2 Hilltop S1200 Hilltop Blvd 149111
3 18 S1200 State Rte 18 149536
4 Hilltop S1200 Hilltop Con 152083
5 Hilltop S1200 Hilltop Blvd 152199
6 Hilltop S1200 Hilltop Blvd 152200
Op 4 touches: 0
Op 5 within: 0
Op 6 envelope_intersects: 0
Op 7 index_intersects: 11
BASENAME MTFCC NAME OBJECTID
0 City Creek S1200 City Creek Rd 105267
1 330 S1200 State Rte 330 148784
2 Hilltop S1200 Hilltop Blvd 149111
3 18 S1200 State Rte 18 149536
4 Ve S1200 Ve Trl 149657
5 173 S1200 State Rte 173 151713
6 Hilltop S1200 Hilltop Con 152083
7 Hilltop S1200 Hilltop Blvd 152199
8 Hilltop S1200 Hilltop Blvd 152200
9 38 S1200 State Rte 38 222160
10 Rim of the World S1200 Rim of the World Hwy 223822

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.

Input
impacted_sec_roads_list = []
impacted_loc_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
    f_buffer_sec_roads = secondary_roads_layer.query(geometry_filter=f_buffer_filter)
    df_f_buffer_sec_roads = f_buffer_sec_roads.sdf
    
    # print results and save
    if df_f_buffer_sec_roads.shape[0] > 0:
        # display(df_f_buffer_sec_roads.shape[0])
        # display(df_f_buffer_sec_roads[["BASENAME","MTFCC","NAME","OBJECTID"]].head())
        impacted_sec_roads_list.append(f_buffer_sec_roads)
        
    # query a feature layer for features that meet filter criteria
    f_buffer_loc_roads = local_roads_layer.query(geometry_filter=f_buffer_filter)
    df_f_buffer_loc_roads = f_buffer_loc_roads.sdf
    
    # print results and save
    if df_f_buffer_loc_roads.shape[0] > 0:
        # display(df_f_buffer_loc_roads.shape[0])
        # display(df_f_buffer_loc_roads[["BASENAME","MTFCC","NAME","OBJECTID"]].head())
        impacted_loc_roads_list.append(f_buffer_loc_roads)

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.

Input
for road in impacted_sec_roads_list:
    map1.draw(road,
              symbol = {
                          "type": "simple-line",
                          "color": "darkblue",
                          "width": "8px",
                          "style": "short-dot"
                        })
Input
for road in impacted_loc_roads_list:
    map1.draw(road,
              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.

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.