Minor Flood Frequency and Duration

2. Minor Flood Frequency and Duration#

../_images/40c0ad9ce275e1931e93aef470ece3518634e90e907067c2e4f9fa2fb14ad025.png

In this notebook we will plot two indicators concerning flooding at the Malakal tide gauge, after first taking a general look at the type of data we are able to plot from the UHSLC. These indicators are based on a ‘flooding’ threshold, using relative sea level.

Download Files: Map | Time Series Plot | Table

Important

This file is taken from https://github.com/jpotemra/PCCM/blob/main/Chapter 5%3A Sea Level/Malakal_SL_FloodFrequency.ipynb to demonstrate how it might work in a juptyer book context. It has been adapted (a lot).

2.1. Setup#

We first need to import the necessary libraries, access the data, and make a quick plot to ensure we will be analyzing the right thing.

2.1.1. Import necessary libraries.#

Next, we’ll establish our directory for saving and reading data, and our output files.

data_dir = Path('../data')
output_dir = Path('../output')  

Then we’ll establish some basic plotting rules for this notebook to keep everything looking uniform.

2.2. Retrieve the Tide Station(s) Data Set(s)#

Next, we’ll access data from the UHSLC. The Malakala tide gauge has the UHSLC ID: 7. We will import the netcdf file into our current data directory, and take a peek at the dataset. We will also import the datums for this location.

and we’ll save a few variables that will come up later for report generation.

2.2.1. Set the Datum to MHHW#

In this example, we will set the datum to MHHW. This can be hard coded, or we can read in the station datum information from UHSLC. I’m not going to do this, because there’s no elegant way to parse the datum table at the moment as far as I can tell. It’s a simple call to the NOAA COOPS API if that’s the data source, though. Instead I have a saved datums_007.csv file that we’ll call from.

# read datums_007.csv
datumtable = pd.read_csv('../data/datums_007.csv', sep=',')

# extract the given datum from the dataframe
datumname = 'MHHW'
datum = datumtable[datumtable['Datum'] == datumname]['Value'].values[0]

# make sure datum is a float64
datum = np.float64(datum)

rsl['datum'] = datum*1000 # convert to mm
rsl['sea_level'] = rsl['sea_level'] - rsl['datum']

# assign units to datum and sea level
rsl['datum'].attrs['units'] = 'mm'
rsl['sea_level'].attrs['units'] = 'mm'

glue("datum", datum, display=False)
glue("datumname", datumname, display=False)

datumtable

Show code cell output Hide code cell output

	Datum	Value	Description
0	Status	14-Nov-2022	Processing Date
1	Epoch	01-Jan-1983 to 31-Dec-2001	Tidal Datum Analysis Period
2	MHHW	2.162	Mean Higher-High Water (m)
3	MHW	2.087	Mean High Water (m)
4	MTL	1.530	Mean Tide Level (m)
5	MSL	1.532	Mean Sea Level (m)
6	DTL	1.458	Mean Diurnal Tide Level (m)
7	MLW	0.974	Mean Low Water (m)
8	MLLW	0.753	Mean Lower-Low Water (m)
9	STND	0.000	Station Datum (m)
10	GT	1.409	Great Diurnal Range (m)
11	MN	1.114	Mean Range of Tide (m)
12	DHQ	0.075	Mean Diurnal High Water Inequality (m)
13	DLQ	0.220	Mean Diurnal Low Water Inequality (m)
14	HWI	Unavailable	Greenwich High Water Interval (in hours)
15	LWI	Unavailable	Greenwich Low Water Interval (in hours)
16	Maximum	2.785	Highest Observed Water Level (m)
17	Max Date & Time	23-Jun-2013 22	Highest Observed Water Level Date and Hour (LST)
18	Minimum	-0.034	Lowest Observed Water Level (m)
19	Min Date & Time	19-Jan-1992 16	Lowest Observed Water Level Date and Hour (LST)
20	HAT	2.567	Highest Astronomical Tide (m)
21	HAT Date & Time	27-Aug-1988 22	HAT Date and Hour (LST)
22	LAT	0.175	Lowest Astronomical Tide (m)
23	LAT Date & Time	31-Dec-1986 16	LAT Date and Hour (LST)
24	LEV	1.611	Switch 1 (m)
25	LEVB	1.529	Switch 2 (m)

2.2.2. Assess Station Data Quality for the POR (1983-2022)#

To do this, we’ll plot all the sea level data to make sure our data looks correct, and then we’ll truncate the data set to the time period of record (POR).

fig, ax = plt.subplots(sharex=True)
fig.autofmt_xdate()
ax.plot(rsl.time.values,rsl.sea_level.T.values/10)

ax.set_ylabel(f'Sea Level (cm, {datumname})') #divide by 10 to convert to cm

2.2.2.1. Identify epoch for the flood frequency analysis#

