nhd_query, use polygon/bbox to select

[glorious1]

nhd_plus_query allows the use of a polygon, but nhd_query does not. Would be nice if it did.

The point + buffer gives me 0 features. It's unclear if I can just use lon/lat for the point or it requires an sf point object broken somehow into lon and lat.

Edit: it looks like my lon/lat point is being converted to UTM in meters, which is why my buffer_dist in degrees is not working. This is way out of line with the documentation.

[jsta]

Happy to help with this. Could you provide some sample code for testing? Do you mean the documentation for the nhd_query function itself (`?nhd_query) or in the README/vignettes?

[glorious1]

Thank you. The documentation I'm referring to is ?nhd_query. As I read it, I provide lon and lat in degrees and buffer_dist, also in degrees. The fact that the default buffer_dist is 0.05 suggests it is degrees.

However, as shown below, I can raise buffer_dist to 2, which would be about 1/4 of Colorado and should contain many thousands of flowlines. Yet I get 0. I can get some if I raise it to many thousand.

The output suggests that, for some reason, the point is being converted to UTM, so I assume that's why buffer_dist is interpreted as meters, not degrees. The flowline object is also in UTM. I don't understand this behavior.

> fl <- nhd_query( lon=-107.2, lat=39.45, dsn="NHDFlowline", buffer_dist=2 )
Reading layer `NHDFlowline' from data source 
  `/Users/jim/Documents/R/nhdR_data/NHD_H_Colorado_State_GDB.gdb' using driver `OpenFileGDB'
Simple feature collection with 1278092 features and 15 fields
Geometry type: MULTILINESTRING
Dimension:     XYZM
Bounding box:  xmin: -110.4218 ymin: 35.71496 xmax: -100.64 ymax: 41.97636
z_range:       zmin: 0 zmax: 0
m_range:       mmin: 0 mmax: 100
Geodetic CRS:  NAD83
> fl
$pnt
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 310701.7 ymin: 4369025 xmax: 310701.7 ymax: 4369025
CRS:           +proj=utm +zone=13 +datum=WGS84
POINT (310701.7 4369025)

$sp
$sp$NHDFlowline
Simple feature collection with 0 features and 15 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
CRS:           +proj=utm +zone=13 +datum=WGS84
 [1] Permanent_Identifier        FDate                       Resolution                 
 [4] GNIS_ID                     GNIS_Name                   LengthKM                   
 [7] ReachCode                   FlowDir                     WBArea_Permanent_Identifier
[10] FType                       FCode                       MainPath                   
[13] InNetwork                   VisibilityFilter            Shape_Length               
[16] Shape                      
<0 rows> (or 0-length row.names)

[jsta]

I'd recommend specifying the buffer_dist as a units object instead of a bare numeric. This was in the nhd_plus_query docs but not nhd_query. I'll push a commit syncing the two shortly. At 40 deg lat, 2 deg is approximately 190 km, so the call would be:

fl <- nhd_query(lon = -107.2, lat = 39.45,
  dsn = "NHDFlowline", buffer_dist = units::set_units(190, "km"))

[glorious1]

Thank you. Speaking of nhd_plus_query, my other question. It allows selecting by a polygon, which is what I would prefer, but nhd_query does not. Could it, or is there a reason why not?

To be honest, I sought out nhdR because I was unable to get st_read to work on NHD gdb's with wkt_filter, which allows selecting with a polygon. This does work with other spatial gdb layers, but for some reason not with NHD. Are you aware of why that might be?

[jsta]

Just not implemented yet. I wasn't using NHD (regular) when I developed the wkt_filter functionality. Giving it a try on nhd_query now...

[jsta]

Should work now, with something like:

library(nhdR)
bbox <- sf::st_bbox(
  c(xmin = -107.45, ymin = 39.2, xmax = -106.95, ymax = 39.7),
  crs = 4326)
qry_lines <- nhd_query(poly = sf::st_as_sfc(bbox), dsn = "NHDFlowline")
nrow(qry_lines$sp$NHDFlowline) # == 17199

[glorious1]

Awesome, thank you.