Finding Watershed Triple Points

Many years ago a friend shared the concept of watershed double points. These occur when adjacent watersheds have their highest point in close proximity to each other. This often occurs on mountain peaks - water that falls on one side goes to one river, the opposite side to another river. For high peaks, this area is often the highest point in both watersheds.

Triple points - where three watersheds come together at their highest point - are much rarer. My friend seemed to think there were a few in Vermont, but he couldn’t name the locations or the watersheds. Today while my kid was napping I decided to do a little analysis in postgis.

Data

I grabbed the NHD dataset from the National Map at this link, then used ogr2ogr to import the HUC10 layer into postgis, reprojecting to VT State Plane while doing so:

ogr2ogr -f Postgresql 'PG:user=charliehofmann' NHD_H_Vermont_State_GPKG.gpkg \
  WBDHU10 -nln huc10_vt -lco GEOMETRY_NAME=geom -overwrite -nlt PROMOTE_TO_MULTI \
  -t_srs EPSG:32145

HUC10 seemed to be the right level of specificity for this - plenty of watersheds with familiar names - Black River, Williams River, three branches of the White, etc.

I then downloaded a 10m DEM from VCGI from this link and brought it into postgis with:

psql -c "create extension postgis_raster;"

raster2pgsql -Y -d -t 256x256 -N '-32768' -I -C -M -n "path" hdr.adf vt_dem_10 | psql

Hat tip to Matt Perry for the canonical post on postgis raster analysis - I reference it whenever I need to deal with raster data.

Analysis

First I needed to find the max elevation for each watershed - running zonal stats for my HUC10 geometries against the DEM.

I’m not super familiar with the postgis raster functions, but again Matt’s blog post was super helpful to get me started:

create table max_elevation_by_watershed as
select
  objectid,
  huc10_vt.name AS name,
  (ST_SummaryStatsAgg(ST_Clip(raster.rast, huc10_vt.geom, true), 1, true)).max as max_elevation
from
  vt_dem_10 as raster
inner join huc10_vt on
  ST_Intersects(huc10_vt.geom, raster.rast)
group by
  name, objectid;

I then had this nice table:

Next I needed to find where these max elevation pixels actually were in each watershed. I decided to use ST_Reclass to reclassify the DEM to only 1 for the max elevation pixel and NODATA for everything else.

After finding each pixels that matched the max for each watershed, I converted them to points, then found the closest point to the watershed geometry in the event that any were just outside. Man, postgis is cool.

create table highest_point_by_watershed as
-- find a single point for each watershed, ordering by the distance to
-- the watershed geom.
-- as we're talking about tiles here, the reclass will likely have some
-- points outside the watershed geometry itself, but at least one that intersects
select distinct on (objectid) objectid, name, pt_geom as geom
from (
  select
    objectid,
    huc10_vt.name,
    huc10_vt.geom,
    (
      -- dump the pixels with data (in this case the max pixels) to points
      ST_PixelAsPoints(ST_Reclass(
        rast,
        1,
        -- dynamically reclassify the raster based on the max value for that watershed
        '0-' || (max_elevation - 1) || ':-32768, ' || max_elevation || ':1,' || (max_elevation + 1) || '-10000:-32768',
        '16BSI',
        -32768), 1, true
      )
    ).geom pt_geom
  from vt_dem_10
  -- filter the huc10_vt raster to the tiles that intersect each watershed
  join huc10_vt on ST_Intersects(geom, rast)
  join max_elevation_by_watershed using (objectid)
) reclassed
order by objectid, ST_Distance(pt_geom, geom);

Finally, cluster the watershed points to find groups of three within 100 ft.

select json_build_object(
  'type', 'FeatureCollection',
  'features', json_agg(ST_AsGeoJSON(clustered.*)::json)
)
from (
  select name, ST_Transform(geom, 4326),
    ST_ClusterDBSCAN(geom, eps := 100, minpoints := 3) over () AS cid
  from highest_point_by_watershed
) clustered
where cid is not null; -- filter out non-clustered points

Results

We did it! There are three triple point locations within the state of Vermont!

There’s the Winooski / Browns / Lamoille River triple point just north of the nose on Mount Mansfield:

The Barton / Clyde / Passumpsic on Bald Mountain:

And maybe my personal favorite, the Waits / Wells / Winooski:

Thanks for playing everyone! Please direct any further Vermont geographical puzzles my way . . . between VCGI’s data and PostGIS we should be able to answer a lot of them!