index.Rmd

---
title: "DenmarkSF"
author: "gntem2"
date: "02/03/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

##
This section deals with downloading NUTS region from eurostat.

Next the geolocation of hospitals is extracted using tmapstools

```{r europe}
library(dplyr)
library(tidyverse)
library(sf)
library(eurostat)
library(leaflet)
library(mapview)
library(tmaptools)
#download geospatial data for NUTS-3 regions

#euro_nuts3_sf <-
# eurostat::get_eurostat_geospatial(output_class = #'sf', resolution = '60', nuts_level = 3) %>%
# sf::st_transform(crs = 3035)
#save(euro_nuts3_sf,file="euro_nuts3_sf.Rda")
#load("euro_nuts3_sf.Rda")

#download geospatial data for NUTS-2 regions
#euro_nuts2_sf <-
# eurostat::get_eurostat_geospatial(output_class = 'sf', resolution = '60', nuts_level = 2) %>%
# sf::st_transform(crs = 3035)
#save(euro_nuts2_sf,file="euro_nuts2_sf.Rda")
load("euro_nuts2_sf.Rda")

#Demark = DK
#^ denotes start with
#$denotes end with
#DKnuts3_sf<- euro_nuts3_sf%>% filter(str_detect(NUTS_ID,"^DK"))

DKnuts2_sf<- euro_nuts2_sf%>% filter(str_detect(NUTS_ID,"^DK"))

##https://ec.europa.eu/eurostat/web/nuts/background
##NUTS 1: major socio-economic regions
##NUTS 2: basic regions for the application of regional policies
##NUTS 3: small regions for specific diagnoses

load("europeLBDK.Rda")
load("europeRGDK.Rda")

#geocode hospital addresses using tmaptools free
#does not need to get key like google maps api

#hosp_addresses <- c(AarhusHospital = "aarhus university hospital, aarhus, Denmark",
#                     AalborgHospital = "aalborg university hospital, aalborg, Denmark",
 #                   HolstebroHospital = "Regionshospitalet Holstebro, Holstebro, Denmark",
#                    VejleHospital="Vejle Sygehus,Beriderbakken 4, Vejle, Denmark",
#                    EsbjergHospital="Esbjerg Sygehus, Esbjerg, Denmark",
                    #SoenderborgHospital="Soenderborg Sygehus,1 Sydvang, Sønderborg Denmark",
#        OdenseHospital="Odense Sygehus, Odense, Denmark",   
#        RoskildeHospital="Roskilde Sygehus,  Roskilde, Denmark",  
#        BlegdamsvejHospital="Rigshospitalet blegdamsvej, 9 Blegdamsvej, København, Denmark",  
#        GlostrupHospital="Rigshospitalet Glostrup, Glostrup, Denmark")

#HospLocations <- tmaptools::geocode_OSM(hosp_addresses, as.sf=TRUE)

#HospLocations <- sf::st_transform(HospLocations,
 #            sf::st_crs(europeRGDK))

#HospLocations$Center<-c("CSC", "PSC", "PSC", "PSC", "PSC", "PSC", "CSC", "PSC", "CSC", "PSC")

#save(HospLocations,file="HospLocations.Rda")
load("HospLocations.Rda")

#helicopter df
load("DenmarkHelicopter.Rda")

dfhelicopter<-st_as_sf(df,
    coords = c("Lon","Latitude"),
      crs=4326)
  
```

