Part 4 - Applying spatial filters¶

Use Case
Import and Data Preparation
arcgis.geometry.filters module
- intersects
contains V.S. within
- crosses V.S. overlaps
  - Visualizing the filtered results
  - Consider the impacted roads
Conclusion

Previously, in 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.

In [92]:

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

In [2]:

gis = GIS('home')

Input Data (Freeway)¶

In [3]:

freeway_item = gis.content.search("USA Freeway System AND type:Feature Layer", outside_org=True)[0]
freeway_item

Out[3]:

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

In [4]:

freeway_lyr = freeway_item.layers[0]
freeway_lyr

Out[4]:

<FeatureLayer url:"https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services/USA_Freeway_System/FeatureServer/1">

Input Data (Fire)¶

In [5]:

wildfire_item = gis.content.search("USA Current Wildfires AND type:Feature Layer", outside_org=True)[0]
wildfire_item

Out[5]:

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

In [6]:

wildfire_lyr = wildfire_item.layers[0]
wildfire_lyr

Out[6]:

<FeatureLayer url:"https://services9.arcgis.com/RHVPKKiFTONKtxq3/arcgis/rest/services/USA_Wildfires_v1/FeatureServer/0">

Input Data (Wilderness Preservation Area)¶

In [7]:

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

Out[7]:

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

In [8]:

wild_areas_lyr = wild_areas_item.layers[0]
wild_areas_lyr

Out[8]:

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

In [9]:

usa = Country.get('US')
type(usa)

Out[9]:

arcgis.geoenrichment.enrichment.Country

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

In [12]:

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">

Out[12]:

In [13]:

sr_sb_county = named_area_sb_county.geometry["spatialReference"]
sr_sb_county

Out[13]:

{'wkid': 4326, 'latestWkid': 4326}

In [15]:

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">

Out[15]:

In [16]:

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">

Out[16]:

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

In [26]:

map1 = gis.map('San Bernardino, California')
map1

ArcGIS API for Python Guide