Mapping the extreme level estimates (A London tidal flood map example)

Introduction 

As an example of how the extreme value analysis, described in the eponymous section, can be used, a step-by-step description of a method to map the 1000-year tidal flood zone in London is the subject of this page. We will run through a process that can be undertaken in a couple of hours and can be done with data found on the internet. The core ingredients we will use are a terrain, a 1000-year level at the mouth of the estuary and a tide table. In fact the code (using ‘R’) runs in less than 20 seconds, so it’s the data manipulation and writing the code that takes the time. For the analysis we will load in four ‘R’ packages: rgdal, UKFE, raster, and rgeos. Each can be loaded into the ‘R’ work environment using library(PackageName) and if not installed yet they require install.packages(“PackageName”). There are numerous ways such mapping can be undertaken. The more common method is the use of hydrodynamic modelling, in which case the free software HECRAS could be used and the necessary data to derive cross sections could be derived from bathymetry or, for a free alternative, nautical charts that can be viewed online (a laborious task). A statistical method could also be used to translate levels up the estuary. Either way, these both require calibration and data that can’t be easily obtained in a couple of hours, hence this expedient method (although reasonable for mapping this indicative flood risk area, as we’ll see). The process is done in three steps: estimate the extreme levels at the mouth of the estuary, calculate the associated levels upstream through London, and then map the levels across the floodplain.

Why Map?

A question we might ask ourselves once we’ve estimated levels or flows associated with a given return period is what area is likely to be flooded by that event? With the answer we can then see what areas need flood risk management consideration and to what extent (accounting for possible depth of flooding). Hence the Environment Agency map 1000-year and 100-year fluvial flood zones and 1000-year and 200-year tidal flood zones across the country, to facilitate flood risk management decisions in the planning process. These flood zones are mapped assuming there are no flood defences so that the area at risk can be determined and the areas benefitting from the defences can be highlighted. The Thames Barrier and associated tidal walls and embankments prevent such flooding, as mapped here, from happening.

The 1000-year level at Southend

We will first estimate the 1000-year level at Southend on Sea, which is at the mouth of the Thames Estuary, and then translate this level up the Estuary and through London. We will then map the level across the London floodplain. The method assumes that there are no flood defences and is therefore only indicative of the flood risk and the areas that need flood risk management consideration (like the Environment Agency’s flood zones). Here we will use ‘R’, but it should be noted that an alternative, which is more user friendly for making pretty maps, is the freely available software QGIS.

At Southend there is a nice long record from which the annual maximum (AM) sample can be extracted. The following script, if pasted into the ‘R’ console creates a data frame with the maximum levels at Southend each year from 1911 to 2017 and we will assume the generalised extreme value distribution to estimate the 1000-year level (this will then be compared to the EA estimates available online):

Levels <- c(3.32, 3.6, 3.44, 3.32, 3.26, 3.44, 3.35, 3.51, 3.35, 3.28, 3.89, 3.23, 3.44, 3.51, 3.44, 3.25, 3.41, 4.15, 3.51, 3.63, 3.63, 3.69, 3.18, 3.48, 3.54, 3.84, 3.41, 4.18, 3.66, 3.75, 3.49, 3.35, 3.84, 3.96, 3.51, 3.29, 3.35, 3.41, 4.24, 3.44, 3.78, 3.68, 4.6, 3.6, 3.72, 3.32, 3.66, 3.6, 3.44, 3.41, 3.9, 3.51, 3.41, 3.41, 4.15, 3.78, 3.81, 3.61, 3.6, 3.49, 3.54, 3.54, 3.92, 3.4, 3.62, 3.8, 3.82, 4.21, 3.55, 3.67, 3.51, 3.48, 3.46, 3.37, 3.48, 3.68, 3.45, 3.85, 3.53, 3.81, 3.43, 3.58, 3.87, 3.77, 3.79, 3.98, 3.79, 3.83, 3.53, 3.73, 3.87, 3.63, 3.68, 3.94, 3.95, 3.77, 3.99, 3.81, 3.74, 3.52, 3.51, 3.54, 4.11, 3.56, 3.55, 3.64, 3.71)

Year <- seq (1911, 2017) #create a sequence of the associated years

SouthendAM <- data.frame(Year, Levels) #combine the levels and dates in a data frame

With this data frame and the following code we can create a plot which should look like figure 1:

plot(SouthendAM, type = "h", main = "Southend annual maximum sea levels", ylab = "Level (mAOD)", xlab = "Year")

One thing to note is that there is a trend in the sea level at Southend (the sea level has been rising over the period of record). This needs to be accounted for in our estimation of the extreme levels (otherwise we would underestimate). Using the UKFE ‘R’ package there are two ways this can be considered. Firstly, we can undertake a trend hypothesis test which provides the probability of the apparent observed trend occurring by chance.

TrendTest(SouthendAM)

The resulting p-value provides the probability of interest and is 0.0000324; i.e. there is a very small probability that the apparent trend is by chance. We will use the GEVAM function to derive our estimate using the trend = TRUE option to account for the trend in the location parameter. We will also detrend the data (to the most recent year) and estimate with the GEVAM again with the default trend = FALSE. Both viable options for accounting for trend and we will use the mean of them as our final estimate.

