r_mapping-notes.R

# # # RESOURCES -----------------------------------------------------------------------------------------------------------------
#
# https://media.readthedocs.org/pdf/a-little-book-of-r-for-time-series/latest/a-little-book-of-r-for-time-series.pdf
# http://brettromero.com/wordpress/data-science-a-kaggle-walkthrough-introduction/
# http://amsantac.co/blog.html
#
# - boundaries
#   * https://stat.ethz.ch/pipermail/r-sig-geo/2011-May/011814.html
#   * http://eco-data-science.github.io/spatial_analysis2_R/
#   * https://mgimond.github.io/Spatial/index.html
#   * https://rpubs.com/BeccaStubbs/bringing_shapefiles_into_R
#   * https://github.com/Robinlovelace/sdvwR
#   *
#
### - mapping general
#   * http://eriqande.github.io/rep-res-web/lectures/making-maps-with-R.html
#   * http://zevross.com/blog/2014/07/16/mapping-in-r-using-the-ggplot2-package/
#   * https://cran.r-project.org/web/packages/maptools/vignettes/combine_maptools.pdf
#   * https://stat.ethz.ch/pipermail/r-sig-geo/2011-May/011814.html
#   * http://www.molecularecologist.com/2012/09/making-maps-with-r/
#   * http://www.geog.uoregon.edu/GeogR/examples/maps_examples01.htm
#   * http://www.kevjohnson.org/making-maps-in-r/
#   * http://www.computerworld.com/article/3038270/data-analytics/create-maps-in-r-in-10-fairly-easy-steps.html
#   * http://clarkdatalabs.github.io/mapping_R/
#   * http://geog.uoregon.edu/bartlein/courses/geog495/lec06.html
#   * http://stackoverflow.com/documentation/r/1372/introduction-to-geographical-maps#t=201702271434306481539
#   * http://www.molecularecologist.com/2016/07/making-maps-in-r-volume-2-ggplots/
#   * https://r-resources.massey.ac.nz/lurn/LURNch16.html
#   * http://axismaps.github.io/thematic-cartography/
#   * https://cran.r-project.org/doc/contrib/intro-spatial-rl.pdf
#   * https://cran.r-project.org/web/packages/maptools/vignettes/combine_maptools.pdf
#   * https://github.com/Robinlovelace/Creating-maps-in-R/blob/master/intro-spatial.Rmd
#   * http://urbanspatialanalysis.com/portfolio/predicting-gentrification-using-longitudinal-census-data/
#   * http://axismaps.github.io/thematic-cartography/articles/classification.html
#   *
#   *
#
# - ggplot (for mapping)
#   * http://tutorials.iq.harvard.edu/R/Rgraphics/Rgraphics.html
#   * http://www.gradaanwr.net/
#   * http://stackoverflow.com/questions/19791210/r-ggplot2-merge-with-shapefile-and-csv-data-to-fill-polygons
#   * http://rmhogervorst.nl/cleancode/blog/2017/01/06/plotting-a-map-with-ggplot2.html
#   * http://r-statistics.co/Top50-Ggplot2-Visualizations-MasterList-R-Code.html
#   *
#
# - ggmap
#   * https://github.com/dkahle/ggmap
#   * https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf
#   * https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/ggmap/ggmapCheatsheet.pdf
#   * https://blog.dominodatalab.com/geographic-visualization-with-rs-ggmaps/
#   * http://www.kevjohnson.org/making-maps-in-r-part-2/
#   * http://bcb.dfci.harvard.edu/~aedin/courses/R/CDC/maps.html
#   * http://urbanspatialanalysis.com/dataviz-tutorial-mapping-san-francisco-home-prices-using-r/
#   *
#   *
#
# - tmap
#   * https://github.com/Arevaju/mapping-tool
#   * https://behrica.github.io/blog/2016/10/food-consumption-tutorial.html
#   * http://www.computerworld.com/article/3175623/data-analytics/mapping-in-r-just-got-a-whole-lot-easier.html
#   *
#
# - leaflet
#   * http://zevross.com/blog/2015/10/14/manipulating-and-mapping-us-census-data-in-r-using-the-acs-tigris-and-leaflet-packages-3/#convert-the-data.frame-back-to-a-spatialpolygonsdataframe
#   * https://lilbigdataboy.wordpress.com/2017/01/14/les-crimes-de-san-francisco-en-temps-reel-shiny-et-leaflet/
#   *
#   *
#

