NCG618 Assignement n°2

Author

Quentin Mangold, Erasmus student (ID : 24120677)

Published

December 11, 2024

Median rainfall in January in Ireland

The goal of this publication is to provide a map of rainfall in Ireland for the 25 weather stations, colour coding the symbol for each station according to its median rainfall level in January. To do that, we will need first some libraries from R :

library(tidyverse)
library(sp)
library(RColorBrewer)
library(tmap)
library(sf)
  • tidyverse: A collection of R packages for data manipulation, visualization, and analysis. We will use the function summarise from this library.
  • sp: Provides classes and methods for spatial data processing, such as handling points, lines, and polygons in spatial datasets.
  • RColorBrewer: A package for creating visually appealing color palettes.
  • tmap: A package for creating thematic maps, allowing both static and interactive visualizations of spatial data.
  • sf: Provides modern tools for working with simple features (spatial vector data).

A brief presentation of the data

We load our data set rainfall. This document includes two objects : stations and rain. The object stations provides the coordinates, the county, the abbreviation, etc.. of 25 weather stations in Ireland. The object rain provides the rainfall for each month and each station since 1850. Here you can have a look at their first values:

load("C:/Users/33642/Documents/rainfall.RData")
head(stations)
# A tibble: 6 × 9
  Station     Elevation Easting Northing   Lat  Long County  Abbreviation Source
  <chr>           <int>   <dbl>    <dbl> <dbl> <dbl> <chr>   <chr>        <chr> 
1 Athboy             87  270400   261700  53.6 -6.93 Meath   AB           Met E…
2 Foulksmills        71  284100   118400  52.3 -6.77 Wexford F            Met E…
3 Mullingar         112  241780   247765  53.5 -7.37 Westme… M            Met E…
4 Portlaw             8  246600   115200  52.3 -7.31 Waterf… P            Met E…
5 Rathdrum          131  319700   186000  52.9 -6.22 Wicklow RD           Met E…
6 Strokestown        49  194500   279100  53.8 -8.1  Roscom… S            Met E…
head(rain)
# A tibble: 6 × 4
   Year Month Rainfall Station
  <dbl> <fct>    <dbl> <chr>  
1  1850 Jan      169   Ardara 
2  1851 Jan      236.  Ardara 
3  1852 Jan      250.  Ardara 
4  1853 Jan      209.  Ardara 
5  1854 Jan      188.  Ardara 
6  1855 Jan       32.3 Ardara

Then, we can plot each station on an Ireland map.
First, we convert the stations data frame into a spatial object (sf class) using the sf package. The coords argument specifies the columns (Long for longitude and Lat for latitude) used as spatial coordinates. The crs = 4326 argument sets the coordinate reference system (CRS) to WGS84.

stations_sf <- stations %>% st_as_sf(coords=c('Long','Lat'),crs=4326)

Then, we activate interactive viewing mode in the tmap package with the function tmap_mode(), allowing for dynamic exploration of the map.

tmap_mode('view')

Finally, we create a map using the spatial data (stations_sf): tm_shape(stations_sf) specifies the spatial object to visualize and tm_bubbles() adds bubble markers to represent the points. We see on the map that the stations are well settled. We can basically click on a station to get all of its information.

tm_shape(stations_sf) + tm_bubbles()

Now that our coordinates are good, we can set up our data.

Filtering the data

So now, we need to filter our data to keep only the January values. Then we calculate the median rainfall (median(Rainfall)) for each Station in the filtered dataset. The .by argument groups the data by the Station column before performing the summarization. This result is stored in the data frame rain_summary.

rain %>% filter(Month=='Jan') %>% summarise(median_rain=median(Rainfall),.by=Station) -> rain_summary

Now that we have our filtered data, we need to link it to the coordinates of each station. To do so, we use a left join. The left join is a join between two tables that returns all records from the left table, even if there is no match with the right table. If there is no match, the missing values are set to NULL. We put all of these information in the data set station_median_sf.

stations_sf %>% left_join(rain_summary) -> station_medians_sf
head(station_medians_sf[, c(1:4, ncol(station_medians_sf))])
Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -8.1 ymin: 52.28 xmax: -6.22 ymax: 53.75
Geodetic CRS:  WGS 84
# A tibble: 6 × 6
  Station     Elevation Easting Northing median_rain      geometry
  <chr>           <int>   <dbl>    <dbl>       <dbl>   <POINT [°]>
1 Athboy             87  270400   261700        87.1  (-6.93 53.6)
2 Foulksmills        71  284100   118400       106.   (-6.77 52.3)
3 Mullingar         112  241780   247765        80.6 (-7.37 53.47)
4 Portlaw             8  246600   115200       123.  (-7.31 52.28)
5 Rathdrum          131  319700   186000       127.  (-6.22 52.91)
6 Strokestown        49  194500   279100       102.   (-8.1 53.75)

