Add tributary catchments to NHD-PRMS crosswalk

[lekoenig]

In reviewing PR #38 I wanted to check that the summed areas for the contributing NHD catchments were approximately equal to the area of the PRMS catchment polygons associated with each reach. This check shows that for many catchments, the PRMS catchment area is much greater than our summed NHD area calculated from the comid's in the crosswalk target:

This is because our crosswalk was focused on finding those NHD reaches that intersect the PRMS segments, which excludes the headwater tributaries (and their catchments) that are not included in the PRMS network:

Therefore, in order to collate land use information as in PR #38 , we will need to add a column to the crosswalk target that includes these headwater catchments that contribute to the PRMS-scale catchment.

[lekoenig]

I've been working on a fix to gather all of the NHD catchments that would directly drain to a given PRMS segment (as described above, this isn't correct in the current xwalk table for larger PRMS catchments that are missing contributing tributary COMID's).

Here's my current solution:

1. Each PRMS segment has a unique value of comid_down, the NHD COMID at the bottom of that segment. The GF also refers to these NHD reaches as points of interest or POIs. For each comid_down, find all upstream COMIDs.
2. Find other POIs in the xwalk table that are located upstream of the target PRMS reach. For each of those, find all upstream COMIDs
3. By difference, identify the COMIDs that are returned in the subset of reaches from 1) but not located in 2). These are the COMIDs that intersect the local PRMS catchment.

I've attempted to display those steps for the example catchment above here:

In general, this really improves our PRMS-area to summed-NHD-area relationship:

However, the PRMS polygons are not always approximately identical to the dissolved polygon from NHD. These two areas may be slightly different in area and/or in land area covered. I suspect, but am not sure, that this is because the GF v1 was created from NHDPlusV1 and we're working with NHDPlusV2. These differences would potentially have implications for attributes that we estimate from raster processing (i.e. mask using PRMS HRU polygons) versus from attributes referenced to NHDPlusV2 such as the MasterList dataset, so I appreciate any comments/questions/or feedback you have here.

[lekoenig]

Here I wanted to dig a little deeper into the patterns described for the fix above. For ~90% of the PRMS segments in DRB that have corresponding HRU's, the two areas (PRMS HRU vs summed NHD areas) are within 10 km2:

I estimate 28 PRMS segments with discrepancies between the two areas > 10 km2 (red points in scatterplot, some of which are labeled by subsegid):

For 3 of those segments, the discrepancy comes from the fact that the PRMS network has been split as part of the delaware-model-prep pipeline but the PRMS HRU's are still the same (e.g. for p1_reaches_sf$subsegseg = "3", which is associated with both subsegid == 3_1 and subsegid == 3_2; this also applies to subsegid == 8_2):

Then there are cases (I count 6) where the corresponding HRU appears to actually include catchments for multiple PRMS segments, not just the "matched" segment. For example, the HRU's labeled as belonging to subsegseg == 157 seem to be comprised of catchments for three different PRMS segments (subsegid's 155_1, 156_1, and 157_1):

Here's another example of that - the HRU's supposedly belonging to subsegseg == 379_1 seem to actually belong to three different PRMS segments (subsegid's 379_1, 378_1, and 382_1). This PRMS segment/HRU was mentioned in a previous issue within delaware-model-prep.

For the remaining 19 segments (from those 28 segments with contributing area discrepancies > 10 km2), there appear to just be boundary discrepancies between the contributing area represented by the HRU polygons and the assigned NHD catchments. For example, for subsegid == 354_1, part of the land area represented by the corresponding HRU gets assigned to a different NHDPlusV2 reach (and thus a different PRMS segment). These discrepancies range from 10 - 67 km2.

The NHD catchments here are derived from the workflow described in the comment above. This function uses helper functions in nhdplusTools to navigate the upstream network. The advantage is it's fast and secondarily, it alerts us to oddities with the HRU polygons, including that ~10% of the PRMS reaches don't have a corresponding HRU, and that some HRU's seem misspecified (see the first two examples above). The disadvantage of this approach is that the specific land area covered by the resulting NHD sometimes differs from the PRMS HRU's.

Alternatively, we could assume that the PRMS HRU's represent the appropriate land area and just find the NHDPlusV2 catchments that intersect the PRMS polygons (e.g. at least 90% of the NHD catchment area overlaps with the PRMS HRU's). The advantage of this approach is that the areas are well-aligned (see below) and therefore the land areas we're using to derive any attributes are pretty much consistent.

The disadvantages are that this adds a spatial processing step to the cross walk target/is slow, we don't address the HRU edge cases identified above, and for the ~40 segments without HRU's, we'd either have to derive the NHD catchments using another method or just forego gathering catchment attributes for those segments. I checked on this and it looks like there are 39 SC sites that overlap PRMS segments without HRU's, including 2 NWIS sites (1,765 observation days, which represents ~1% of our NWIS days).