# - projects
#   * http://urbanspatialanalysis.com/portfolio/predicting-gentrification-using-longitudinal-census-data/
#   * http://blog.revolutionanalytics.com/2017/02/forecasting-gentrification.html
#   * http://urbanspatialanalysis.com/dataviz-tutorial-mapping-san-francisco-home-prices-using-r/
#   * http://zevross.com/blog/2016/04/19/r-powered-web-applications-with-shiny-a-tutorial-and-cheat-sheet-with-40-example-apps/
#   * http://www.computerworld.com/article/2497464/business-intelligence/business-intelligence-60-r-resources-to-improve-your-data-skills.html
#   *
#

# --------------------------------------------------------------------------------
mytext <- "
    The following is the R code to create arbitrary aggregation boundaries, both in shapefile (for tmap and leaflet) and dataframe (for ggplot and ggmap) format, starting from some base shapefile layer.
    Specifically, in the UK the reference is the union of Output Areas in England and Wales (EW), Output Areas in Scotland (SC), and Small Areas in Northern Ireland (NI).
    There is no official global UK map for these census 2011 small areas, so it has to be created binding together the three different maps, each available from the corresponding Statistical Office as for the following links:
    EW: Output Areas (Generalised Clipped) http://geoportal.statistics.gov.uk/datasets/09b8a48426e3482ebbc0b0c49985c0fb_2
    SC: Output Areas https://www.nrscotland.gov.uk/statistics-and-data/geography/our-products/census-datasets/2011-census/2011-boundaries
    NI: Small Areas http://www.nisra.gov.uk/geography/SmallAreas.htm
    Once downloaded, rename each group of files with the corresponding two characters acronym. It's necessary to keep only the files with the four following extensions: shp, shx, prj, dbf.
    One small problem to be aware of is the different projection of the boundaries:
    EW: '+init=epsg:27700'  # http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
    SC: '+init=epsg:27700'   # http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
    NI: '+init=epsg:29902'   # http://spatialreference.org/ref/epsg/29902/
    but all of the above will be transformed in the standard WGS84 for web mapping.
    The following are the number of features in each set of boundaries:
    EW: 181,408 (E 171,372 + W 10,036)
    SC:    46,351
    NI:       4,537
    for a total UK of 232,296. This is important to know because when simplifying the polygons overall, some of the features could dissolve.
    To successfully aggregate, the script needs also a lookup table that gives back a unique aggregated area code for every output/small area (from now on just OA).
    This table could be built for many types of aggregations using, for example, the Postcodes Directories
    NHS Postcode Directory: http://ons.maps.arcgis.com/home/item.html?id=4e63d9b589604c0ba10fbbd3bef08e6f
    ONS Postcode Directory: http://ons.maps.arcgis.com/home/item.html?id=dfa0ff74981b4b228d2030d852f0b14a
    These are maintained and published every quarter by the ONS, that relates both current and terminated postcodes in the UK to a range of current and past statutory administrative, electoral, health and other geographies. The downside to this approach is that some OA could belong to different aggregations, and it becomes somewhat arbitrary, without more knowledge, decide which is the right destination area.
    When saving to the database, there is an added column 'type' that is part of the foreign key ('id', 'type') connecting to the primary key ('boundary_id', 'type') in the table *boundaries_ids*, used to retrieve the correct id of the feature to link to external data.
    I saved all files in s:\ICS\Projects\Data Management\R_Shiny\Resources\boundaries_aggregations\
    It is finally important to understand that when reproducing ONS content, the output must include a source accreditation to ONS.
