mappER - Thiessen Polygons • mappER

Overview

Currently the mappER only has the functionality to support Thiessen (Voronoi) polygons. To access this functionality you can load flode or mappER directly.

library(mappER)

library(flode)

Loading gauges

The best way to get gauge data into mappER is to use the riskyData functions such as loadAPI().

 library(riskyData)
 crowle <- loadAPI(ID = '457592',
                  measure = 'rainfall',
                  period = 900, 
                  type = 'total',
                  datapoints = 'earliest')

To save time we will use the pre-loaded data from the package;

data(crowle); data(bickley); data(barnhurst); data(hollies); data(ledbury); 
data(bettwsYCrwyn)

To compile these sites into a spatial dataset we can use the getCoords() function;

gauges <- getCoords(crowle,
                    bickley,
                    barnhurst,
                    hollies,
                    ledbury,
                    bettwsYCrwyn)
gauges
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -3.172591 ymin: 52.0317 xmax: -2.097173 ymax: 52.7993
#> Geodetic CRS:  WGS 84
#>      stationName  WISKI Easting Northing                   geometry
#> 1         Crowle 457592  393455   255757     POINT (-2.097173 52.2)
#> 2        Bickley 445031  363130   271330 POINT (-2.542565 52.33879)
#> 3      Barnhurst 091266  389910   301620  POINT (-2.15044 52.61225)
#> 4        Hollies 091862  381586   322452  POINT (-2.274545 52.7993)
#> 5        Ledbury 459793  370233   237123  POINT (-2.435301 52.0317)
#> 6 Bettws-Y-Crwyn 441011  320360   281360  POINT (-3.172591 52.4244)

If you wished to plot these data, you can do so using leaflet, an open-source JavaScript library for interactive maps.

library(leaflet)
#> Warning: package 'leaflet' was built under R version 4.2.3
leaflet(gauges) %>%
  addTiles() %>%
  addCircleMarkers(color = "orangered",
                   radius = 6,
                   fillOpacity = 0.8,
                   stroke = "black",
                   label = ~stationName,
                   labelOptions = labelOptions(permanent = TRUE,
                                               noHide = TRUE,
                                               direction = "bottom"))

Developing the Thiessen plots

A Thiessen polygon is a Voronoi Diagram that is also referred to as the Dirichlet Tessellation. Given a set of points, it defines a region around each point. A Thiessen polygon divides a plane such that each point is enclosed within a polygon and assigns the area to a point in the point set.

In forecasting we use these to determine gauge proportion. With mappER they are incredibly easy to develop. We can use the teeSun() function with the above gauge data to develop the tessellation.

voronoi <- teeSun(gauges)

Intersecting the Thiessen with a catchment

To load a catchment shapefile we can use the loadCatchment() funtion.

loadCatchment(filepath = "c://abcd/efg/hi/Catchment.shp")

Now we have the catchment loaded, it can be intersected by the Thiessen polygons. This can be done using the intersectPoly() function.

int <- intersectPoly(coords = gauges,
                     voronoi = voronoi,
                     catchment = bewdCatch)

To plot with a custom colour scheme;

colors <- colorFactor(palette = c("blue", "green", "yellow", "red"),
                      domain = (int$stationName))
leaflet(int) %>%
  addTiles() %>%
  addPolygons(color = colors(int$stationName))

We can obtain more details on the gauges when we use the gaugeProp() function, this displays the proportion of catchment that each rain gauge covers. It uses the intersected polygon as the basis for its measurements.

gaugeProp(int)

Station	WISKI	Area	Proportion
Bettws-Y-Crwyn	441011	2465.18	57.127827
Hollies	091862	1050.89	24.353217
Bickley	445031	632.26	14.651928
Barnhurst	091266	166.87	3.867028

These outputs can then be applied in riskyData to generate catchment averaged rainfall.

Putting it all together

To plot all the elements simultaneously we can use the teeSunPlot() function;

teeSunPlot(gaugeCoords = gauges, intersectPoly = int)