Monitoring hydrologic water quality in pasturelands through spatial overlay analysis

Introduction

Much of the grazing in the state of Oregon occurs on federal lands. Grazing areas are divided into allotments by federal agencies. They issue permits or leases to ranchers for individual grazing allotments.

Studies by the state of Oregon Department of Environment Quality (DEQ) indicate that streams located in grazing areas are majorly polluted by sediments and animal waste. This is a substantial concern as it causes degradation of water quality, threatening human and ecological health. The department thus wants to inspect the effect of livestock grazing on the state’s water quality.

While federal agencies manage the grazing lands by allotments, the state's biologists monitor water quality by watersheds, or hydrologic basins (as the hydrologists refer to them). If a basin has water quality issues, then biologists who monitor water quality for watersheds or hydrologic basins could identify all grazing allotments that are in that basin. They can then work with federal agencies who manages the grazing allotments to ensure that permit holders are conforming to best practices.

Grazing allotments and hydrologic basin boundaries. Many allotments fall in more than one hydrologic basin.

Since grazing allotments were not created with basin boundaries in mind, an allotment can fall completely within a hydrologic basin, or can cross basin boundaries, falling in two or more basins.

This sample uses ArcGIS API for Python to find out which watershed, or watersheds, each grazing allotment falls in, for water quality monitoring.

It demonstrates using tools such as overlay_layers to identify allotments in a particular basin. This will assign each allotment to the hydrologic basins it falls within.

Moreover, in order to successfully identify the source of pollution in each basin, each basin is assigned the grazing allotment name and the number of streams within. This will help identify which allotment each segment of each stream passes through. If field tests find a water quality issue with a particular stream, biologists can link back to the federal database and get a report on each suspect allotment (the type and number of livestock, the owner information, the administrating office, and so on). The information will help them determine the source of pollution.

Workflow

Necessary Imports

Input
%matplotlib inline

import pandas as pd
pd.set_option('mode.chained_assignment', None)
import matplotlib.pyplot as plt
from IPython.display import display
from datetime import datetime as dt

from arcgis.gis import GIS
from arcgis.features.manage_data import overlay_layers

Connect to your GIS

Connect to our GIS via an existing profile or creating a new connection by e.g. gis = GIS("https://www.arcgis.com", "arcgis_python", "P@ssword123")

Input
gis = GIS('home')

Get data for analysis

Search for grazing allotments and watersheds layer in ArcGIS Online.

Input
items = gis.content.search('grazing allotments and watersheds owner: api_data_owner', 
                           'Feature Layer')

The code below displays the items.

Input
for item in items:
    display(item)
grazing allotments and watersheds
Overlay Layers Use Case - grazing allotments and watershedsFeature Layer Collection by api_data_owner
Last Modified: May 05, 2019
0 comments, 3 views

We will use the first item for our analysis. Since the item is a Feature Layer Collection, accessing the layers property will give us a list of FeatureLayer objects.

Input
data = items[0]

The code below cycles through the layers and prints their names.

Input
for lyr in data.layers:
    print(lyr.properties.name)
hydro_units
counties
grazing_allotments
streams

Let us now get the layers and assign a variable to them.

Input
hydro_units = data.layers[0]
Input
counties = data.layers[1]
Input
grazing_allotments = data.layers[2]
Input
streams = data.layers[3]
Input
map1 = gis.map("Oregon")
map1
Output
Input
map1.add_layer(hydro_units)
map1.add_layer(grazing_allotments)

Assigning basin information to allotments

In order to find out which hydrologic basins each grazing allotment is in, we will use overlay_layers tool. It combines two layers, an analysis layer and an overlay layer, into a new layer, creating new features and combining attributes from the two layers according to the overlay method selected. Overlay operations supported are Intersect, Union, and Erase.

We will overlay grazing allotments with hydrologic basins using the Intersect option.

Input
basin_overlay = overlay_layers(grazing_allotments,
                               hydro_units,
                               overlay_type='Intersect',
                               output_name="OverlayAllotmentWithBasin" + str(dt.now().microsecond))