"


# # # STEP 0a - PREPARATION OF BOUNDARIES BY COUNTRIES  -------------------------------------------------------------------

# Boundaries from Census 2011
# EW: Output Area (Generalised Clipped) http://geoportal.statistics.gov.uk/datasets/09b8a48426e3482ebbc0b0c49985c0fb_2
# SC: Output Area https://www.nrscotland.gov.uk/statistics-and-data/geography/our-products/census-datasets/2011-census/2011-boundaries
# NI: Small Area http://www.nisra.gov.uk/geography/SmallAreas.htm
# IE: Small Area http://www.cso.ie/en/census/census2011boundaryfiles/  (=> they lack the projection file)
# Once downloaded, rename each group of files with the corresponding two characters acronym

# Original projections
# ewgrid = '+init=epsg:27700'  # http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
# scgrid = '+init=epsg:27700'  # http://spatialreference.org/ref/epsg/osgb-1936-british-national-grid/
# nigrid = '+init=epsg:29902'  # http://spatialreference.org/ref/epsg/29902/
# iegrid = '+init=epsg:29902'  # http://spatialreference.org/ref/epsg/29902/

# Number of features:
#  - EW: 181408 (171372+10036)
#  - SC:  46351
#  - NI:   4537
#  UK => 232296

# Reference Country codes
# E92000001 = England;
# W92000004 = Wales;
# S92000003 = Scotland;
# N92000002 = Northern Ireland;

# # # STEP 0b - LOAD PACKAGES ---------------------------------------------------------------------------------------------------
pkg <- c('broom', 'classInt', 'data.table', 'ggplot2', 'ggmap', 'maptools', 'RColorBrewer', 'rgeos', 'rgdal', 'rmapshaper', 'RMySQL', 'sf', 'tmap')
invisible( lapply(pkg, require, character.only = TRUE) )

# Additional linux libraries needed to install geospatial packages
# - CentOS 7:
#   * sudo yum install epel-release
#   * sudo yum install gdal
#   * sudo yum install gdal-devel
#   * sudo yum install geos-devel
#   * sudo yum install proj-epsg
# - Ubuntu 14.04 LTS:
#   * sudo aptitude install libproj-dev
#   * sudo aptitude install libgdal-dev
#   * sudo apt-get install libv8-dev


# # # STEP 0c - SET VARIABLES ---------------------------------------------------------------------------------------------------

# set the directory of the boundaries shapefiles
boundaries.path <- 'C:/cloud/OneDrive - University College London/data/boundaries/shp'
# set the projection string for WGS84
proj.wgs <- '+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0'
# set UK centroid and bounded box
UK_centre <- c(lon = -2.421976, lat = 53.825564)
UK_bounds <- c(lng1 = 1.8, lat1 = 49.9, lng2 = -8.3, lat2 = 59.0 )


# # # STEP 1 - LOAD AND MERGE EW, SC AND NI OAs MAP INTO A UNIQUE UK MAP  -------------------------------------------------------

# England and Wales
# read polygons from shapefile (using rgdal, projection is automatically applied if the prj file is found)
shp.ew <- readOGR(boundaries.path, layer = 'EW')
# check the projection, and read the field to keep as future id
summary(shp.ew)
# transform the shapefile to WGS84 projection
shp.ew <- spTransform(shp.ew, CRS(proj.wgs))
# keep in the data slot only the ONS Output Area id, renaming it as 'id'
shp.ew <- shp.ew[, 'oa11cd']
colnames(shp.ew@data) <- c('id')
# reassign the polygon IDs
shp.ew <- spChFIDs(shp.ew, as.character(shp.ew$id))