Now, we’ll calculate trend starting from the beginning of the tidal datum analysis period epoch to the last time processed. The epoch information is given in the datums table.

#get epoch start time from the epoch in the datumtable
epoch_times = datumtable[datumtable['Datum'] == 'Epoch']['Value'].values[0]

#parse epoch times into start time
epoch_start = epoch_times.split(' ')[0]
epoch_start = datetime.datetime.strptime(epoch_start, '%d-%b-%Y')

# and for now, end time the processind end time
epoch_end = datumtable[datumtable['Datum'] == 'Status']['Value'].values[0]

# convert to datetime object
epoch_end = datetime.datetime.strptime(epoch_end, '%d-%b-%Y')

hourly_data = rsl.sel(dict(time=slice(epoch_start, epoch_end)))
hourly_data

glue("startEpochDateTime",epoch_start.strftime('%Y-%m-%d'), display=False)
glue("endEpochDateTime",epoch_end.strftime('%Y-%m-%d'), display=False)

and plot the hourly time series

fig, ax = plt.subplots(sharex=True)
fig.autofmt_xdate()
ax.plot(hourly_data.time.values,hourly_data.sea_level.T.values/10) #divide by 10 to convert to cm

ax.set_ylabel(f'Sea Level (cm, {datumname})')

glue("TS_full_fig",fig,display=False)

../_images/a62cff3b63a956a0e1c1b4781341e98a1d89c356d1726e01ecfefefdb504ee91.png — Fig. 2.1 Full time series at Malakal,Palau tide gauge for the entire record from 1969-05-18 to 2023-10-31. Note that the sea level is plotted in units of cm, relative to MHHW.#

2.4. Calculate the percent change over time#

Next we’ll calculate the percent change over the POR at the tide station/s, for both Frequency and Duration.

The next code cell fits a trend line to the flood days per year data and calculates the trend line. It then plots the data and the trend line in a figure, and calculates the percent change in flood days per year using the trend line. The same process is repeated for flood hours per year.

"""

"""

slope, intercept, r_value, p_value, std_err = stats.linregress(flood_days_per_year['year_storm'], flood_days_per_year['flood_days_count'])
trend = intercept + slope * flood_days_per_year['year_storm']

fig, axs = plt.subplots(2,1)
axs[0].plot(flood_days_per_year['year_storm'], flood_days_per_year['flood_days_count'], label='Data')
axs[0].plot(flood_days_per_year['year_storm'], trend, label='Trend')
axs[0].set_ylabel('Flood Days')
axs[0].legend()
axs[0].set_xticklabels([])

glue("trend_fig", fig, display=False)

percent_change_days = (trend[-1:] - trend[0]) / trend[0] * 100
glue("percent_change_days", percent_change_days, display=False)
percent_change_days

slope, intercept, r_value, p_value, std_err = stats.linregress(flood_hours_per_year['year_storm'], flood_hours_per_year['flood_hours_count'])
trend = intercept + slope * flood_hours_per_year['year_storm']

axs[1].plot(flood_hours_per_year['year_storm'], flood_hours_per_year['flood_hours_count'], label='Data')
axs[1].plot(flood_hours_per_year['year_storm'], trend, label='Trend')
axs[1].set_xlabel('Storm Year')
axs[1].set_ylabel('Flood Hours')
axs[1].legend()

glue("trend_fig", fig, display=False)

percent_change_hours =  (trend[-1:]-trend[0])/trend[0]*100
glue("percent_change_hours", percent_change_hours, display=False)
percent_change_hours

39    113.806306
Name: year_storm, dtype: float64

../_images/47595c0ed9db42210b81c4937d8ad5df85b9bfa262d36b27542979e1b4c0e046.png

Now we’ll generate a table with this information, which will be saved as a .csv in the output directory specified at the top of this notebook.

# make a dataframe with the percent change in flood days and hours per year, with given threshold
percent_change_df = pd.DataFrame({'percent_change_days': percent_change_days.values[0], 'percent_change_hours': percent_change_hours.values[0], 'threshold': threshold}, index=[0])

# add the station name and country
percent_change_df['station'] = station
percent_change_df['country'] = country

# reorder the columns
percent_change_df = percent_change_df[['station', 'country', 'threshold', 'percent_change_days', 'percent_change_hours']]

# Define your attributes
attributes = {
    'station': 'Station name',
    'country': 'Country name',
    'threshold': 'Threshold in cm above MHHW',
    'percent_change_days': 'Percent change in flood days per year',
    'percent_change_hours': 'Percent change in flood hours per year'
}

# Open the file in write mode
with open(output_dir / 'SL_FloodFrequency_percent_change.csv', 'w') as f:
    # Write the attributes as comments
    for column, attribute in attributes.items():
        f.write(f'# {column}: {attribute}\n')

    # Write the DataFrame to the file
    percent_change_df.to_csv(f, index=False)

