How to use Python to delineate watersheds in the continental US with USGS web services

I recently tried my hand at using NLDI, the new water data web service from the U.S. Geological Survey. The service provides access to water data from the National Hydrography Dataset (NHD), version 2. I was specifically interested in using the NHD geodata to delineate watersheds. This data was digitized from 1:100,000-scale maps. (There is newer, more detailed data, called NHDPlus HR, for high-resolution, but no web service for these data yet.)

Snake River watershed upstream of the Lower Granite Dam in eastern Washington state (just a location I picked at random for this demo!)

It took me a couple days of tinkering to get the results I wanted, but the results look impressive and seem accurate. You can run the code yourself by downloading one of these files from GitHub Gists. (Or follow along below.)