# Scotland
shp.sc <- readOGR(boundaries.path, layer = 'SC')
summary(shp.sc)
shp.sc <- spTransform(shp.sc, CRS(proj.wgs))
shp.sc <- shp.sc[, 'code']
colnames(shp.sc@data) <- c('id')
shp.sc <- spChFIDs(shp.sc, as.character(shp.sc$id))

# N.Ireland
shp.ni <- readOGR(boundaries.path, layer = 'NI')
summary(shp.ni)
shp.ni <- spTransform(shp.ni, CRS(proj.wgs))
shp.ni <- shp.ni[, 'SA2011']
colnames(shp.ni@data) <- c('id')
shp.ni <- spChFIDs(shp.ni, as.character(shp.ni$id))

# UK as a merge of all previous boundaries
shp.uk <- spRbind(spRbind(shp.ew, shp.sc), shp.ni)
# reduce the complexity of the boundaries
# shp.uk <- ms_simplify(shp.uk, keep = 0.05)
# save merged polygons as a unique shapefile
try( do.call(file.remove, as.list(paste0(boundaries.path, '/', list.files(path = boundaries.path, pattern = 'UK') ) ) ), silent = FALSE )
writeOGR(shp.uk, dsn = boundaries.path, layer = 'UK', driver = 'ESRI Shapefile')


# # # STEP 2 - CREATE AGGREGATION MAPS USING UK MAP AND LOOKUP TABLE ------------------------------------------------------------

# Define the area type that will be the base for aggrega
base.loc <- 'OA'
# read polygons from shapefile and apply projection
shp.uk <- readOGR(boundaries.path, layer = base.loc)
# connect to database
db_conn <- dbConnect(MySQL(), user = 'root', password = 'root', dbname = 'geography')
# retrieve lookup table
lookups <- dbReadTable(db_conn, 'lookups')
# close db connection
dbDisconnect(db_conn)
# list of aggregations to process
# areas <- c('LSOA', 'MSOA', 'LAD', 'CTY', 'RGN', 'CTRY')
areas <- c('CCG', 'LAT', 'NHSR', 'CCR', 'CTRY')
for(area in areas){
    # join shapefile data slot and lookup table on the area code
    shp.uk.tmp <- merge(shp.uk, lookups[, c(base.loc, area)], by.x = 'id', by.y = base.loc)
    # Build the list of subareas
    subareas <- sort(unique(shp.uk.tmp[[area]]))
    # Define first area
    subarea <- subareas[1]
    # Print a processing message
    print(paste('Processing', area, 'subarea', subarea, '- number 1 out of', length(subareas)))
    # select all OA contained in first subarea
    shp.area <- subset(shp.uk.tmp, shp.uk.tmp[[area]] == subarea)
    # dissolve submap
    shp.area <- ms_dissolve(shp.area)
    # define object id
    shp.area$id <- subarea
    shp.area <- spChFIDs(shp.area, as.character(shp.area$id))
    # proceed in the same way for all other subareas, attaching every results to previous object
    for(idx in 2:length(subareas)){
        subarea <- subareas[idx]
        print(paste('Processing', area, 'subarea', subarea, '- number', idx, 'out of', length(subareas)))
        shp.tmp <- subset(shp.uk.tmp, shp.uk.tmp[[area]] == subarea)
        shp.tmp <- ms_dissolve(shp.tmp)
        shp.tmp$id <- subarea
        shp.tmp <- spChFIDs(shp.tmp, as.character(shp.tmp$id))
        shp.area <- spRbind(shp.area, shp.tmp)
    }
    print(paste0('Simplifying ', area, 's'))
    # reduce the details of the boundaries for less space (it's a statistical map, not an OS Explorer Map!)
    shp.area <- ms_simplify(shp.area, keep = 0.05)
    # Print a processing message
    print(paste('Saving', area, 'in dataframe format'))
    # save lookup (rmapshaperid, id) adding the area type
    lkp <- shp.area@data
    lkp$type <- area
    colnames(lkp) <- c('boundary_id', 'id', 'type')
    # connect to database
    db_conn <- dbConnect(MySQL(), user = 'root', password = 'root', dbname = 'geography')
    # save lookup to database
    dbSendQuery(db_conn, paste0("DELETE FROM boundaries_ids WHERE type = '", area, "'") )
    dbWriteTable(db_conn, 'boundaries_ids', lkp, row.names = FALSE, append = TRUE)
    # close db connection
    dbDisconnect(db_conn)
    # create dataframe suitable for ggplot
    df.area <- tidy(shp.area)
    # change names in dataframe to coordinates to avoid possible mismatching with programming languages keywords
    setnames(df.area, c('long', 'lat'), c('X_lon', 'Y_lat'))
    # add "type" column
    df.area$type <- area
    # connect to database
    db_conn <- dbConnect(MySQL(), user = 'root', password = 'root', dbname = 'geography')
    # save dataframe to database
    dbSendQuery(db_conn, paste0("DELETE FROM boundaries WHERE type = '", area, "'") )
    dbWriteTable(db_conn, 'boundaries', df.area, row.names = FALSE, append = TRUE)
    # close db connection
    dbDisconnect(db_conn)
    # Print a processing message
    print(paste('Saving', area, 'in shapefile format'))
    # delete the rmapshaperid from Polygons
    shp.area <- shp.area[, 'id']
    # save Polygons as shapefile
    try( do.call(file.remove, as.list(paste0(boundaries.path, '/', list.files(path = boundaries.path, pattern = area) ) ) ), silent = TRUE )
    writeOGR(shp.area, dsn = boundaries.path, layer = area, driver = 'ESRI Shapefile')
    # some polishing...
    rm(df.area)
    rm(shp.area)
    gc()
}


