In RobinKohrs/rainfallR: Extract and manipulate gridded rainfall measurements

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
) 

Load the packages

• First we load some of the necessary packages

library(tidyverse) # For much of the data preprocessing
library(grid)
library(rainfallR) # The access and processing of the rainfall data
library(here) # Project paths
library(gganimate) # For a little animation
library(iffitoR) # The landslide data
library(sf) # Handling of the spatial data
library(glue) # Concattening strings
library(purrr) # Mapping functions
library(crayon)  # Coloured consolde output 
library(raster) # Raster data handling

Get the slope units

• In order to find the slope units on the hard drive in differs between the Linux and the Microsoft path

# which os to automatically set the paths
os = Sys.info()["sysname"]

if(os == "Linux"){
  path_ncdf = "/mnt/CEPH_PROJECTS/Proslide/PREC_GRIDS_updated/"
  su_path = "/mnt/CEPH_PROJECTS/Proslide/Envir_data/SlopeUnits/su_opt_16_TAA/su_16_TAA.shp"
}else if(os == "Windows"){
  path_ncdf = "\\\\projectdata.eurac.edu/projects/Proslide/PREC_GRIDS_updated/"
  su_path = "\\\\projectdata.eurac.edu/projects/Proslide/Envir_data/SlopeUnits/su_opt_16_TAA/su_16_TAA.shp"
}else{
  stop(call. = F, "what the hell are you working on...")
}

What have we got

• Lets understand a bit about the data that we got

• The slope units not only cover South Tyrol, but also the trentino Region. So lets mask it

• We can use the iffitoR::get_shape_southtyrol-function in order to get the outline of South Tyrol

# get vector of South Tyrol
st = iffitoR::get_shape_southtyrol() 

ggplot(st) +
  geom_sf() +
  theme_minimal()

• Now we get the vector data for the slope units

su_orig = st_read(su_path)
# and verify directly if they are in the same crs
st_crs(su_orig) == st_crs(st)

• As they are not in the same crs, lets reproject the ouline of South Tyrol. The difference, to my knowledge lies in the different geodetic reference systems they use (ETRS89 in the case of South Tyrol and GRS 80 for the slope units).

st_reproj = st_transform(st,st_crs(su_orig))
st_crs(st_reproj) == st_crs(su_orig)

Clip the slopunits to South Tyrol

• In order to save some computing time we clip the extent of the slope units to the extent of South Tyrol

su = st_crop(su_orig, st_reproj)
# are the two bounding boxes equal?
st_bbox(su) == st_bbox(st_reproj) # for some reason the xmax is not the same

• How many slope units are left now?

dim(su_orig) # how many were there
dim(su)

Get the rainfall for the slope Units

• For each Slope Unit we will now get the dates with landslides in them

• We can use the function slide_dates_in_polygon in order to extract all the points (landslides-initiation-points) that fall in each polygon (slope unit)

• As we are using the landsld-object we are considering landslides of all types. Here is where we should apply some kind of filtering in order to select the slides that are of intereset to us.

• Each polygon (could be catchements, slope units, ...) will have a column in the output that is called poly_id. Each observation (each row) is uniquely identifiable by the date and the poly_id.

• If there are more than one row for one poly_id this means that in the same polygon on several dates slides happened

• We thus only extract the rainfall once for each polygon that saw a movement on a date

• In order to not neglect the fact that potentially multiple slides could have happened, the column slides_per_poly_date present the number of slides that happened at that day in that polygon

landsld = landsld %>% filter(str_detect(second_level, "translational|rotational")) %>% slice(1:50)

su_points_in_poly = rainfallR::slide_dates_in_polygon(poly = su_path
                                                      ,landsld = landsld)

• Now we have all the dates of all slides that happened within each polygon of su

• One interesting question would be to ask what the day and slope unit was that saw most movements?

su_points_in_poly %>% 
  dplyr::select(c(poly_id, date, month.int, second_level, slides_per_poly_date)) %>%  
  arrange(desc(slides_per_poly_date)) %>% 
  head()

Extract the rainfall data for each su with landslides for each day

• We only are interested in the rainfall of the slope units that actually experienced landslides. Therefore we loop over all the slope units and only take the ones with slides.

What about the rainfall Events that happened in slope units and did not cause landslides?

• Then we look in the dates column and extract the rainfall for each of the date for that slope unit

• We could do this in sequence, but lets do it in parallel as the extraction of the rainfall for each slope unit are independent tasks

Use the results from the new slide_dates_in_polygon-function

• The function get_rainfall_for_polygons iterates over the dataframe returned by slides_dates_in_polygon

• This is done in parallel as each row is one spatial object and is idependent of the others

• In the end they are just stacked up to build one big dataframe

# filte out the slides before the start of the rainfall data
su_points_in_poly = su_points_in_poly[su_points_in_poly$date > as.Date("1980-01-01"), ]

# get the antecedent rainfall for all the dates in all the slope units with landslides
slopeunits_rainfall = rainfallR::get_rainfall_for_polygons(su_points_in_poly, nc_var = "precipitation", days_back = 14, data_path = path_ncdf)

What have we got back?

Within the get_rainfall_for_polygons-function:

• The function ex_rainfall which is called as many times as rows in the su_points_in_poly-object returns a dataframe

• We used the .rbind-parameter for the .combine.argument, slopeunits_rainfall is one large dataframe

# verfiy that each slope unit has 6 entries for all the six days before the slide
slopeunits_rainfall %>% 
  group_by(poly_id, PIFF_ID) %>% 
  summarise(n = n()) %>% 
  ungroup() %>% 
  count(n)

look at the antecedent rainfall

g =slopeunits_rainfall %>% 
  group_by(poly_id, PIFF_ID) %>% 
  mutate(unique_id = dplyr::cur_group_id()) %>% 
  ggplot() +
  geom_path(mapping = aes(x = days_before_event,
                  y = cumsum,
                  group = unique_id,
                  color = second_level)) +
  scale_x_reverse() +
  scale_color_manual(values = c("translational slide"="#D4f0f0", "rotational slide" = "#ffc8a2")) +
  theme_light() +
  labs(
    title = "Five day antecedent rainfall", 
   y = "cumulative sum [mm]",
   x = "Days before landslide"

  ) +
theme(
 panel.background = element_rect(fill="black"),
 plot.background = element_rect(fill="black"),
 panel.grid = element_blank(),
 plot.margin = margin(20, 20, 20, 20),
 axis.text = element_text(colour = "white")
)

g = ggplotGrob(g)
bg = g$grobs[[1]] 
round_bg = roundrectGrob(x=bg$x, y=bg$y, width=bg$width, height=bg$height,
                          r=unit(0.1, "snpc"),
                          just=bg$just, name=bg$name, gp=bg$gp, vp=bg$vp)
g$grobs[[1]] <- round_bg
plot(g)

table(slopeunits_rainfall$days_before_event)

RobinKohrs/rainfallR documentation built on Oct. 3, 2021, 1:42 a.m.