# INBO CODING CLUB
28 March 2020

Welcome!

## Share your code snippet

If you want to share your code snippet, copy paste your snippet within a section of three backticks (```):

As an **example**:

library(tidyverse)

(*you can copy paste this example and add your code further down*)

### Deelnemers

Name | Challenges
--- | ---
Joost Vanoverbeke | * *
Dirk Maes | *
Laurian Van Maldeghem | * * *
Patrik Oosterlynck | * * *
Ruben Perez Perez | * * *
Marc Portier | * * o
Els Lommelen | * *
Arno Thomaes | * * *
Leen Govaere | * * *
Jim Casaer |
Loïc van Doorn |
Wim Mertens | *
Anneleen Rutten | 
Raisa Carmen | * * * *
An Vanden Broeck | * *
Luc De Bruyn |
Gizem Poffyn | ***
Els De Bie |
Peter Desmet | * *
Marijke Thoonen |
Damiano Oldoni |
Salvador Fernandez | * * *
Alexander Van Braeckel |**
Anton Van de Putte |
Hans Van Calster | 
Pieterjan Verhelst | * *

No yellow sticky notes online. Just put a "|" if you are present, a space and add a "*" each time you solve a challenge.
_E.g. Damiano Oldoni is present and has solved the 1st challenge):
Damiano Oldoni | *_

> suggestion:
> * change | into - if you leave
> * change * into o if you skip a challenge?
> * eg Marc Portier | * o * *

### Screen Tip

Maybe you could divide the screen in two to tip in R while following Damiano's live coding at the same time.

### Another tip

Patrik: i find the mapview package very handy (read easy ;-)) for a quick and dirty view of geospatial stuff. Check it out: https://www.youtube.com/watch?v=p2mu15exMSY

### Challenge 1

map<-tm_shape(provinces) +
  tm_borders(col = 'black', lwd = 0.3) +
  tm_fill(col = 'SURFACE.GIS.km2', title = 'Area (km2)') +
  tm_facets(by="TX_RGN_DESCR_NL",nrow=2)
tmap_mode('view')
map

### add text with province names to the map (Dirk)

map <- tm_shape(provinces) +
	tm_borders(col = 'black',
						 lwd = 1) +
	tm_fill(col = 'SURFACE.GIS.km2',
					title = 'Area (km²)') +
	tm_facets(by="TX_RGN_DESCR_NL",
						nrow = 2) +
	tm_text("TX_PROV_DESCR_NL") 
tmap_mode('view')
map

With `tm_text()` you can add labels to the shapes, e.g. the name of the provinces.

### Wim - another way

tmap_mode('view') 
a <- provinces %>% filter(TX_RGN_DESCR_NL == "Vlaams Gewest") %>% tm_shape() +
  tm_fill(col = 'SURFACE.CAD.km2', title = 'Oppervlakte (km²)') + 
  tm_borders() +
  tm_layout(title = "Vlaams gewest",
    legend.outside = TRUE, legend.outside.position="right")
b <- provinces %>% filter(TX_RGN_DESCR_NL == "Waals Gewest") %>% tm_shape() +
  tm_fill(col = 'SURFACE.CAD.km2', title = 'Oppervlakte (km²)') + 
  tm_borders() +
  tm_layout(title = "Waals gewest",
    legend.outside = TRUE, legend.outside.position="right")

tmap_arrange(a, b)

### Intermezzo: ggplot and geospatial data

Did you know that you could use ggplot to plot geospatial data as well?

ggplot will transform your data in 4326 (WGS 84 projection) but will maintain the original orientation. Check the code here below

# Plot provinces in original projection
ggplot(provinces) +
  geom_sf()

# Provinces in CRS: epsg code 3130
provinces_3130 <- st_transform(provinces, 3130)
st_crs(provinces_3130)

ggplot(provinces_3130) +
  geom_sf()

### Intermezzo: tmap view mode

You can switch between static and interactive mode. Two options:

1. `ttm()`: it switches between the modes.
2. `tmap_mode('plot')` or `tmap_mode('view')`

### sf data.frame and leaflet

A nice property of leaflet is that it recognizes authomatically the latitude/longitude from the column `geometry`. Of course, you have to provide data in WGS84. If don't, you get this warning:
> Warning messages:
1: sf layer is not long-lat data
2: sf layer has inconsistent datum (+proj=tmerc +lat_0=0 +lon_0=23 +k=1 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ).
Need '+proj=longlat +datum=WGS84'

And your points are put in totally worng positions.

### Challenge 2

Raïsa:

m <- leaflet(data=obs_butterflies) %>% 
     setView(lng = 4.89, lat = 51.01, zoom = 7)%>%
     addTiles() %>% 
     addMarkers(~decimal_longitude, ~decimal_latitude,clusterOptions = markerClusterOptions(),
                 label = ~as.character(species), popup= ~sprintf('Date: %s',obs_butterflies$date))

Peter:

map <- leaflet(data = obs_butterflies) %>%
  setView(lng = 4.89, lat = 51.01, zoom = 7) %>%
  addProviderTiles(providers$Stamen.Toner) %>%
  addMarkers(
    ~decimal_longitude,
    ~decimal_latitude,
    popup = ~paste("Date:", as.character(date)),
    label = ~species,
    clusterOptions = markerClusterOptions()
)
map

leaflet(obs_butterflies) %>%
  setView(lng = 4.89, lat = 51.01, zoom = 7) %>%
  addTiles() %>%
  addMarkers(clusterOptions = markerClusterOptions(),
             popup = ~htmlEscape(paste0("Date: ", paste(year,
                                                       month,
                                                       day,
                                                       sep = "-"))),
             label = ~htmlEscape(species))

### Challenge 3
Raïsa:

#Create a color palette
pal <- colorNumeric(
  palette = "YlGnBu",
  domain = provinces$SURFACE.GIS.km2
)

#create the map of belgium
m <- leaflet(provinces) %>%
  addTiles() %>%  # Add default OpenStreetMap map 
  addPolygons(weight=2, 
            color='black', 
            fillOpacity = 0.5,
            fillColor = ~pal(SURFACE.GIS.km2)) %>%
  addMarkers(~lng, ~lat, label = ~as.character(TX_PROV_DESCR_NL)) %>% 
  addLegend("bottomright", pal = pal, values = ~SURFACE.GIS.km2,
            title = "Surface area (km2)",
            opacity = 1)

Patrik:

bins <- c(1000,1500 , 2000, 2500, 3000, 3500, 4000, 4500, Inf)
pal <- colorBin("viridis", domain = provinces$SURFACE.GIS.km2, bins = bins)


leaflet(provinces) %>%
  addTiles() %>%
  addPolygons(fillColor = ~pal(SURFACE.GIS.km2), color = "white", fillOpacity = 1) %>%
  addMarkers(~lng, ~lat, label = ~TX_PROV_DESCR_NL)
  addLegend(pal = pal,
            values=~SURFACE.GIS.km2,
            opacity = 0.7, title = "Area (km2)",
            position = "bottomright")

Hans:

leaflet(provinces) %>%
  addTiles() %>%
  addPolygons(fillColor = ~pal(SURFACE.GIS.km2),
              label = ~htmlEscape(TX_PROV_DESCR_NL),
              popup = ~htmlEscape(paste("Area:", round(SURFACE.GIS.km2)))) %>%
  addLegend(pal = pal,
            values = ~SURFACE.GIS.km2,
            title = "Area (km<sup>2</sup>)")

###Bonus challenge

Raïsa:

#create a color palette for the butterfly species and the provinces
pal1 <- colorFactor(c("yellow", 'green','black',"red"),
        domain = c("Citroenvlinder","Grote vos","Kleine vos","Oranjetipje"))
pal2 <- colorNumeric(
  palette = "YlGnBu",
  domain = provinces$SURFACE.GIS.km2
)
# make the map
m <- leaflet(obs_butterflies_selection) %>%
  addTiles(group = "OSM (default)") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "Esri") %>%
  addPolygons(data=provinces, weight=2, color='black', fillOpacity = 0.5,
              fillColor = ~pal2(SURFACE.GIS.km2),group = "Province") %>%
  addCircles(~decimal_longitude, ~decimal_latitude, radius= ~coord_uncertainty_in_meters , label = ~species, color = ~pal1(species),group = "Butterflies") %>%
  # Layers control
  addLayersControl(
    baseGroups = c("OSM (default)", "Esri"),
    overlayGroups = c("Butterflies", "Province"),
    options = layersControlOptions(collapsed = FALSE)
  )%>% 
  addLegend("bottomright", pal = pal1, values = ~species,
            title = "Butterfly species",
            #labFormat = labelFormat(prefix = "$"),
            opacity = 1)
m