# # # STEP 2b - ADDING IRELAND ---------------------------------------------------------------------------------------------------------

shp.ie <- readOGR(boundaries.path, layer = 'IE_OA')
shp.ie <- spTransform(shp.ie, CRS(proj.wgs))
summary(shp.ie)
shp.ie <- shp.ie[, 'SMALL_AREA']
colnames(shp.ie@data) <- c('id')
shp.ie <- spChFIDs(shp.ie, as.character(shp.ie$id))
try( do.call(file.remove, as.list(paste0(boundaries.path, '/', list.files(path = boundaries.path, pattern = 'IE_OA') ) ) ), silent = TRUE )
writeOGR(shp.ie, dsn = boundaries.path, layer = 'IE_OA', driver = 'ESRI Shapefile')
shp.ie <- ms_dissolve(shp.ie)
shp.ie <- ms_simplify(shp.ie, keep = 0.05)
colnames(shp.ie@data) <- c('id')
shp.ie@data$id <- 'IE'
try( do.call(file.remove, as.list(paste0(boundaries.path, '/', list.files(path = boundaries.path, pattern = '^IE\\.') ) ) ), silent = TRUE )
writeOGR(shp.ie, dsn = boundaries.path, layer = 'IE', driver = 'ESRI Shapefile')
df.ie <- tidy(shp.ie)
setnames(df.ie, c('long', 'lat'), c('X_lon', 'Y_lat'))
df.ie$type <- 'IE'
db_conn <- dbConnect(MySQL(), user = 'root', password = 'root', dbname = 'geography')
dbSendQuery(db_conn, paste0("DELETE FROM boundaries WHERE type = 'IE'") )
dbWriteTable(db_conn, 'boundaries', df.ie, row.names = FALSE, append = TRUE)
dbDisconnect(db_conn)

# HOW TO PLOT ALL COUNTRIES
shp.ctry <- readOGR(boundaries.path, layer = 'CTRY')
shp.ctry <- spRbind(shp.ie, shp.ctry)
shp.ctry
# OR
shp.ctry <- dbGetQuery(db_conn, "SELECT * FROM boundaries WHERE type IN ('CTRY', 'IE')")
ggplot(shp.ctry, aes(X_lon, Y_lat, group = group)) + geom_polygon(fill = 'darkred', colour = 'yellow') + coord_map()
# END


