<jupyter_start><jupyter_text>Web Services<jupyter_code>from pygeoogc import ArcGISRESTful, WFS, WMS, ServiceURL
import pygeoutils as geoutils
from pynhd import NLDI
import numpy as np<jupyter_output><empty_output><jupyter_text>[PyGeoOGC](https://github.com/cheginit/pygeoogc) and [PyGeoUtils](https://github.com/cheginit/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 a watershed geometry using NLDI and use it for subsetting the data.<jupyter_code>basin_geom = NLDI().getfeature_byid(
    "nwissite",
    "USGS-11092450",
    basin=True
)<jupyter_output><empty_output><jupyter_text>PyGeoOGC has ``NamedTuple`` called ``ServiceURL`` that contains URLs of the some of the popular web services. Let's use it to access Water Boundary Dataset RESTful service and get the HUC 12s within our basin and use ``pygeoutils.json2geodf`` to convert it into a GeoDataFrame.<jupyter_code>wbd12 = ArcGISRESTful(f"{ServiceURL().restful.wbd}/6")
wbd12.get_featureids(basin_geom.geometry[0])
resp = wbd12.get_features()
huc12 = geoutils.json2geodf(resp)
huc12.plot(figsize=(8, 8));<jupyter_output><empty_output><jupyter_text>Next, let's get wetlands using the [National Wetlands Inventory](https://www.fws.gov/wetlands/) WMS service. First we need to connect to the service using ``WMS`` class.<jupyter_code>wms = WMS(
    ServiceURL().wms.fws,
    layers="0",
    outformat="image/tiff",
    crs="epsg:3857",
)
wms<jupyter_output><empty_output><jupyter_text>Then we can get the data using the ``getmap_bybox`` function. Note that this function only accpets a bounding box, so we need to pass a bounding box and mask the returned data later on usin ``pygeoutils.gtiff2xarray`` function.<jupyter_code>r_dict = wms.getmap_bybox(
    basin_geom.geometry[0].bounds,
    1e3,
    box_crs="epsg:4326",
)
wetlands = geoutils.gtiff2xarray(r_dict, basin_geom.geometry[0], "epsg:4326")<jupyter_output><empty_output><jupyter_text>There are three bands in the returned data. Using the information provided for the layer [here](https://www.fws.gov/wetlands/arcgis/rest/services/Wetlands/MapServer/0) we can find the color code of emergent wetlands and extract them.<jupyter_code>emergent = np.array([127, 195, 28])
wl_emergent = wetlands.where((wetlands.isel(band=0) == emergent[0])
                        & (wetlands.isel(band=1) == emergent[1])
                        & (wetlands.isel(band=2) == emergent[2]))<jupyter_output><empty_output><jupyter_text>For plotting the wetlands note that in our request we set the CRS to EPSG:3857 so we need transform the watershed geometry to this CRS so they can be superimposed.<jupyter_code>basin_3857 = basin_geom.to_crs("epsg:3857")
ax = basin_3857.plot(facecolor="none", edgecolor="k", lw=0.5, figsize=(8, 8))
wl_emergent.isel(band=0).plot(ax=ax, add_colorbar=False);<jupyter_output><empty_output><jupyter_text>Next, let's getthe base flood elevation data from [FEMA](https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHL/MapServer/16) and use ``pygeoutils.json2geodf`` to convert it into a GeoDataFrame.<jupyter_code>wfs = WFS(
    ServiceURL().wfs.fema,
    layer="public_NFHL:Base_Flood_Elevations",
    outformat="esrigeojson",
    crs="epsg:4269",
)
wfs
r = wfs.getfeature_bybox(basin_geom.geometry[0].bounds, box_crs="epsg:4326")
flood = geoutils.json2geodf(r.json(), "epsg:4269", "epsg:4326")
flood.plot(figsize=(8, 8));<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup