Opportunity Analysis Methods
Overview
This document describes the methods used to create spatial input layers for the Tampa Bay Coastal Master Plan (TBCMP) habitat opportunity analysis. All processing is implemented in R/funcs.R and executed in R/01_inputs.R.
Inputs are created for seven counties in the Tampa Bay region: Citrus, Hernando, Hillsborough, Manatee, Pasco, Pinellas, and Sarasota. County boundaries are derived from the 2025 TIGER/Line census subdivision files and projected to EPSG:3087 (NAD83(HARN) / Florida GDL Albers, meters) for all spatial analysis.
Coastal Stratum
The coastal stratum defines the zone of interest for reservation space in the Tampa Bay Habitat Master Plan. It represents the land area between Mean Lower Low Water (MLLW) and 5 feet (1.524 m) NAVD88 most likely to be impacted by sea level rise.
The delineation follows three main steps:
- Load a topobathymetric DEM covering the study area
- Derive a spatially varying MLLW offset surface using the NOAA VDatum API
- Mask the DEM to cells within the stratum bounds and convert to a polygon
Digital Elevation Model
Source
The DEM used is the NOAA Continuously Updated DEM (CUDEM), accessed from the TBCMP ESRI FileGeodatabase:
T:/05_GIS/TBEP/TBCMP/TBCMP_TBEP_DATASETS.gdb :: NOAA_CUDEM_tbcmp_extent_mCUDEM is a topobathymetric product that seamlessly integrates LiDAR-derived land elevation with bathymetric soundings in nearshore areas. This makes it suitable for delineating the intertidal zone. The vertical datum is NAVD88, units are meters.
Processing
The raw CUDEM layer is 2 m resolution and covers the full TBCMP seven-county extent. To make it tractable for analysis:
- Crop and mask to the seven-county study area boundary
- Aggregate by a factor of 5 (mean) to produce a 10 m resolution raster
The processed raster is saved as a Cloud Optimized GeoTIFF (COG) and hosted on Amazon S3 for reproducible access without requiring the source GDB:
https://tbcmp.s3.amazonaws.com/cudem_3087.tifTidal Datum Offset Surface
Rationale
The lower bound of the coastal stratum is MLLW, not a fixed NAVD88 elevation. MLLW varies spatially due to differences in tidal range, bay geometry, and freshwater influence. Across the seven-county study area, MLLW ranges from roughly 0.12 m to −0.75 m NAVD88. Using a spatially varying surface rather than a single constant ensures the stratum boundary is physically meaningful everywhere.
VDatum API
Tidal datum offsets are obtained from the NOAA VDatum web API, which converts between vertical datums at arbitrary geographic locations. For each sample point, the API converts 0 m MLLW to its equivalent NAVD88 elevation, providing the local MLLW offset.
API endpoint:
https://vdatum.noaa.gov/vdatumweb/api/convertKey parameters:
| Parameter | Value |
|---|---|
| Source vertical datum | MLLW |
| Target vertical datum | NAVD88 |
| Horizontal frame | NAD83_2011 |
| Geoid | GEOID18 |
| Source value | 0 m |
Points outside VDatum coverage (e.g., far offshore or inland) return NOAA’s sentinel value (−999999) and are treated as NA.
Sample Grid
A regular grid of sample points at ~5 km spacing (0.05° in geographic coordinates) is generated over the study area bounding box and clipped to the county boundaries. Requests are sent at ~2 per second to stay within NOAA’s API rate limits.
Interpolation
Valid offset values are interpolated to a continuous surface using a thin-plate spline (TPS) from the fields package. TPS is appropriate here because:
- Tidal datum variation is smooth and gradual
- It exactly interpolates at observed points
- It extrapolates smoothly to unsampled areas within the study extent
The spline is predicted onto a 0.005° (~500 m) grid in EPSG:4326, then projected to EPSG:3087 and resampled to exactly match the DEM grid using bilinear interpolation.
Coastal Stratum Mask
The stratum is defined as all DEM cells satisfying:
\[\text{MLLW}_{xy} \leq \text{elevation}_{xy} \leq 1.524 \text{ m NAVD88}\]
where \(\text{MLLW}_{xy}\) is the spatially varying MLLW surface and 1.524 m is the 5 ft NAVD88 upper bound. Cells outside this range are set to NA.
The raster mask is then vectorized (dissolved) to a single polygon using terra::as.polygons() and validated with sf::st_make_valid().
Output attributes
| Field | Value |
|---|---|
stratum |
"coastal" |
lower_datum |
"MLLW (spatially varying, NAVD88)" |
upper_elev_ft |
5 |
upper_elev_m |
1.524 |
crs |
"EPSG:3087" |
The final polygon is saved to data/coastal_stratum.RData.
Key Functions
| Function | Purpose |
|---|---|
query_vdatum() |
Single-point VDatum API call |
build_vdatum_grid() |
Generate clipped sample grid |
batch_query_vdatum() |
Loop over grid with rate limiting |
build_mllw_surface() |
TPS interpolation + resample to DEM |
make_coastal_stratum() |
Mask DEM and polygonize |
Soils Classification
The soils layer classifies SSURGO map units across the seven-county study area into three hydrologic categories — Xeric, Mesic, and Hydric (Ries and Scheda 2014). These categories reflect the soil moisture regime and are used to characterize habitat restoration potential.
Source Data
Soil map unit polygons are obtained from the NRCS SSURGO database via the Soil Data Access (SDA) web service, queried using the soilDB R package. Polygons are fetched by county to avoid API timeouts, with each county geometry subdivided into a regular tile grid (~0.1° tiles) and queried tile by tile.
Classification Method
Map units are classified by matching the SSURGO map unit name (muname) against lists of Florida soil series associated with each hydrologic category. Matching uses case-insensitive regex with word-boundary anchors to avoid partial string collisions.
The classification is applied in the following order of precedence:
Xeric (gridcode 100) — matched first. Well-drained, upland sandy soils (e.g., Tavares, Astatula, Candler, Paola). Xeric is evaluated before Mesic so that complex map unit names with a leading Xeric series (e.g., “Tavares-Myakka complex”) resolve correctly.
Hydric keyword override (gridcode 300) — map unit names containing the terms
muck,depressional, ormostly wetlandare classified as Hydric regardless of series name. This handles series such as Pompano and Basinger that appear in both Mesic and Hydric contexts.Hydric by series (gridcode 300) — organic and permanently saturated mineral soils (e.g., Sellers, Kesson, Wulfert, Terra Ceia, Samsula, Weekiwachee, Hydraquents).
Mesic (gridcode 200) — seasonally saturated flatwoods soils that are wet but not permanently inundated (e.g., Myakka, Immokalee, Pineda, Wauchula, EauGallie, Basinger).
Map units not matching any series list (primarily water bodies) receive NA and are excluded from the output.
Output
Classified polygons are clipped to the study area boundary and dissolved by gridcode, producing three polygons representing the spatial extent of each soil type.
| gridcode | Descrip |
|---|---|
| 100 | Xeric |
| 200 | Mesic |
| 300 | Hydric |
The final layer is saved to data/soils.RData.
Key Functions
| Function | Purpose |
|---|---|
fetch_soils_tiled() |
Fetch SSURGO polygons via tiled SDA queries |
classify_soils_muname() |
Query munames from SDA and classify by series |
build_soils_layer() |
Combine fetch, classify, clip, and dissolve |
Salinity
The salinity layer classifies the study area into three salinity zones based on long-term mean surface water salinity. The zones are used to characterize estuarine habitat type and inform restoration potential.
Source Data
Salinity observations are retrieved from the USEPA Water Quality Portal (WQP). Data are queried for the characteristic name "Salinity" over a multi-year period using the WQP REST API:
https://www.waterqualitydata.us/data/Result/searchQueries are issued county by county using WQP county codes to limit request sizes. Results are cached locally as zip archives in data-raw/wqp_cache.
Units psu, PSU, ppt, and ppth are treated as equivalent at estuarine salinities and pooled before computing station means.
Station Aggregation
For each monitoring station, the long-term mean salinity is computed across all observations in the query period. Stations are converted to an sf point layer in EPSG:3087.
Spatial Interpolation
Station means are interpolated to a continuous surface using a thin-plate spline (TPS) from the fields package, fitted in geographic coordinates (EPSG:4326). The spline is predicted onto a regular grid (default spacing 0.001°, ~100 m), limited to a minimum of 0 psu, then projected to EPSG:3087 using bilinear resampling.
Classification
The continuous salinity surface is reclassified into three zones following the original Tampa Bay Habitat Master Plan salinity layer schema and the 18 psu cutoff of Eleuterius and Eleuterius (1979):
| Class | Descrip | Range |
|---|---|---|
| 0 | Fresh | < 0.5 psu |
| 1 | Oligohaline/Mesohaline | 0.5 – 18 psu |
| 3 | Polyhaline/Euhaline | ≥ 18 psu |
Each class is vectorized from the classified raster using terra::as.polygons() (dissolved), clipped to the county boundary, and unioned to a single MULTIPOLYGON feature.
Boundary Smoothing
Raster-to-vector conversion produces staircase artifacts along zone boundaries. To remove these, the combined three-class polygon layer is simplified using rmapshaper::ms_simplify() with topology preservation enabled. This method applies the Visvalingam-Whyatt algorithm across all features simultaneously, ensuring that shared boundaries between adjacent zones are continuous.
The simplify_keep parameter (default 0.05) controls the fraction of vertices retained. Lower values produce smoother boundaries and higher values stay closer to the original raster edges.
Output
The final layer contains one feature per salinity class with the actual minimum and maximum predicted salinity values within each zone.
| Field | Description |
|---|---|
Classes |
Integer class code (0, 1, 3) |
Value_Min |
Minimum predicted salinity within zone (psu) |
Value_Max |
Maximum predicted salinity within zone (psu) |
Descrip |
Zone label |
crs |
EPSG:3087 |
The final layer is saved to data/salinity_layer.RData.
Key Functions
| Function | Purpose |
|---|---|
build_salinity_layer() |
Fetch WQP data, interpolate, classify, and vectorize |
wqp_download_chunk() |
Download and cache a single WQP query chunk |
wqp_read_zip_csv() |
Extract and parse columns from a WQP zip archive |