This page was generated from webservices.ipynb. Interactive online version:
Web Services Queries#
import pygeoutils as geoutils from pygeoogc import WFS, WMS, ArcGISRESTful, ServiceURL from pynhd import NLDI
import warnings warnings.filterwarnings("ignore", message=".*Content metadata for layer.*")
PyGeoOGC and PyGeoUtils can be used to access any web services that are based on ArcGIS RESTful, WMS, or WFS. It is noted that although, all these web service have limits on the number of objects (e.g., 1000 objectIDs for RESTful) or pixels (e.g., 8 million pixels) per requests, PyGeoOGC takes care of dividing the requests into smaller chunks under-the-hood and then merges them.
Let’s get started by retrieving a watershed geometry using NLDI and use it for subsetting the data.
basin = NLDI().get_basins("11092450") basin_geom = basin.geometry
PyGeoOGC has a
ServiceURL that contains URLs of the some of the popular web services. Let’s use it to access NHDPlus HR Dataset RESTful service and get all the catchments that our basin contain and use
pygeoutils.json2geodf to convert it into a GeoDataFrame.
hr = ArcGISRESTful(ServiceURL().restful.nhdplushr, 10, outformat="json") oids = hr.oids_bygeom(basin_geom, "epsg:4326", spatial_relation="esriSpatialRelContains") resp = hr.get_features(oids) catchments = geoutils.json2geodf(resp)
oids_bygeom has an additional argument for passing any valid SQL WHERE clause to further filter the data on the server side. For example, let’s only keep the ones with an area of larger than 0.5 sqkm.
oids = hr.oids_bygeom(basin_geom, geo_crs=4326, sql_clause="AREASQKM > 0.5") catchments = geoutils.json2geodf(hr.get_features(oids))
ax = catchments.plot(figsize=(8, 8)) ax.figure.savefig("_static/sql_clause.png", bbox_inches="tight", facecolor="w")
We can also submit a query based on IDs of any valid field in the database. If the measure property is desired you can pass
True to the
get_features class method:
oids = hr.oids_byfield("NHDPlusID", [5000500013223, 5000400039708, 5000500004825]) resp = hr.get_features(oids, return_m=True) catchments = geoutils.json2geodf(resp)
Next, let’s get DEM using the 3D Elevation Program WMS service. First we need to connect to the service using
wms = WMS( ServiceURL().wms.nm_3dep, layers="3DEPElevation:None", outformat="image/tiff", crs=5070, )
Then we can get the data using the
getmap_bybox function. Note that this function only accepts a bounding box, so we need to pass a bounding box and mask the returned data later on using
bbox = basin_geom.bounds r_dict = wms.getmap_bybox( bbox, 100, box_crs=4326, ) dem = geoutils.gtiff2xarray(r_dict, bbox, 4326)
<matplotlib.collections.QuadMesh at 0x1b4259630>
Next, let’s use WaterData service to get HUC8 using
WFS and CQL filter.
layer = "wmadata:huc08" wfs = WFS( ServiceURL().wfs.waterdata, layer=layer, outformat="application/json", version="2.0.0", crs=4269, ) resp = wfs.getfeature_byfilter("huc8 LIKE '13030%'") huc8 = geoutils.json2geodf(resp, 4269, 4326)