The tables are well linked.

Ploting the map

Once we have our data settled, we can create our map. The col argument in tm_bubbles() sets how the points are colored. We also set the colour in blue gradient palette = "Blues", the size at 1 and we put a title to the legend title.col = "Median rain levels". tm_layout() is used to put a name to the map and to specify that we want the legend outside the map :

tm_shape(station_medians_sf) + 
  tm_bubbles(col = 'median_rain', 
             palette = "Blues", # blue bitmap    
             size = 1, 
             title.col = "Median rainfall levels" #change the title of the legend
             ) +
  tm_layout(
    title = "Median rainfall levels in January for each station",
    legend.outside = TRUE        # legend outside the map
  )

This map is interesting.
Stations along the southwest coast generally show higher median rainfall (177.7 mm for Killarney and 166 mm for Valentia). This trend could suggest that the southwest is wetter compared to other regions.
The east coast and especially Dublin stations (67.6 mm at the Phoenix Park and 63 mm at Dublin Airport) have lower rainfall levels compared to the west. This suggests that the west tends to receive more rain due to its exposure to Atlantic weather systems.
To be more accurate and have a better perception, we can plot a field map to “fill the gap” in between the stations.

To go further : using field map

So basically, we will need two more libraries :

library(fields)
library(raster)
  • fields: Provides tools for spatial data analysis, including functions for interpolation, spatial smoothing, and working with thin plate splines.
  • raster: A package to handle pixel-based geographical objects.

We will also need the Ireland coast document.

coast <- st_read('ireland_coast.geojson', quiet=TRUE)

We create a spatial smoothing model using the Thin Plate Spline (TPS) method provided by the fields package in R. We also use our previous data set where we have all of the medians linked to the stations. The REML method stands for Restricted Maximum Likelihood. This method is often used for robust estimation of model parameters, particularly in spatial contexts like here.

station_medians_sf2 <- st_transform(station_medians_sf, crs = 29902)
tps_model <- Tps(st_coordinates(station_medians_sf2), station_medians_sf$median_rain, method='REML')

To create our prediction, we use the same code as in lectures :
- raster: A package designed to work with maps represented as grids of pixels.
- extent: Creates a rectangular boundary (bounding box) for the grid, based on the area covered by the coast object.
- resolution: Defines the size of each pixel in the grid (in this case, 1 km wide).
- interpolate: Uses the TPS model to calculate predicted values for the center of every pixel in the grid.
- crs: Assigns a coordinate reference system (CRS) to the raster.
- rasterize: Converts the coast object (originally a vector map) into a pixel-based map that matches the grid.
- mask: Applies the rasterized coast as a filter to remove pixels outside the land area.

rex <- extent(st_bbox(coast))
cover_grid <- raster(rex, resolution=1000)
rain_grid <- interpolate(cover_grid, tps_model)
crs(rain_grid) <- CRS("+init=epsg:29902")
m <- rasterize(coast, rain_grid)
rain_grid <- mask(rain_grid, m)

And then, we plot our field map. tm_shape(rain_grid) specifies the spatial object to visualize: here rain_grid is the object we want to visualize. tm_raster() function is used to visualize raster data. It takes rain_grid (the raster object) and displays it with specific visual styling: style='cont'sets the style for displaying continuous data, then we set the title and the colour is still blue. We also take the coastline of Ireland with tm_shape(coast) and tm_borders()adds the borders of the coast object to the map.

tmap_mode('plot')
tm_shape(rain_grid) +
  tm_raster(
    style='cont',
    title='Median rainfall levels in January (prediction)',
    palette = "Blues") + 
  tm_shape(coast) + 
  tm_borders()

The west and southwest coasts show significantly higher rainfall levels, and this is consistent with the trends observed from the stations.
Rainfall levels appear to decrease steadily as we move inland and towards the east coast, with the eastern region receiving less rain.
Note: regions in central Ireland may have less precise predictions due to a lack of nearby stations. The accuracy of predictions in these areas could be improved with more data points.

Conclusion

To conclude, we can see in both maps that the median rainfall level in January can be linked to a geographical aspect. We see that levels are higher on the west and south coast. But, as go more to the north and the east of the island, the rainfall levels decrease (mainly when we go to the east). This suggests that, regions without nearby land, such as the south and west of Ireland directly exposed to the Atlantic, are more prone to heavy rainfall. In contrast, the northeast regions, sheltered by Great Britain, benefit from a barrier effect that reduces their exposure to rain.