GEVAM(SouthendAM$Levels, RP = 1000, trend = TRUE) 

Which provides an estimated sea level of 4.95mAOD

GEVAM(DeTrend(SouthendAM$Levels), RP = 1000)

Which provides an estimated sea level of 5.2mAOD and the average of these is 5.08mAOD

We can compare this to the recent estimates from the Environment Agency (EA) which can be found on environment.data.gov.uk site. A shapefile can be downloaded. https://environment.data.gov.uk/DefraDataDownload/?mapService=EA/CoastalDesignSeaLevels&Mode=spatial

The coastal flood boundaries (CFB) shapefile provides points around the coast of the UK, each of which has a range of levels associated with return periods. We can also download a coastal shapefile from the Ordnance survey (OS). The following link provides an option to download numerous ESRI shapefiles including one called ‘coastal’ which we will use. https://osdatahub.os.uk/downloads/open/Strategi

Assuming we have saved the shapefiles to a folder "C:/Users/MyMappingProject” the following code homes in on the data and one of the CFB points can be seen at the end of Southend Pier (the location of the gauge recording the tidal data).

CFB <- readOGR("C:/Users/MyMappingProject/EA_CoastalDesignSeaLevels_SHP_Full/data/Coastal_Design_Sea_Levels_Coastal_Flood_Boundary_Extreme_Sea_Levels.shp")

Coast <- readOGR("C:/Users/MyMappingProject/OSData/strtgi_essh_gb/strtgi_essh_gb/data/coastline.shp")

plot(Coast, xlim = c(514954, 600000), ylim = c(170948, 186600))

plot(CFB, add = TRUE, pch = 19, col = "red")

The British national grid reference of Southend pier can be obtained using the OS’ interactive mapping webpage (link below) and it’s approximately 589069, 183020:

https://osmaps.ordnancesurvey.co.uk

We can use this information to get the data from the CFB shapefile by subsetting it.

CFBSub <- subset(CFB@data, x_bng > 585000 & x_bng < 590000)

CFBSub <- subset(CFBSub, y_bng > 182000 & y_bng < 185000)

CFBSub$t1000

This provides a 1000-year return period level of 5.02mAOD which is only 6cm lower than our central estimate of 5.08mAOD.

The estimate from the previous 2011 coastal flood boundaries was 4.95mAOD (the same as our lower estimate). A much earlier estimate (and there have been many in between) derived during the plans to construct a Thames flood barrier, was 5.16mAOD (Suthons, 1963) which is quoted in ‘The Thames Barrier’ by Gilbert & Horner (1984). It’s very similar to our upper estimate given that we’re using an extra 54 years of data. We can also get an impression of our estimation uncertainty due to sampling error. The below provides a bootstrapped central estimate with 68% intervals (intervals are approximately 4.8 - 5.5mAOD)

UncSS(DeTrend(SouthendAM$Levels), GEVAM, RP = 1000, conf = 0.68) 

The 1000-year levels upriver through London

The requirement is the level-to-level relationship between Southend and sites through London. Therefore, we’ll use the highest astronomical tides (HAT) quoted in the PLA tide tables (see link below) on page 33 (pdf page 35) to provide a difference between the upriver locations and Southend. It should be noted that this method, like hydrodynamic modelling, is deterministic, providing a single outcome in each case, whereas there is a distribution of possible levels at upstream locations for a given level at Southend (air pressure, wind, wave action, and fluvial flow can act to vary from the average difference by up to 0.3m or more)

http://www.pla.co.uk/assets/pla-tide-tables-2020.pdf

Choosing a few locations we can use the last column in the PLA’s table (subtracting the Level of Chart Datum below Ordnance Datum in the second column) to create our own table (table 1).

The second column is the chart datum HAT levels converted to mAOD. The third column provides the difference between Southend and the location associated with each row. The final column provides our Southend 1000-year estimate and adds the difference column to provide the 1000-year levels upstream.

Mapping the flood levels onto the floodplain

To map the levels we’ll first need to create a shapefile of our locations and the associated 1000-year levels in order to associate them with the floodplain.

Using the ordnance survey interactive mapping here:

https://osmaps.ordnancesurvey.co.uk

We can zoom in using the plus minus bar on the right-hand side and find our locations. We can then right click and the British national grid reference is provided. All our locations are in the TQ square so to convert to the coordinate reference we put five in front of the first part of the number and one in front of the second part. Table 2 provides the associated coordinates we’ll use to create a shapefile.

We’ll create a data frame and turn it into a shapefile with the following code.

Easting <- c(589069, 548498, 543186, 533399, 527403, 517004)

Northing <- c(183020, 181621, 179531, 180389, 177520, 175012)

Location <- c("Southend", "Dagenham", "Woolwich", "Tower Pier", "Albert Bridge", "Richmond")

RP1000 <- c(5.08, 5.9, 6.03, 6.08, 6.19, 6.57)

Results <- data.frame(Location, RP1000, Easting, Northing)

#We can now look at our data frame by entering ‘Results’