Input
basin_overlay
Output
OverlayAllotmentWithBasin172198
Feature Layer Collection by arcgis_python
Last Modified: June 24, 2019
0 comments, 2 views
Input
map2 = gis.map('Oregon')
map2
Output
Input
map2.add_layer(basin_overlay.layers[0])

The new features have all the attributes of the features in the input layers. In this case, the new allotment features are assigned the attributes-including the name and ID-of the hydrologic basin they fall within. Allotments that fall in two or more basins are split at the basin boundaries and the corresponding attributes assigned to each portion of the allotment.

Mapping and exploring basin overlay results

We will explore overlay results using both matplotlib and plot method of Spatially Enabled Dataframe.

Let's read the overlay layer as Spatially Enabled DataFrame.

Input
sdf = pd.DataFrame.spatial.from_layer(basin_overlay.layers[0])
Input
slctd_cols = ['BASIN_NAME', 'allot_name', 'allot_no', 'REGION', 'HUC', 'SHAPE']
Input
basin_overlay_df = sdf[slctd_cols]
Input
basin_overlay_df.head()
Output
BASIN_NAME allot_name allot_no REGION HUC SHAPE
0 COQUILLE WICKENS/BAXTER 20001 17 17100305 {"rings": [[[-13812226.4688, 5342337.4032], [-...
1 UMPQUA KELLOGG 10007 17 17100303 {"rings": [[[-13752329.3678, 5398012.2305], [-...
2 COQUILLE PRINCEOTTO/CAWTHON 20001 17 17100305 {"rings": [[[-13803327.8105, 5347044.5693], [-...
3 UMPQUA BULLOCK 20006 17 17100303 {"rings": [[[-13752015.6692, 5396937.5365], [-...
4 LOST JOHNSON 00833 18 18010204 {"rings": [[[-13504148.0167, 5160990.0443], [-...

We will group the dataframe by 'BASIN_NAME'. To get the number of basins in each group we will use the size() method.

Input
grp_basin = basin_overlay_df.groupby('BASIN_NAME').size()
Input
grp_basin.head()
Output
BASIN_NAME
 ALVORD LAKE           51
 APPLEGATE              2
 BEAVER-SOUTH FORK     52
 BROWNLEE RESERVOIR    52
 BULLY                 29
dtype: int64
Input
grp_basin.nlargest(10)
Output
BASIN_NAME
 POWDER                  156
 BURNT                   146
 LOWER JOHN DAY          128
 UPPER MALHEUR           122
 LOST                     84
 LOWER CROOKED            72
 HARNEY-MALHEUR LAKES     71
 SUMMER LAKE              69
 UPPER CROOKED            65
 WILLOW                   62
dtype: int64
Input
grp_basin.nlargest(10).plot(kind='barh')
Output
<matplotlib.axes._subplots.AxesSubplot at 0x1984f3484a8>

We see that Powder Basin intersected the largest number of grazing allotments, followed by Burnt and Lower John Day basins.

We will now map grazing allotments by hydrologic basin.

Input
map3 = gis.map('oregon')
map3
Output
Grazing allotments are color coded by the basin they are in. Allotment features are split where a basin boundary crosses them. Selecting an allotment feature displays the allotment and basin information--the basin information is associated with the feature.
Input
basin_overlay_df.spatial.plot(kind='map',
                              map_widget=map3,
                              renderer_type='u',
                              col='BASIN_NAME') 
Output
True

Display grazing allotments within a particular basin

Let's get a list of all basin names.

Input
basin_overlay_df.BASIN_NAME.str.strip().unique()
Output
array(['COQUILLE', 'UMPQUA', 'LOST', 'GOOSE LAKE', ..., 'MIDDLE SNAKE-SUCCOR', 'THOUSAND-VIRGIN', 'SOUTH FORK OWYHEE',
       'KLICKITAT'], dtype=object)

We will apply a filter to visualize grazing allotments within a particular basin.

Input
john_day_df = basin_overlay_df[basin_overlay_df['BASIN_NAME'] == ' MIDDLE FORK JOHN DAY']
Input
john_day_df
Output
BASIN_NAME allot_name allot_no REGION HUC SHAPE
58 MIDDLE FORK JOHN DAY MIDDLE FORK 04014 17 17070203 {"rings": [[[-13234233.0817, 5586868.7282], [-...
238 MIDDLE FORK JOHN DAY SIDEHILL 04026 17 17070203 {"rings": [[[-13238795.622, 5579105.8225], [-1...
280 MIDDLE FORK JOHN DAY GIBSON CREEK 04135 17 17070203 {"rings": [[[-13263280.5074, 5602968.1923], [-...
552 MIDDLE FORK JOHN DAY PASS CREEK 04184 17 17070203 {"rings": [[[-13268669.5057, 5592601.4069], [-...
754 MIDDLE FORK JOHN DAY NORTH FORK 04029 17 17070203 {"rings": [[[-13278826.6115, 5607309.3595], [-...
802 MIDDLE FORK JOHN DAY JINKS CREEK 04050 17 17070203 {"rings": [[[-13282326.8528, 5597674.5336], [-...
862 MIDDLE FORK JOHN DAY WEST FORK BURNT RIVER 15324 17 17070203 {"rings": [[[-13170488.1642, 5557278.5571], [-...
871 MIDDLE FORK JOHN DAY MUD SPRINGS 04015 17 17070203 {"rings": [[[-13269728.7602, 5614038.7], [-132...
985 MIDDLE FORK JOHN DAY NORTH FORK BURNT RIVER (USFS) 15329 17 17070203 {"rings": [[[-13179809.2132, 5565524.2035], [-...
989 MIDDLE FORK JOHN DAY DOHERTY P JOE SHEEP 04193 17 17070203 {"rings": [[[-13251986.3822, 5614102.6014], [-...
1152 MIDDLE FORK JOHN DAY SLICKEAR MTN. 04003 17 17070203 {"rings": [[[-13269794.5899, 5593722.8408], [-...

Let's plot the filtered results on map.

Input
map4 = gis.map('oregon', zoomlevel=8)
map4
Output
Grazing allotments are filtered by a particular hydrologic basin (Middle Fork John Day). The map shows all the allotments (or portions of allotments) that are in the basin, and the table above lists them with the associated information.
Input
map4.center = [44.88, -118.83]
Input
map4.add_layer(hydro_units)
Input
john_day_df.spatial.plot(map_widget=map4)
Output
True

Display a particular grazing allotment to see which basins intersect it

We will filter grazing allotments using allotment number 04003.

Input
allot_df = basin_overlay_df[basin_overlay_df['allot_no'] == '04003']
Input
allot_df
Output
BASIN_NAME allot_name allot_no REGION HUC SHAPE
1152 MIDDLE FORK JOHN DAY SLICKEAR MTN. 04003 17 17070203 {"rings": [[[-13269794.5899, 5593722.8408], [-...
1153 NORTH FORK JOHN DAY SLICKEAR MTN. 04003 17 17070202 {"rings": [[[-13276660.1969, 5581616.3978], [-...
Input
map5 = gis.map('oregon', zoomlevel=8)
map5
Output

Grazing allotments are filtered using a particular allotment number (04003). The map shows in which basins the allotment is.

Input
map5.center = [44.88, -118.83]
Input
map5.add_layer(hydro_units)
Input
allot_df.spatial.plot(map_widget=map5)
Output
True

Assigning allotment information to streams

We will again use overlay_layers tool, this time overlaying streams with grazing allotments (area features can be overlaid with line or point features as well as other area features).

We will overlay streams with grazing allotments using the Intersect option. The output layer contains only those stream segments that cross a grazing allotment.

Input
stream_overlay = overlay_layers(streams,
                                grazing_allotments,
                                overlay_type='Intersect',
                                output_name="StreamOverlay" + str(dt.now().microsecond)
                               )
Input
stream_overlay
Output
StreamOverlay58569
Feature Layer Collection by arcgis_python
Last Modified: June 24, 2019
0 comments, 2 views

Read the overlay layer as Spatially Enabled DataFrame.

Mapping and exploring stream overlay results

Input
stdf = pd.DataFrame.spatial.from_layer(stream_overlay.layers[0])
Input
cols = ['allot_name', 'allot_no', 'HUC', 'PNAME', 'SHAPE' ]
Input
stream_overlay_df = stdf[cols]
Input
stream_overlay_df.head()
Output
allot_name allot_no HUC PNAME SHAPE
0 BIG SUMMIT WEST 02580 17070304 JOHNSON CR {"paths": [[[-13383245.2218, 5520910.6927], [-...
1 WEIGAND 00075 17070304 CROOKED R {"paths": [[[-13401832.2407, 5485427.1285], [-...
2 WEST PINE CREEK 00076 17070304 CROOKED R {"paths": [[[-13401868.5519, 5485480.4708], [-...
3 NORTH FORK 04029 17070202 MALLORY CR {"paths": [[[-13279226.9627, 5616433.4654], [-...
4 BIG MUDDY 02512 17070204 *G {"paths": [[[-13419994.254, 5593063.7317], [-1...
Input
st_grp = stream_overlay_df.groupby('PNAME').size()
Input
st_grp.nlargest(10).plot(kind='barh')
Output
<matplotlib.axes._subplots.AxesSubplot at 0x1984be31f98>

John Day River is the third longest free flowing river in contiguous US. The plot shows that this river is mostly used for ranching.

Input
stream_overlay_df.groupby('allot_name').size().nlargest(10)
Output
allot_name
THREE FINGERS           36
JACKIES BUTTE SUMMER    23
QUARTZ MOUNTAIN         23
SADDLE BUTTE            23
SUMMIT PRAIRIE          22
WALLROCK                21
G.I.                    16
BIG BUTTE               15
BOARD CORRALS           15
COYOTE LAKE             14
dtype: int64

Thirty-six (36) stream segments from 17 different streams pass through the 'THREE FINGERS' grazing allotment.

Input
map6 = gis.map('oregon')
map6
Output
Map of streams overlaid by grazing allotments.
Input
map6.add_layer(grazing_allotments)
Input
stream_overlay_df.spatial.plot(kind='map',
                               map_widget=map6,
                               renderer_type='u',
                               col='PNAME') 
Output
True

Display a stream and grazing allotments it intersects

We will filter the analysis_layer for a specific stream—BRIDGE CR.

Input
creek_df = stream_overlay_df[stream_overlay_df['PNAME'] == 'BRIDGE CR']
Input
creek_df.head()
Output
allot_name allot_no HUC PNAME SHAPE
155 BURNT RANCH 02624 17070204 BRIDGE CR {"paths": [[[-13392864.8203, 5579859.7598], [-...
156 SUTTON MTN. 02533 17070204 BRIDGE CR {"paths": [[[-13392552.7024, 5577996.559], [-1...
771 UPPER BRIDGE CREEK 00701 17120005 BRIDGE CR {"paths": [[[-13488884.9907, 5323673.5952], [-...
772 BUCK CREEK-BRIDGE CR 00702 17120005 BRIDGE CR {"paths": [[[-13485835.5035, 5327100.6087], [-...
875 RING BUTTE 10208 17050116 BRIDGE CR {"paths": [[[-13160937.9964, 5491764.6228], [-...
Input
map7 = gis.map('oregon')
map7
Output

This map shows in which allotment each stream segment falls.

Input
map7.add_layer(grazing_allotments)
Input
creek_df.spatial.plot(map_widget=map7)

Conclusion

Biologists can now identify which hydrologic basin and stream(s) intersect with which grazing allotments in an effort to identify sources of chronic water quality issues.

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