# # # STEP 3 - CALCULATE "METRICS" (using rgeos) --------------------------------------------------------------------------------

# Extract subareas codes
y <- shp.area@data

# Add (mass) CENTROIDS
y <- cbind(y, as.data.frame(gCentroid(shp.area, byid = TRUE) ) )

# Add PERIMETER (Length)
y <- cbind(y, as.data.frame(gLength(shp.area, byid = TRUE) ) )

# Add AREA
y <- cbind(y, as.data.frame(gArea(shp.area, byid = TRUE) ) )

# Other interesting functions:





# # # STEP 4 - PLOT MAPS  -------------------------------------------------------------------------------------------------------

# function to read and polish dataframe boundaries
load_boundaries <- function(cut.out.shetland = TRUE){
    db_conn <- dbConnect(MySQL(), user = 'root', password = 'root', dbname = 'geography')
    if(cut.out.shetland){
        bnd <- data.table( dbGetQuery(db_conn, 'SELECT * FROM boundaries WHERE Y_lat < 58.74') )
    } else {
        bnd <- data.table( dbReadTable(db_conn, 'boundaries') )
    }
    bnd_lkp <- data.table( dbReadTable(db_conn, 'boundaries_ids') )
    dbDisconnect(db_conn)
    setkeyv(bnd, c('id', 'type'))
    setkeyv(bnd_lkp, c('boundary_id', 'type'))
    bnd <- bnd[bnd_lkp]
    bnd[, id := NULL]
    setnames(bnd, 'i.id', 'id')
    bnd[, `:=`(
        order = as.integer(order),
        hole = as.logical(hole),
        piece = as.factor(piece),
        group = as.factor(group),
        id = as.factor(id),
        type = as.factor(type)
    )]
}
boundaries <- load_boundaries()

# read dataframe boundaries and create spatialPolygonDataframes (only for one-piece polygons)
create.poly <- function(area_bid, area_type){
    poly <- boundaries[id == area_bid & type == area_type & piece == 1, .(X_lon, Y_lat) ]
    Polygons( list(Polygon(poly)), bnd_lkp[boundary_id == area_bid & type == area_type, id] )
}
area.types <- c('CCG', 'LAT', 'NHSR', 'CCR', 'CTRY')
sp.polys <- list()
for (area.type in area.types){
    areas <- bnd_lkp[type == area.type & boundary_id %in% unique(boundaries[type == area.type, id]), .(boundary_id, id)]
    sp.poly <- lapply(areas$boundary_id, function(x) create.poly(x, area.type) )
    sp.poly <- SpatialPolygons(sp.poly, proj4string = CRS(proj.wgs) )
    sp.polys[[area.type]] <- SpatialPolygonsDataFrame(sp.poly, data = data.frame(row.names = areas$id, id = areas$id))
}

# # # STEP 4a - PLOT MAPS USING TMAP AND SHAPEFILES -----------------------------------------------------------------------------

#




tm_shape(
        shp = shape.name                    # defines the object to use for subsequent layer, there could be more than one call
        projection = proj.name              # allows to swap projections (as with coord_map() in ggplot2)
    ) +
        tm_borders(col = 'colname') +
    tm_fill(
        col = 'varname',                    # define which variable / metric to use to color-in the polygons and how
        style = 'quantile'
    ) +
    tm_bubbles(size = 'varname') +
    tm_raster(
        col = 'layer.name',                 # needed when the object to plot is a raster
        palette = 'pal.name',
    ) +
    tm_dots() +
    tm_text() +
    tm_polygons() +
    tm_lines() +
    tm_grid(n.x = , n.y = ) +               # adds equispaced longitude and latitude lines
    tm_legend(position = 'pos.name') +
    tm_compass() +
    tmap_style('style.name')  +
    tm_style_classic()


# # # STEP 4b - PLOT MAPS USING GGPLOT/GGMAP AND CSV/DB -------------------------------------------------------------------------