percent_change_df  

	station	country	threshold	percent_change_days	percent_change_hours
0	Malakal	Palau	30	61.836836	113.806306

2.5. Create a map#

In this next section we’ll create a map with an icon depicting the percent change over the Period of Record at the tide station/s, for both Frequency and Duration.

The following cell defines a function that creates a “zebra stripe” pattern for the map borders.

The following cell defines a function to implement the zebra stripes on a given axis.

The following cell defines a function to add arrow icons denoting percent change.

And here is our final plotting code:

xlims = [128, 138]
ylims = [0, 13]
vmin, vmax = 0,100

# fig, ax, crs,cmap = plot_map(vmin,vmax,xlims,ylims)
crs = ccrs.PlateCarree()
fig, axs = plt.subplots(1, 2, figsize=(10, 6), subplot_kw={'projection': crs})


for i, ax in enumerate(axs):
    ax.set_xlim(xlims)
    ax.set_ylim(ylims)

    ax.coastlines()

    # Fill in water
    ax.add_feature(cfeature.LAND, color='lightgrey')
    # add a) b) labels
    ax.text(0.95, 0.95, f'({chr(97 + i)})', 
            horizontalalignment='right', verticalalignment='top', transform=ax.transAxes,
            fontsize=16)
    
    ax.add_feature(cfeature.OCEAN, color='lightblue')

    # add map grid
    gl = ax.gridlines(draw_labels=True, linestyle=':', color='black',
                  alpha=0.5,xlocs=ax.get_xticks(),ylocs=ax.get_yticks())
    
    #remove labels from top and right axes
    gl.top_labels = False
    gl.right_labels = False

    if ax == axs[1]:
        gl.left_labels = False

add_arrow(axs[0], rsl['lat'], rsl['lon'], percent_change_days, crs, vmin, vmax)
add_arrow(axs[1], rsl['lat'], rsl['lon'], percent_change_hours, crs, vmin, vmax)

# Plot zebra frame
# plot_zebra_frame(axs[0], lw=5, segment_length=2, crs=crs)
# plot_zebra_frame(axs[1], lw=5, segment_length=2, crs=crs)

# Add text for percent change
for i in range(len(rsl['lon'])):
    axs[0].text(rsl['lon'][i] + 0.25, rsl['lat'][i] + 0.25, '{:.1f}%'.format(percent_change_days.values[0]), fontsize=14,
                color=adjusted_heatmap_palette(percent_change_days.values / vmax))
    axs[1].text(rsl['lon'][i] + 0.25, rsl['lat'][i] + 0.25, '{:.1f}%'.format(percent_change_hours.values[0]), fontsize=14,
                color=adjusted_heatmap_palette(percent_change_hours.values / vmax))


glue("mag_fig", fig, display=False)

# Save the figure
output_file_path = output_dir / 'SL_FloodFrequency_map.png'
fig.savefig(output_file_path, dpi=300, bbox_inches='tight')

../_images/87a33c406b547903a5b20fcef3f11ef7dcc17c7e062a162055ebaa036e50cbed.png — Fig. 2.5 Map of the percent change in average flood (a) days and (b) hours per year above a threshold of 30.0 cm above MHHW per year at Malakal,Palau tide gauge from 1983-01-01 to 2022-11-14#

2.6. Citations#

SPM+14: William Sweet, Joseph Park, John Marra, Chris Zervas, and Stephen Gill. Sea level rise and nuisance flood frequency changes around the United States. Technical Report, NOAA, 2014. URL: https://repository.library.noaa.gov/view/noaa/30823.
TWM+19: Philip R. Thompson, Matthew J. Widlansky, Mark A. Merrifield, Janet M. Becker, and John J. Marra. A Statistical Model for Frequency of Coastal Flooding in Honolulu, Hawaii, During the 21st Century. Journal of Geophysical Research: Oceans, 124(4):2787–2802, 2019. URL: https://onlinelibrary.wiley.com/doi/abs/10.1029/2018JC014741, doi:10.1029/2018JC014741.

Minor Flood Frequency and Duration

Contents

2. Minor Flood Frequency and Duration#

2.1. Setup#

2.1.1. Import necessary libraries.#

2.2. Retrieve the Tide Station(s) Data Set(s)#

2.2.1. Set the Datum to MHHW#

2.2.2. Assess Station Data Quality for the POR (1983-2022)#

2.2.2.1. Identify epoch for the flood frequency analysis#

2.2.3. Adjust the data from calendar year to storm year#

2.2.4. Assign a Threshold#

2.3. Calculate and Plot Flood Frequency#

2.3.1. Plot Flood Frequency Counts#

2.3.2. Plot Flood Duration#

2.4. Calculate the percent change over time#

2.5. Create a map#

2.6. Citations#