## different types of plotting methods
```{r plots}

#######################
#plot using simple features by region names
plot(DKnuts2_sf["NUTS_NAME"])


#######################
#ggplot with simple features
#colour by NUTS
ggplot(data=DKnuts2_sf,aes(fill=NUTS_NAME))+geom_sf()

g<-ggplot() + # set up the framework
  geom_sf(data=DKnuts2_sf,aes(fill=NUTS_NAME))+
  #geom_point(data=df_sf,aes(x=Lon, y=Latitude), fill="orange", pch=21, alpha=0.7, size=2)+
  labs(x="Longitude (WGS84)", y="Latitude", title="Map of Points") + 
  theme_bw() # change this to sans if it doesn't plot
g



```
```{r euope}
#https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/communes
#shapefile for europe
#use sf
#europeBN<-st_read("./COMM-01M-2013-SH/COMM_01M_2013_SH/data/COMM_BN_01M_2013.shp")


#########
#europeLB<-st_read("./COMM-01M-2013-SH/COMM_01M_2013_SH/data/COMM_LB_2013.shp")

#filter starts with DK
#europeLBDK<-europeLB %>% filter(str_detect(COMM_ID,"^DK"))

#save(europeLBDK,file="europeLBDK.Rda")
#load("europeLBDK.Rda")
#######

#europeRG<-st_read("./COMM-01M-2013-SH/COMM_01M_2013_SH/data/COMM_RG_01M_2013.shp")


#europeRGDK<-europeRG %>% filter(str_detect(COMM_ID,"^DK"))
#save(europeRGDK,file="europeRGDK.Rda")
#load("europeRGDK.Rda")
mapview(europeRGDK["Shape_Area"])


```
##convert pdf file to xcel file using pdftables
##extract data from xcel file
##errors occurred in conversion with Danish names 
```{r extraction}
dk<-read.csv("denmarkstrokepdf.csv")

#extract only data on large regions=NUTS2
dk2<-dk[c(4:8),]

#clean up column X.1 containing stroke data
#remove numerator before back slash
#https://github.com/STAT545-UBC/Discussion/issues/394

#remove number before slash sign
#then remove slash sign
dk2$strokenum<-str_replace(dk2$X.1,"[0-9]*","") %>%
  str_replace("/","\\")

dk2$Uoplyst<-str_replace(dk2$Uoplyst,"Sjælland","Sjælland")
```

```{r combine}
#merge sf file for DK nuts2 with stroke number
DKnuts2_sf2<-right_join(DKnuts2_sf,dk2,by=c("NUTS_NAME"="Uoplyst"))

#add hospital location data to mapview using+
#label different types of hospitals
#n<-mapview(DKnuts2_sf2["strokenum"])+mapview(HospLocations, zcol="Center")

#include helicopter data
n<-mapview(DKnuts2_sf2["strokenum"])+mapview(dfhelicopter, zcol="Center")

#make pics using mapshot
mapshot(n, url = paste0(getwd(),file="/denmark_stroke_nuts2.html"), file = paste0(getwd(), "/denmark_stroke_nuts2.png"))

n



```

