NLDI.get_basins() returns no features for almost all gages

[wkeenan21]

What happened?

I have a list of 216 USGS gage IDs. Most of them are in Colorado. This code errors out on all of them due to "returns no features".

What did you expect to happen?

I expected it to grab the delineated watershed for each gage.

Minimal Complete Verifiable Example

gage_ids = ['07203000', '07216500', '07218000', '08247500', '08248000', '08255500', '08265000', '08266820', '08267500', '08275500', '08276300', '08277470', '08279000', '08284100', '08289000', '08291000', '08294210', '08315480', '08316000', '08377900', '08378500', '08380500', '08386505', '08400000', '08401200', '08401900', '08405150', '08405450', '08405500', '08406000', '08477110', '09386900', '09430600', '09266500', '09277500', '09279000', '09288000', '09291000', '09292000', '09292500', '09296800', '09310500', '09312600', '09337000', '09378170', '09378630', '09403600', '09404900', '09405500', '09408000', '09408400', '09409100', '09413900', '10011500', '10105900', '10109000', '10113500', '10128500', '10130500', '10131000', '10134500', '10137500', '10140100', '10146400', '10149000', '10154200', '10157500', '10164500', '10172700', '10172800', '10172952', '10173450', '10194200', '10205030', '10215900', '10234500', '06614800', '06710385', '06712000', '06715000', '06716100', '06718550', '06719505', '06730200', '06741510', '06746110', '06751150', '06751490', '07079300', '07083000', '07096250', '07103700', '07103703', '07103797', '07103800', '07105000', '07108900', '07124200', '08245000', '08246500', '09010500', '09024000', '09025000', '09025300', '09026500', '09032000', '09032100', '09034900', '09035700', '09035900', '09036000', '09037500', '09038500', '09041400', '09046490', '09047700', '09050100', '09050700', '09051050', '09059500', '09063900', '09064600', '09065100', '09065500', '09066000', '09066200', '09066300', '09066325', '09066510', '09067020', '09067200', '09073300', '09073400', '09080400', '09081600', '09109000', '09110000', '09113980', '09118450', '09124500', '09126000', '09132500', '09143500', '09146200', '09147025', '09147500', '09165000', '09166950', '09237450', '09239500', '09246200', '09253000', '09255000', '09306242', '09342500', '09353800', '09354500', '09358000', '09358550', '09359010', '09359020', '09371000', '09371520', '09372000', '401723105400000', '401733105392404', '06309200', '06622700', '06623800', '06625000', '06632400', '09188500', '09196500', '09210500', '09223000', '09229500', '10032000', '13010065', '13011900', '13016305', '13016450', '13018300', '06814000', '06853800', '06863500', '06871500', '06876700', '06878000', '06879650', '06885500', '06888500', '06889200', '06889500', '06890100', '06891500', '06892000', '06910800', '06911490', '06911900', '06912500', '06917000', '07144910', '07145700', '07147070', '07167500', '07169800', '07170700', '07172000', '07179500', '07179795', '09119000', '09304500', '06620000', '08227000', '09057500', '09105000']


basins_gdf = gpd.GeoDataFrame(columns=['gage', 'geometry'], crs="EPSG:4326")

# Loop through the gage IDs and get the basin for each
for gage_id in gage_ids:
    try:
        print(gage_id)
        # Retrieve the watershed basin using PyNHD
        basin = NLDI().get_basins(f"USGS-{gage_id}", split_catchment=True)

        # Append to the main GeoDataFrame
        basins_gdf = gpd.GeoDataFrame(pd.concat([basins_gdf, basin], ignore_index=True), crs=basins_gdf.crs)
        print(f'got: {gage_id}')
    except Exception as e:
        print(f"Failed to get basin for gage {gage_id}: {e}")

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

Anything else we need to know?

No response

Environment

Details

SYS INFO

commit: None
python: 3.12.8 | packaged by Anaconda, Inc. | (main, Dec 11 2024, 16:48:34) [MSC v.1929 64 bit (AMD64)]
python-bits: 64
OS: Windows
OS-release: 11
machine: AMD64
processor: Intel64 Family 6 Model 151 Stepping 2, GenuineIntel
byteorder: little
LC_ALL: None
LANG: None
LOCALE: ('English_United States', '1252')
PACKAGE VERSION

async-retriever 0.19.3
pygeoogc 0.19.3
pygeoutils 0.19.5
py3dep N/A
pynhd 0.19.3
pygridmet N/A
pydaymet N/A
hydrosignatures 0.19.0
pynldas2 N/A
pygeohydro 0.19.0
tiny-retriever N/A
aiodns 3.2.0
aiofiles 24.1.0
aiohttp 3.11.14
aiohttp-client-cache 0.12.4
aiosqlite 0.21.0
brotli 1.1.0
cytoolz 1.0.1
orjson 3.10.16
numpy 2.2.4
pandas 2.2.3
scipy 1.15.2
xarray 2025.3.0
numba N/A
numbagg N/A
click 8.1.8
geopandas 1.0.1
rasterio 1.4.3
rioxarray 0.18.2
shapely 2.0.7
netcdf4 1.7.2
pyproj 3.7.1
defusedxml 0.7.1
folium 0.19.4
h5netcdf 1.5.0
matplotlib 3.10.0
planetary-computer N/A
pystac-client N/A
joblib 1.4.2
multidict 6.2.0
owslib 0.33.0
requests 2.32.3
requests-cache 1.2.1
typing-extensions 4.13.0
url-normalize 1.4.3
urllib3 2.3.0
yarl 1.18.3
networkx 3.4.2
pyarrow 19.0.1
py7zr N/A
flox N/A
opt-einsum N/A

[BaptisteFrancois]

Hello,
I have the same problem. Out of 3,246 basins spread across the USGS, more than half returns with no features. I have used the exact same script a few months ago with the exact same basins, and I could get features for all of them.

I don't post a script snippet, because it would be the exact same as the one posted already.

[cheginit]

@wkeenan21 the issue is with split_catchment=True. The NLDI web service has been having some issues with this option recently. So, unfortunately, until it's fixed it cannot be used.

@BaptisteFrancois and @wkeenan21 note that NLDI recently has set a rate limit on their web service (3600 requests per IP per hour), so if your request is large you will get an error.

[cheginit]

The NLDI web service, under-the-hood, uses nldi-flowtools for splitting catchment. So you can use this package directly if this option is important for your use case.

[BaptisteFrancois]

Thanks for your answer. I did not know about the 3,600 requests per IP per hour.

As an update, when I tried again to run my script over the weekend, I could get all the 3,246 basins but 42.
When I tried again this morning, I could get all of them but 3.

I assume the issue was likely coming from the NLDI web service and someone must have done something to fix whatever the issue was.

[wkeenan21]

Thanks for your answer. I did not know about the 3,600 requests per IP per hour.

As an update, when I tried again to run my script over the weekend, I could get all the 3,246 basins but 42. When I tried again this morning, I could get all of them but 3.

I assume the issue was likely coming from the NLDI web service and someone must have done something to fix what whatever the issue was.

@BaptisteFrancois Could you share your code? My script is still failing for "Service returned no features." on 95% of my gages.

[BaptisteFrancois]

It is literally the same as yours, except that I do not use the option spilt_catchment=True.
Cheginit mentions the option is broken currently and can't be used.

Have you tried we this option set to False?

[wkeenan21]

Oh, I thought you were getting good responses with split_catchment=True.

It does work with the option set to False, but this option is important for my workflow.

Thanks

[cheginit]

@wkeenan21 have you recently tried split_catchment, it seems that they solved the issue on the server side.