# retrieve boundaries, points. and and data
db_conn <- dbConnect(MySQL(), user = 'root', password = 'root', dbname = 'geography')
boundaries <- dbReadTable(db_conn, 'boundaries')
hospitals <- dbReadTable(db_conn, 'hospitals')
# dataset <- dbReadTable(db_conn, 'hospitals')
dbDisconnect(db_conn)

# initialize hospitals table for mapping
coordinates(hospitals) <- ~X_lon+Y_lat
# impose a projection on hospitals
proj4string(hospitals) <- CRS(proj.wgs)
# go back to dataframe/datatable
hospitals <- as.data.frame(hospitals)

# get some map
UK_map <- get_map(location = UK_centre, zoom = 6, scale = 2, source = 'stamen', maptype = 'toner')

# extract boundaries for requested area = 'XXX'
y <- boundaries[type = 'XXX']

# join boundaries and data to calculate var.fill
y <- merge(boundaries, dataset[, c('XXX', 'metric')], by.x = 'id', by.y = 'XXX')

# build tooltip
y$ttip <- ''

# build first layer
g <- ggmap(ggmap = UK_map, base_layer = ggplot(data = boundaries, aes(x = X_lon, y = Y_lat, group = group)) )

# add boundaries layer with fill reference, hover tooltip and reaction to click
# the alpha argument in the call allows the map to show through the polygons
g <- g + geom_polygon_interactive(aes(tooltip = ttip, data_id = ttip, fill = var.fill), colour = 'black', size = 0.2, alpha = 0.8)

# apply breaks to boundaries (equal: bin with same width, quantile: bin with same number of units, pretty: clear integer limits, fixed ?)


# apply palette style to boundaries (usually sequential, some times diverging, but never qualitative - unless there is nothing to display)


# add additional layer drawing thicker borders for outer areas


# join hospitals and data to calculate var.size or/and var.color


# build tooltip


# add hospitals layer
g <- g + geom_point_interactive(data = hospitals, aes(x = X_lon, y = Y_lat, size = var.size, colour = var.color, tooltip = ttip, data_id = ttip))

# apply breaks to hospitals


# apply palette style to hospitals


# g <- g + coord_map()   ???



# 2- MAPPING

# 2A- GGPLOT
# - operates on data frames, therefore polygons objects (shapefiles) need to be translated into a data frame format beforehand.
#   Maps in this format can then be plotted using the polygon geom, geom_polygon() draws lines between subsequent points and closes them up.
#   Of course, x = long and y = lat are the main aesthetics, but the "group" column have to be mapped to the "group" aesthetic to let ggplot understand when a polygon has to be closed, and start with another one
# - ggplot provides the map_data('mapname') function that loads a map found in the maps or mapdata package as a dataframe.
#   The structure of the resulting data frames is prepared as wanted with the following:
#    * long/lat are the coordinates: long is the longitud, and lat is the latitude
#    * order. This just shows in which order ggplot should connect the dots
#    * region and subregion tell what region or subregion a set of points surrounds.
#    * group. nt! ggplots functions can take a group argument which controls (amongst other things) whether adjacent points should be connected by lines. If they are in the same group, then they get connected, but if they are in different groups then they don’t.
#      Essentially, having to points in different groups means that ggplot lifts the pen when going between them.



ggplot() +
    geom_polygon(data = internal, aes(x = X_lon, y = Y_lat, group = group, fill = Y), color = 'black', size = 0.25) +
    geom_polygon(data = external, aes(x = X_lon, y = Y_lat, group = group), fill = NA, color = 'black', size = 0.40) +
    geom_point(data = points, aes(X_lon, Y_lat, group = id, fill = Z1, size = Z2), colour = 'black', size = 0.25) +
    coord_map(xlim = UK_bounds[c('lng2', 'lng1')], ylim = UK_bounds[c('lat2', 'lat1')])