##estimate distance from each commune to hospital
```{r distance}
#determine distance hospital to centroid
dist_to_loc <- function (geometry, location){
    units::set_units(st_distance(st_centroid (geometry), location)[,1], km)
}

#set distance 10 km
#change to 30 km
dist_range <- units::set_units(30, km)


##
europeRGDK <- mutate(europeRGDK,
       DirectDistanceToAarhus = dist_to_loc(geometry,HospLocations["AarhusHospital", ]),
       DirectDistanceToAalborg     = dist_to_loc(geometry,HospLocations["AalborgHospital", ]),
       DirectDistanceToHolstebro     = dist_to_loc(geometry,HospLocations["HolstebroHospital", ]),
       DirectDistanceToVejle     = dist_to_loc(geometry,HospLocations["VejleHospital", ]),
       DirectDistanceToEsbjerg     = dist_to_loc(geometry,HospLocations["EsbjergHospital", ]),
       DirectDistanceToSoenderborg = dist_to_loc(geometry,HospLocations["SoenderborgHospital", ]),
       DirectDistanceToOdense     = dist_to_loc(geometry,HospLocations["OdenseHospital", ]),
       DirectDistanceToRoskilde     = dist_to_loc(geometry,HospLocations["RoskildeHospital", ]),
       DirectDistanceToBlegdamsvej     = dist_to_loc(geometry,HospLocations["BlegdamsvejHospital", ]),
       DirectDistanceToGlostrup     = dist_to_loc(geometry,HospLocations["GlostrupHospital", ]),
       #
       DirectDistanceToNearest   = pmin(DirectDistanceToAarhus,
                                      DirectDistanceToAalborg,DirectDistanceToHolstebro,
                                 DirectDistanceToVejle,DirectDistanceToEsbjerg, DirectDistanceToSoenderborg,DirectDistanceToOdense,DirectDistanceToRoskilde,DirectDistanceToBlegdamsvej,DirectDistanceToGlostrup
                                 
                                      )
    )

#reduce from 2178 communes to 992
StrokeHosp <- filter(europeRGDK,
                                 DirectDistanceToNearest < dist_range) %>%
        mutate(Postcode = as.numeric(COMM_ID)) %>%
        select(-starts_with("DK"))


p<-mapview(StrokeHosp,map.type="OpenStreetMap.HOT", alpha.regions=0.5)+mapview(HospLocations, zcol="Center")+mapview(DKnuts2_sf2["strokenum"])
p
mapshot(p, url = paste0(getwd(),file="denmark_stroke_nuts2_catchment.html"), file = paste0(getwd(), "/denmark_stroke_nuts2_catchment.png"))

```


Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
```{r helicopter} 
if (file.exists("HeliLocations.Rda")) {
  load("HeliLocations.Rda")
} else {
#geocode helicopter addresses using tmaptools free
#does not need to get key like google maps api
heli_addresses <- c(AalborgHeli = "ny lufthavnsvej 100, 9400 Nørresundby, Denmark", 
                    SkiveHeli = "Killesmosevej 1, 7800 Skive, Denmark", 
                    BillundHeli = "Firhøjevej 20, 7190 Billund, Denmark", 
                    RingstedHeli = "Haslevvej 56, 4100 Ringsted, Denmark")

HeliLocations <- tmaptools::geocode_OSM(heli_addresses, as.sf=TRUE)

HeliLocations <- sf::st_transform(HeliLocations,
           sf::st_crs(europeRGDK))


save(HeliLocations,file="HeliLocations.Rda")
}

#determine distance helicopter to centroid
dist_to_loc_Heli <- function (geometry, location){
    units::set_units(st_distance(st_centroid (geometry), location)[,1], km)
}
#set distance 10 km
#change to 200 km
dist_range <- units::set_units(200, km)


##
europeRGDK <- mutate(europeRGDK,
       DirectDistanceToAalborgHeli = dist_to_loc(geometry,HeliLocations["AalborgHeli", ]),
       DirectDistanceToSkiveHeli     = dist_to_loc(geometry,HeliLocations["SkiveHeli", ]),
       DirectDistanceToBillundHeli     = dist_to_loc(geometry,HeliLocations["BillundHeli", ]),
       DirectDistanceToRingstedHeli     = dist_to_loc(geometry,HeliLocations["RingstedHeli", ]),
              #
       DirectDistanceToNearestHeli   = pmin(DirectDistanceToAalborgHeli,           DirectDistanceToSkiveHeli, DirectDistanceToBillundHeli,                           DirectDistanceToRingstedHeli))

#reduce from 2178 communes to 992
StrokeHeli <- filter(europeRGDK,
                                 DirectDistanceToNearestHeli < dist_range) %>%
        mutate(Postcode = as.numeric(COMM_ID)) %>%
        select(-starts_with("DK"))


p<-mapview(StrokeHeli,map.type="OpenStreetMap.HOT", alpha.regions=0.5)+mapview(HeliLocations)+mapview(DKnuts2_sf2["strokenum"])
p
mapshot(p, url = paste0(getwd(),file="/denmark_stroke_nuts2_catchment.html"), file = paste0(getwd(), "/denmark_stroke_nuts2_catchment.png"))
```