To turn our resulting levels and locations into a shapefile we need to attribute the coordinates and provide it with a coordinate reference system and projection.

coordinates(Results) <- ~ Easting + Northing

#We can attribute the same projection as our coast shapefile.

proj4string(Results) <- CRS(proj4string(Coast))

We can now plot them. In the same download as the coastal shapefile comes an admin line shapefile, which I read in using readOGR naming it Admin. We’ll add that to the map for some extra context. The plot is zoomed into our London locations, which are highlighted by red dots.

plot(Coast, xlim = c(514954, 557396), ylim = c(170948, 186600))

plot(Admin, add = TRUE)

plot(Results, pch = 19, col = "red", add = TRUE)

Providing figure 3:

The land levels (henceforth called DTM) that we need can be downloaded from the OS website as gridded (raster) data with a resolution of 50 metres:

https://osdatahub.os.uk/downloads/open/Terrain50

Using the raster ‘R’ package we can load in the merged DTM for the TQ British national grid square which covers our area of interest and plot it. I have called the DTM “TQMerged.tif”. First, we’ll set our working directory so we don’t need to put in a file path every time

setwd("C:/Users/MyMappingProject")

TQDEM <- raster("TQMerged.tif")

#We’ll also add the coast shapefile for context

plot(TQDEM)

plot(Coast, add = TRUE)

We’ll crop the DTM so that our mapping is more efficient, using only the area of interest. First we need to create a data frame of the bounding coordinates and turn it into a shapefile.

xbound <- c(514000, 514000, 557350, 557350)

ybound <- c(188800, 166000, 188800, 166000)

Bounds <- data.frame(xbound, ybound)

coordinates(Bounds) <- ~ xbound + ybound

proj4string(Bounds) <- CRS(proj4string(Coast))

CroppedDTM <- crop(TQDEM, Bounds)

plot(CroppedDTM)

The mapping bit

The following code, creates a points shapefile version of the DTM, then measures the distance from each point to each of our locations, it then attributes the closest 1000-year level. Then creates a raster from the point file which means we have a raster of land level and a raster of flood level. We subtract the land level from the flood levels to derive flood depths (attributing anything with a negative value as zero flood depth).

TQGrid <- rasterToPoints(CroppedDTM, spatial = TRUE)

dist <- gDistance(Results, TQGrid, byid = TRUE)

nearest.gauge <- apply(dist, 1, which.min)

TQGrid$nn <- Results$RP1000[nearest.gauge]

TQGridRast <- rasterize(TQGrid, CroppedDTM, "nn")

FloodGrid <- TQGridRast-CroppedDTM

FloodGrid [FloodGrid <= 0] <- 0

We now have a grid of flood depths for the 1000-year tidal flood zone which can be interrogated in any GIS software such as QGIS. We can map it in ‘R’ too:

BreakPoints <- c(0, 0.029, 18.03)

Colours <- c(rgb(0,0,0,0), rgb(0,0.3,0.7, 0.5), rgb(0,0.3,0.7, 0.5))

plot(Coast, xlim = c(514954, 557396), ylim = c(170948, 186600), main = "London 1000-year indicative flood map")

plot(FloodGrid, breaks = BreakPoints, col = Colours, add = TRUE, legend = FALSE)

plot(Admin, add = TRUE)

The map can be refined of course, with some post processing (such as removing “flood islands”) or could be improved with finer resolution DTM (any mapping such as this or hydrodynamic mapping, is sensitive to the DTM) and additional locations. As specified earlier, more data from the tide gauges would also be beneficial to better estimate the extreme levels as a function of Southend and the associated uncertainty. However, for an indicative strategic scale flood map such as this, the flooded area with perhaps a precautionary buffer (based on uncertainty) and filling in some “dry islands” provides an idea of the areas requiring flood risk management (and should be precautionary because the method doesn’t account for friction on the floodplain). We could also re-do the analysis but change the levels using the upper uncertainty to obtain a more precautionary outline, or an uncertainty layer for the one we have. Uncertainty of the upstream levels as a function of Southend could also be considered. Projected Southend levels for future years could be used to see the flood map for the year 2100. Figure 7 provides this by applying the same process as above but with a Southend level of 5.8mAOD (rounded up based on the median RCP8.5 scenario from the UK climate projections 2018 - Palmer, M. et al (2018)). Quite a similar map in that the extents haven’t increased significantly but the gaps are filled in to a greater extent.

We can compare these maps to similar ones on the internet, although the majority include fluvial flood risk too. For another fun comparison, figure 8 provides a tidal flood zone map from the 1984 book ‘The Thames Barrier’ (photographed).

References:

Environment Agency (2019). Coastal flood boundary conditions for the UK: update 2018. SC060064/TR6. Bristol

Gilbert. S, Horner. R, (1984). The Thames Barrier. Thomas Telford Limited. London

Palmer, M., Howard, T., Tinker, J., et al., (2018) UKCP18 Marine Projections: Science Report. Met Office Hadley Centre, Exeter, UK.

Suthons. C, (1963). Frequency of the occurrence of abnormally high sea levels on the east and south coast of England. Proc. Instn Civ. Engrs, 25, 433-440