[jds485]

Thanks for your investigative work here! This is becoming more challenging a problem than expected, but will ultimately be useful for modeling.

For 3 of those segments, the discrepancy comes from the fact that the PRMS network has been split as part of the delaware-model-prep pipeline but the PRMS HRU's are still the same

For these, it sounds like we should update the PRMS HRU area to the NHD area because that's more representative of the actual area for that sub-segment. Would that change apply to all reaches with _n and n > 1?

there are cases (I count 6) where the corresponding HRU appears to actually include catchments for multiple PRMS segments

Similarly, we can update the HRU area with only the area for the corresponding PRMS segment.

For the remaining 19 segments...there appear to just be boundary discrepancies between the contributing area represented by the HRU polygons and the assigned NHD catchments.

It's surprising to see such large differences. I have a compromise suggestion merging your 2 approaches to get a faster processing time and (hopefully) good accuracy. Could we take only these 19 segments, get their HRUs and neighboring HRUs, then use the NHD intersection method for only those HRUs?

For the ~40 segments without HRU's, we'd either have to derive the NHD catchments using another method or just forego gathering catchment attributes for those segments.

Wouldn't we have this problem with the 1st method, too? Or does the nhdplusTools compute catchments for these segments?

[lekoenig]

Thanks, Jared. For the majority of segments it won't really matter whether we use the spatial join or gather the upstream COMID's, so I've taken the compromise approach you describe above. Therefore, for most segments we gather the upstream COMID's but for the segments with larger area discrepancies between the corresponding HRU's and the summed NHD areas that were identified by visual inspection, we'll derive the contributing COMID's by a spatial join. To answer your questions above:

For these, it sounds like we should update the PRMS HRU area to the NHD area because that's more representative of the actual area for that sub-segment. Would that change apply to all reaches with _n and n > 1?

Yeah, I do think that for all reaches with _n and n>1, the NHD area is more representative of the actual area that drains to the segment in question. There are not many, but I think it's subsegseg = 3, 8, and 12.

Wouldn't we have this problem with the 1st method, too? Or does the nhdplusTools compute catchments for these segments?

I'm not really sure of the reason why a PRMS segment would lack an HRU (aside from a handful of segments where I think their catchment just got lumped in with another HRU (e.g. 155_1, 377_1). Generally speaking those segments still have NHDPlusV2 catchments (that can be retrieved using nhdplusTools), so we still get back catchments using the first method even if there is no HRU polygon associated with that segment as part of the Geospatial Fabric. An exception is 287_1, which Margaux has pointed out before - it has one NHDPlusV2 catchment and the area for that catchment is given as zero within the value added attributes table.

[jds485]

Thanks! Just to understand the final plan:

There are not many, but I think it's subsegseg = 3, 8, and 12

Are you going to use the NHD catchment areas for these, then? I guess that would correct areas for all _n reaches for those subsegsegs.

so we still get back catchments using the first method even if there is no HRU polygon associated with that segment as part of the Geospatial Fabric

Are you also using the NHD areas for these ~40? Except 287

[lekoenig]

Are you going to use the NHD catchment areas for these, then? I guess that would correct areas for all _n reaches for those subsegsegs.

Yep, in the PR that is currently open the xwalk target will return those COMIDs that fall within the subset areas rather than the full HRU. So that means we'd be using the NHD catchment areas for processing the VarsOfInterest, for example.

However, we don't currently replace any catchments within the spatial data, so Margaux's raster processing presumably still uses the HRU polygons. I think we could download and dissolve the NHD catchment polygons for these few segments if we decide that's the way to go.

Are you also using the NHD areas for these ~40? Except 287

Yeah, the xwalk target will return the upstream COMIDs that drain directly to those segments. So even though we don't have corresponding HRU polygons for those 40, we still have COMID's and therefore areas that we'll use for processing the VarsOfInterest. I'll also note that we still return a COMID in the crosswalk table for 287, but because the lone returned COMID has AREASQKM = 0, it will functionally (or explicitly) get omitted when we process the area-weighted-mean VarsOfInterest.

[jds485]

However, we don't currently replace any catchments within the spatial data, so Margaux's raster processing presumably still uses the HRU polygons. I think we could download and dissolve the NHD catchment polygons for these few segments if we decide that's the way to go.

I think it would make sense to update the polygons for these segments. Do you think that should be a new issue or be part of one of the open PRs/issues?

[lekoenig]

I think it would make sense to update the polygons for these segments. Do you think that should be a new issue or be part of one of the open PRs/issues?

For modularity, I would propose this be a new issue that potentially gets addressed along with issue #60. I think the open NLCD PR just uses the COMIDs (and not the catchment polygons), but the open FORE-SCE PR does use the catchment polygons.

[jds485]

Good suggestion to address with #60. I'll link the comment to that issue