Sunny Day Flooding

Edmund B Robinson

Workflow Knowledge Domains

The Wolfram Language places a large number of knowledge domains at your fingertips.  This workflow only touches a small fraction of what is available.  Some of these domains are included below in order of appearance.

Data import and navigation

URL, Import, Dataset, Query, ...

Package Development

.wl Files

Geographic visualisation and modeling

DynamicGeoGraphics, GeoElevationData, ...

Wolfram Knowledge Base

Entity Framework and Administrative Divisions

Probability and Statistics

TimeSeriesModelFit, GeneralizedLinearModelFit, NonlinearModelFit, ...

Machine Learning

Predict

Data visualisation and presentation

DateListPlot, Histogram, Plot, Grid, Show, ...

User interface construction

Manipulate and DynamicModule

Sunny Day Flooding?

There are increasing occurrences of tidal heights  exceeding their historic ingress inland.  These occurrences are attributed both the subsidence and change in sea level.

First@StringSplit[WikipediaData["Sunny day flooding"],"\n"]

Many recorded occurrences.

dsSunnySearch=WebSearch["Sunny day flooding",MaxItems->10];
dsSunnySearch[TableForm,
Tooltip[
Hyperlink[#"PageTitle",#"Hyperlink"],
Panel[#"Snippet",BaseStyle->{FontSize->Scaled[.02]}]
]&]

Economical and social impacts.

dsRefugees=WebSearch["USA climate refugees",MaxItems->10];
dsRefugees[TableForm,
Tooltip[
Hyperlink[#"PageTitle",#"Hyperlink"],
Panel[#"Snippet",
BaseStyle->{FontSize->Scaled[.02]}]
]&]

Industry and governmental impacts

Hyperlink["Shell and Exxon's secret 1980s climate change warnings","https://www.theguardian.com/environment/climate-consensus-97-per-cent/2018/sep/19/shell-and-exxons-secret-1980s-climate-change-warnings"]

dsLawSuit=WebSearch["City sues oil industry",MaxItems->10];
dsLawSuit[TableForm,
Tooltip[
Hyperlink[#"PageTitle",#"Hyperlink"],
Panel[#"Snippet",
BaseStyle->{FontSize->Scaled[.02]}]
]&]

Data Discovery and Exploration

What data is needed to explore this issue?

Where can this data be found?

How can it be turned into information?

Set out to locate and interpret this information.

Locating Data

NOAA Tides & Currents

NWS Food Stage Thresholds

Accessing Data

NOAA Data API

NWS Data

Data Wrangling

Tidal Station Meta

Flood Stage Data

Aligning Data Sets

Selecting Tidal Data

Analysis and Communication

Selecting Tidal Data

Modeling Sunny Day Flooding

Observations and Further Analysis

1. Locating Data

1.1   NOAA Tides & Currents

The (HyperlinkURL["National Oceanic and Atmospheric Administration", https://www.noaa.gov/]) (NOAA) maintains a network of tidal stations.

Current station meta data and recent sensor reading can be obtained from their (HyperlinkURL["Tides & Current", https://tidesandcurrents.noaa.gov/]) website.

1.2  NWS Flood Stage Thresholds

The (HyperlinkURL["National Weather Service", https://www.weather.gov/]) (NWS) is responsible for setting flood stage thresholds for NOAA’s tidal stations.

The three threshold types of interest are :

Action Stage: Tidal height when local authorities should assess take action to respond to a flood threat.

Flood Stage: Water is present beyond its normal path of flow.  Height of nuisance flooding.

Moderate Stage: Height of water sufficient to create material damage.

Current station meta data and flood stage thresholds can be obtained from station their pages on their (HyperlinkURL["National Water Observations", https://water.weather.gov/]) website.

2. Accessing Data

2.1 NOAA Data API

NOAA Tides & Currents exposes historical data via its (HyperlinkURL["CO-OPS API For Data Retrieval", https://tidesandcurrents.noaa.gov/api/]).

Various options and some conditions on volume of data returned for each call.

NOAAData Package

Currently only holds the NOAATidesAndCurrents function.

SetDirectory[NotebookDirectory[]];
Needs["NOAAData`"]
ResetDirectory[];

?NOAATidesAndCurrents

NOAATidesAndCurrents[
{
"station"-> "8658120",
"begin_date"->DateObject@{2019,8,1},
"end_date"->DateObject@{2019,8,2},
"product"->"high_low",
"datum"->"MLLW",
"interval"->"h",
"units"->"metric",
"time_zone"->"lst_ldt",
"format"->"json",
"application"->"Wolfram Tech Conference 2019 "<>ToString@$SessionID
}
]

Dataset@NOAAResultToAssociation@%

2. Accessing Data - NWS Data

2.2 NWS Data

The NWS  National Water Observations exposes current flood stage thresholds on each station’s hydrograph page. An example is the (HyperlinkURL["Cape Fear River at Wilmington hydrograph", https://water.weather.gov/ahps2/hydrograph.php?wfo=ilm&gage=wlon7]).  However, there are no active APIs to extract this data.  After a few calls and emails a data file was provided with some of the information needed.

3. Data Wrangling

A large effort always falls to preparing the data for analysis.

3.1 Tidal Station Meta

XML import of tidal station meta data from (HyperlinkURL["Stations and Their Sensors Listing", https://opendap.co-ops.nos.noaa.gov/stations/index.jsp]).

xmlStationMeta=Import["https://opendap.co-ops.nos.noaa.gov/stations/stationsXML.jsp","XML"];

xmlStationMeta[[2,3,1]]

dsStationMeta=stationMetaData/@xmlStationMeta[[2,3]]//Dataset

Water level sensor filter

Not all stations record water level as indicated by the absence of a "Water Level" sensor in the parameters.

dsStationMeta=dsStationMeta[Select[MemberQ[#[["parameter",All,"name"]],"Water Level"]&]]

Costal station locations

dsStationMeta[GeoPosition/*Point/*(DynamicGeoGraphics[{PointSize[Large],Purple,#},ImageSize->Large]&),"metadata","location",{"lat","long"}/*Values/*ToExpression]

3. Data Wrangling

3.2 Flood Stage Data

NWS and NOAA station identifiers do not match.

No published link between NOAA station and NWS gauge even though they are the same location and instruments.

Contacted NOAA and NWS and and received data to link the two.

dsFloodStage=BinaryDeserialize@ByteArray[]

3. Data Wrangling

3.3 Aligning Data Sets

JoinAcross (SQL  JOIN equivalent) meta data and flood stage datasets the common "ID" column and highlight the stations with non-zero "Flood Stage" entries.

JoinAcross[dsStationMeta,dsFloodStage,"ID","Inner"][
(DynamicGeoGraphics[{PointSize[Large],#},ImageSize->Large]&),
({If[#["Flood Stage"]==0,Red,Purple],Point@GeoPosition[ToExpression/@{#["metadata","location","lat"],#["metadata","location","long"]}]}&)
]

3. Data Wrangling

3.4 Selecting Tidal Data

Select Tidal Station

DynamicModule[{selectedGeoPosition={}},
JoinAcross[dsStationMeta,dsFloodStage,"ID","Inner"][
(DynamicGeoGraphics[{
PointSize[.025],Blend[{RGBColor["#CC6600"],Orange}],Dynamic[Point@selectedGeoPosition],
PointSize[Large],#},
GeoProjection->"Equirectangular",
ImageSize->Large]&),
({If[#3==0,Darker@Red,RGBColor["#003366"]],
Tooltip[Button[Point[GeoPosition@{##4}],stationID=#1;selectedGeoPosition=GeoPosition@{##4}],TableForm[{{#1},{#2}},TableHeadings->{{"ID","name"}, None}]]}&
)@@Flatten@{Lookup[{"ID","name","Flood Stage"}]@#,
Query["metadata","location",Lookup[{"lat","long"}]/*ToExpression]@#}&
]]

Dynamic@stationID

Collect Station Data

Some conditions to access and collecting data where discovered.

Need to be mindful of the calls to the API.

There are restrictions on the volume of data each call can return.

Define a methodology to manage the number of calls and volume of data returned.

(HyperlinkURL["Station Datum", https://tidesandcurrents.noaa.gov/datum_options.html])

With[{fn=FileNameJoin[{NotebookDirectory[],"Tide Datums",StringTemplate["`1` - Datum.wdx"][stationID]}]},
If[FileExistsQ[fn],
dsStationDatum = Import[fn,"WDX"]
,
dsStationDatum=
With[{res=NOAATidesAndCurrents[
{
"station"->stationID,
"product"->"datums",
"datum"->"MLLW",
"units"->"metric",
"format"->"json",
"application"->"Wolfram Tech Conference 2019 "<>ToString@$SessionID
}]},
GeneralUtilities`ToAssociations@ImportString[res["Body"],"JSON"]//Dataset
];
Export[fn,dsStationDatum,"WDX"];
dsStationDatum
]]

3. Data Wrangling

3.4 Selecting Tidal Data cont’d.

Historical tidal dataset

With[{fn=FileNameJoin[{NotebookDirectory[],"Tide Datums",StringTemplate["`1` - Tide.wdx"][stationID]}]},
If[FileExistsQ[fn],
(*Import existing*)
dsTide = Import[fn,"WDX"];
(*Append new data*)
dsTide =
Module[{assoc=Normal@dsTide},
assoc[["data"]]=
Join[
assoc[["data"]],
With[{max=dsTide["data",Max,"t"]},
Normal@fetchTideData[stationID,max]["data",Select[#["t"]>max&]]
]
];
assoc//Dataset
]
,
(*Fetch new station data*)
dsTide=fetchTideData[stationID];
];
(*Export for later use*)
Export[fn,dsTide,"WDX"];
dsTide]

What are some higher high tide occurrences?

dsTide["data",Select[#["ty"]=="HH"&]/*TakeLargestBy[#["v"]&,20]]

Further Adjustments

Adjustments can be made for tropical storm events to remove tidal surge series in the data.  There is no API for NOAA Hurricane Center .SHP files.  However, if this data was readily available it would possible to remove data from stations over time periods they were in some distance and region of a storm.  Can see some techniques that would assist this effort in the (HyperlinkURL["Estimating Windstorm Losses", https://www.youtube.com/watch?v=TKpaiXz_law]) talk from WTC 2018.

4. Analysis and Communication

4.1 Flood Stages

What are the flood stage thresholds?

dsFloodStage[SelectFirst[#["ID"]== stationID&]]

floodStage=dsFloodStage[SelectFirst[#["ID"]== stationID&],"Flood Stage"]
actionStage=dsFloodStage[SelectFirst[#["ID"]== stationID&],"Action Stage"]

In previous slide all of the top 20 were above flood stage.  However, this is on 40 years of data.

dsTide["data",Select[#["ty"]=="HH"&]/*TakeLargestBy[#["v"]&,20]]

Action Days

Number of action days?

actionStage=dsFloodStage[SelectFirst[#["ID"]== stationID&],"Action Stage"];
dsTide["data",GroupBy[CurrentDate[#["t"],"Day"]&]/*Select[actionStage≤Max[#[[All,"v"]]]<floodStage&]/*Length]

How many action days per year?

dhActionDaysPerYear=
dsTide["data",
GroupBy[CurrentDate[#["t"],"Day"]&]/*
Select[actionStage≤Max[#[[All,"v"]]]<floodStage&]/*
Keys/*
(DateHistogram[#,"Year",
ChartStyle->LightGray,PlotRange->{{Automatic,DatePlus[NextDate[Max[#1],"Year"],{5,"Year"}]},Automatic},GridLines->{None,Automatic},Frame->True,PlotLabel->"Tidal Action Days Per Year (as at "<>DateString[dsTide["data",Max,"t"],{"DayShort"," ","MonthNameShort"," ","Year"}]<>")",
ImageSize->400,
Epilog->{Inset[Column[{"Flood Stage: NOAA National Weather Service","Tidal Datum: NOAA Tides & Currents"}],Scaled[{0.02`,0.98`}],{-1,1},BaseStyle->{FontSize->Small,FontFamily->"Verdana",FontColor->Darker[Blue]}]}
]&)
]

Flood Days

Number of flood days?

floodStage=dsFloodStage[SelectFirst[#["ID"]== stationID&],"Flood Stage"];
dsTide["data",Select[#["v"]≥floodStage&]/*DeleteDuplicates/*Length,"t"/*(CurrentDate[#,"Day"]&)]

How many flood days per year?

dhFloodDaysPerYear=
dsTide["data",
Select[#["v"]≥floodStage&]/*
DeleteDuplicates/*
(DateHistogram[#,"Year",
ChartStyle->LightGray,PlotRange->{{Automatic,DatePlus[NextDate[Max[#1],"Year"],{5,"Year"}]},Automatic},GridLines->{None,Automatic},Frame->True,PlotLabel->"Tidal Flood Days Per Year (as at "<>DateString[dsTide["data",Max,"t"],{"DayShort"," ","MonthNameShort"," ","Year"}]<>")",ImageSize->Medium,Epilog->{Inset[Column[{"Flood Stage: NOAA National Weather Service","Tidal Datum: NOAA Tides & Currents"}],Scaled[{0.02`,0.98`}],{-1,1},BaseStyle->{FontSize->Small,FontFamily->"Verdana",FontColor->Darker[Blue]}]}
]&),
"t"/*
(CurrentDate[#,"Day"]&)]

When in the year do the flood days occur?

maxDaysPerMonth=
dsTide["data"/*Max,GroupBy[DateObject[DateValue[#["t"],{"Year","Month"}],"Month"]&],Select[#["v"]≥floodStage&]/*DeleteDuplicates/*Length,"t"/*(CurrentDate[#,"Day"]&)]

Is it possible for a month to have 25 days in a flood stage state?  I called NOAA back to enquire and it was confirmed that the area around this station was regularly flooding.

dsFloodDaysPerMonthPlots=
dsTide["data",
GroupBy[DateValue[#["t"],"Year"]&]/*
KeyValueMap[
Function[{key,value},
key->Query[Select[#["v"]≥floodStage&]/*
DeleteDuplicates/*
(floodDaysPerYearMonthChart[#,
PlotLabel->StringTemplate["Tidal Flood Days Per Month (`1`)"][key],
ChartStyle->LightGray,
ImageSize->Medium,
Frame->True]&),
"t"/*
(CurrentDate[#,"Day"]&)]@value
]]/*
Apply[Association]];
laFloodPerMonth=
ListAnimate[
Normal@Values@dsFloodDaysPerMonthPlots,
AnimationRunning->False,
FrameMargins->Tiny,
ImageMargins->None]

4. Analysis and Communication

With data collected and aligned we can now start to look at the data.

4.1 Behavior over time

Several approaches to examining the past behavior and predicting future behaviour.  The TimeSeries object is known to all of these Wolfram Language methods.

tsAnnualFloodDays=dsTide["data",GroupBy[DateObject[{DateValue[#["t"],"Year"]},"Year"]&],Select[#["v"]≥floodStage&]/*DeleteDuplicates/*Length,"t"/*(CurrentDate[#,"Day"]&)]//TimeSeries@*Normal

A useful function by Anton Antonov, soon to be in the Wolfram Function Repository, is QuantileRegression.

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"];
DynamicModule[{qs,qFuncs},
Manipulate[
qs=DeleteDuplicates@Join[Subdivide[0.05,.5,percentileCount],Subdivide[.5,.95,percentileCount]];
qFuncs=QuantileRegression[Most@tsAnnualFloodDays["Path"],knotCount,qs,InterpolationOrder->interpOrder];
Show[
DateListPlot[tsAnnualFloodDays,Joined->False,PlotTheme->"Detailed",GridLines->{False,True}],
Plot[Evaluate[Through@qFuncs@x],{x,tsAnnualFloodDays["FirstTime"],tsAnnualFloodDays["LastTime"]}]
],
{{percentileCount,1,"Curves"},1,3,1,Appearance->"Labeled"},
{{knotCount,3,"Knot Count"},1,5,1,Appearance->"Labeled"},
{{interpOrder,2,"Interpolation Order"},1,4,1,Appearance->"Labeled"},
SynchronousUpdating->False,
ContinuousAction->False,
SaveDefinitions->True
],
SaveDefinitions->True
]

4. Analysis and Communication

4.1 Behavior over time cont’d.

Time Series Analysis / Forecasting

tsmAnnualFloodDays=TimeSeriesModelFit[TimeSeriesWindow[tsAnnualFloodDays,{Automatic,DatePlus[tsAnnualFloodDays["LastDate"],{-1,"Day"}]}],"ARIMA"]

tsaPlot=plotFittedModel[tsmAnnualFloodDays]

4. Analysis and Communication

4.1 Behavior over time cont’d.

Generalised Linear Models

Application of GLM methods.

glmModels=
AssociationThread[{"ModelFunction","Form","Parameters","Options","Label"}->#]&/@
{{GeneralizedLinearModelFit,\[FormalX],{None},ExponentialFamily->"Poisson","GLM[Poisson]"},
{GeneralizedLinearModelFit,\[FormalX],{None},ExponentialFamily->"Gaussian","GLM[Gaussian]"}
}//Dataset

Fit the models.

With[{tsw=
TimeSeriesWindow[
tsAnnualFloodDays,
{Automatic,
DateObject@{DateValue[tsAnnualFloodDays["LastDate"],"Year"]-1}}]},
glmFits=Normal@
glmModels[Association,
#"Label"->#"ModelFunction"[tsw,Sequence@@ReplaceAll[{None}->Nothing]@{#"Form"/.\[FormalT]->tsw["FirstTime"],#"Parameters",\[FormalX],#"Options"}]&
]
];
glmAllPlots= plotFittedModel/@glmFits

Select a model using the  (HyperlinkURL["Akaike information criterion", https://en.wikipedia.org/wiki/Akaike_information_criterion]) to minimises the estimated information loss.

glmAICLossEst=
Module[{aic=#["AIC"]&/@glmFits,min},
min=Min@aic;
(E^((min-#)/2))&/@aic
]
glmSelected=Part[glmFits,Keys@TakeLargest[glmAICLossEst,1]]
glmPlot=glmAllPlots@@Keys@glmSelected

4. Analysis and Communication

4.1 Behavior over time cont’d.

Model to known forms

Sources of known forms can come from internal or external research.

models=
AssociationThread[{"ModelFunction","Form","Parameters","Options","Label"}->#]&/@
{
{NonlinearModelFit,\[FormalA] \[FormalX]+\[FormalB] (\[FormalX]^2)+\[FormalC] ,{\[FormalA],\[FormalB],\[FormalC]},{None},"(x^2)"},
{NonlinearModelFit,\[FormalA] \[FormalX]+\[FormalB] (\[FormalX]^2)+\[FormalC] (\[FormalX]^3)+\[FormalD] ,{\[FormalA],\[FormalB],\[FormalC],\[FormalD]},{None},"(x^3)"},
{NonlinearModelFit,\[FormalA] (E^(\[FormalB]((\[FormalX]/\[FormalT]))))+\[FormalC],{\[FormalA],\[FormalB],\[FormalC]},{None},"(E^x)"} 
}//Dataset

Fit the models.

With[{tsw=
TimeSeriesWindow[
tsAnnualFloodDays,
{Automatic,
DateObject@{DateValue[tsAnnualFloodDays["LastDate"],"Year"]-1}}]},
kmFits=Normal@
models[Association,
#"Label"->#"ModelFunction"[tsw,Sequence@@ReplaceAll[{None}->Nothing]@{#"Form"/.\[FormalT]->tsw["FirstTime"],#"Parameters",\[FormalX],#"Options"}]&
]
];
kmAllPlots= plotFittedModel/@kmFits

Select a model using the  (HyperlinkURL["Akaike information criterion", https://en.wikipedia.org/wiki/Akaike_information_criterion]) to minimises the estimated information loss.

aicLossEst=
Module[{aic=#["AIC"]&/@kmFits,min},
min=Min@aic;
(E^((min-#)/2))&/@aic
]
fmSelected=Part[kmFits,Keys@TakeLargest[aicLossEst,1]]
kmPlot=kmAllPlots@@Keys@fmSelected

4. Analysis and Communication

4.2 Interactive Exhibits

Modelling Summary

Combine with previous elements.

rptSummary=
Grid[{
{ggStationLoc,TabView@{"Time Series"->tsaPlot,"GLM"->glmPlot,"Known Models"->kmPlot}},
{laFloodPerMonth,dhActionDaysPerYear,}
},
BaseStyle->{FontFamily->"Verdana"}]

Sea Level Simulation

GeoElevationData[GeoBoundsRegion[dsStationMeta[SelectFirst[#["ID"]==stationID&],"metadata","location"/*GeoPosition,Lookup[{"lat","long"}]/*ToExpression],Quantity[1,"Kilometers"]],Automatic,"Region",GeoZoomLevel->12]

Does not include gravitational affects.

5. Observations and Further Analysis

Observations

There is an increase over time of action days and flood days

This increase is not observed at all tidal stations.

There is no dominant trend model across stations.

Further Analysis

Increase station coverage of the flood stage threshold dataset

Extend the historical data capture beyond 1980.

Analyse change in higher high tides.

Compare distribution of rate of change in higher high tides, action days, and flood days across the set of tidal stations.

Include subsidence affects.