Visualising Maps in R

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.2.3

## Warning: package 'ggplot2' was built under R version 4.2.3

## Warning: package 'tibble' was built under R version 4.2.3

## Warning: package 'tidyr' was built under R version 4.2.3

## Warning: package 'readr' was built under R version 4.2.3

## Warning: package 'purrr' was built under R version 4.2.3

## Warning: package 'dplyr' was built under R version 4.2.3

## Warning: package 'forcats' was built under R version 4.2.3

## Warning: package 'lubridate' was built under R version 4.2.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(here)

## Warning: package 'here' was built under R version 4.2.3

## here() starts at D:/john_snow

library(sf)

## Warning: package 'sf' was built under R version 4.2.3

## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(HistData)

## Warning: package 'HistData' was built under R version 4.2.3

library(devtools)

## Warning: package 'devtools' was built under R version 4.2.3

## Loading required package: usethis

## Warning: package 'usethis' was built under R version 4.2.3

devtools::install_github("yutannihilation/ggsflabel")

## Skipping install of 'ggsflabel' from a github remote, the SHA1 (846ede01) has not changed since last install.
##   Use `force = TRUE` to force installation

library(leaflet)

## Warning: package 'leaflet' was built under R version 4.2.3

#load data. assign using alt+- = <-

gbr <- st_read(here("data", "GBR_adm0.shp"))

## Reading layer `GBR_adm0' from data source `D:\john_snow\data\GBR_adm0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 70 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.69139 ymin: 49.86542 xmax: 1.764168 ymax: 61.52708
## Geodetic CRS:  WGS 84

gbr1 <- st_read(here("data", "GBR_adm1.shp"))

## Reading layer `GBR_adm1' from data source `D:\john_snow\data\GBR_adm1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.69139 ymin: 49.86542 xmax: 1.764168 ymax: 61.52708
## Geodetic CRS:  WGS 84

gbr2 <- st_read(here("data", "GBR_adm2.shp"))

## Reading layer `GBR_adm2' from data source `D:\john_snow\data\GBR_adm2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 192 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13.69139 ymin: 49.86542 xmax: 1.764168 ymax: 61.52708
## Geodetic CRS:  WGS 84

deaths <- HistData::Snow.deaths 
streets <- HistData::Snow.streets
pumps <- HistData::Snow.pumps


#plot the data
gbr %>% ggplot() + geom_sf()

gbr1 %>% ggplot() + geom_sf()

gbr2 %>% ggplot() + geom_sf()

wm <- gbr2 %>% filter(NAME_2 == "Westminster")
wm %>% ggplot() + geom_sf()

plot(gbr1)

plot(gbr2)

## Warning: plotting the first 10 out of 11 attributes; use max.plot = 11 to plot
## all

#located westminster and plotted that next
plot(wm)

## Warning: plotting the first 9 out of 11 attributes; use max.plot = 11 to plot
## all

st_crs(wm)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

st_crs(gbr1)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

st_crs(gbr2)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

st_crs(gbr)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

#before making it as spatial data, trying to understand and modify data as per use
names(deaths)

## [1] "case" "x"    "y"

names(streets)

## [1] "street" "n"      "x"      "y"

names(pumps)

## [1] "pump"  "label" "x"     "y"

names(wm)

##  [1] "ID_0"      "ISO"       "NAME_0"    "ID_1"      "NAME_1"    "ID_2"     
##  [7] "NAME_2"    "TYPE_2"    "ENGTYPE_2" "NL_NAME_2" "VARNAME_2" "geometry"

#converting as sf
deaths_sf <- st_as_sf(deaths,coords = c("x","y"),crs = 4326)
streets_sf <- st_as_sf(streets,coords = c("x","y"),crs = 4326)
pumps_sf <- st_as_sf(pumps,coords = c("x","y"),crs = 4326)

deaths_sf %>% ggplot() + geom_sf()

pumps_sf %>% ggplot() + geom_sf()

streets_sf %>% ggplot() + geom_sf()

#streets are also plotted as points, so converting them to lines using linestring

streets_sf1<- streets_sf%>% 
  group_by(street) %>% 
  summarise(n = mean(n)) %>% 
  st_cast('LINESTRING')

#plotting the new line strret map
streets_sf1 %>% ggplot() + geom_sf()

st_crs(streets_sf1)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

#crs was 4326, same as others

ggplot() +
  geom_sf(data = streets_sf1, color='orange') +
  geom_sf(data = pumps_sf,color='black', alpha = 1 ,cex=3) +
geom_sf(data = deaths_sf,color='red', alpha=0.5,size = 2) +
  ggsflabel::geom_sf_label(data = pumps_sf, aes(label = label),
                           nudge_x=0.025,
                            nudge_y=-1)

## Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
## give correct results for longitude/latitude data

pumps_sf1 <- pumps_sf %>% st_transform(32453)
pumps_sf1_buf <- st_buffer(pumps_sf1,dist = 500000)
streets_sf2 <- st_transform(streets_sf1, crs = 32643)

streets_sf2_buf <- st_buffer(streets_sf2, dist = 90000)
streets_sf1_buf <- st_buffer(streets_sf1, dist = 40000)

ggplot() +
  geom_sf(data = pumps_sf1,color='blue') +
  geom_sf(data = pumps_sf1_buf, color = 'red', alpha = 0.5, size = 2)

ggplot() +
  geom_sf(data = pumps_sf1,color='blue') +
  geom_sf(data = pumps_sf1_buf%>% 
            filter(label=="Broad St"), color = 'red', alpha = 0.2)

ggplot() +
  geom_sf(data = streets_sf1) +
  geom_sf(data = deaths_sf, color = 'red', alpha = 0.5, size = 2)

ggplot() +
  geom_sf(data = streets_sf1) +
  geom_sf(data = deaths_sf, color = 'red', alpha = 0.5, size = 2) +
  geom_sf(data = pumps_sf , shape = 22, size = 4, fill = 'blue', color = 'blue')+
    ggsflabel::geom_sf_label(data = pumps_sf, aes(label = label),
                             nudge_x=0.025,
                             nudge_y=-1)+
  labs(title = "John Snow map , Broad St highlighted",
       subtitle = "with streets , pumps and deaths plotted")+
  ggspatial::annotation_scale(location = "bl")+
  ggspatial::annotation_north_arrow(location = "tr")+
  xlab("Longitude")+
  ylab("Latitude")+
  geom_text(aes(x = 2.0, y = 10.0, label = "Credits: Histdata"),
             color = "black", check_overlap = F, size = 3)

## Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
## give correct results for longitude/latitude data

#other modifications tried

   ggplot()+
 geom_sf(data = streets_sf1, color='orange') +
    geom_sf(data = pumps_sf,color='black', alpha = 1 ,cex=3) +
    geom_sf(data = deaths_sf,color='red', alpha=0.5,size = 2) +
    geom_sf(data = pumps_sf1_buf%>% 
              filter(label=="Broad St"), color = 'black', alpha = 0.3)+
    ggsflabel::geom_sf_label(data = pumps_sf, aes(label = label),
                             nudge_x=0.025,
                             nudge_y=-1)+
  labs(title = "Broad street map",
       subtitle = "with streets , pumps and deaths plotted")+
  ggspatial::annotation_scale(location = "bl")+
  ggspatial::annotation_north_arrow(location = "tr")+
  xlab("Longitude")+
  ylab("Latitude")+
  geom_text(aes(x = 2.0, y = 10.0, label = "Credits: Histdata"),
             color = "red", check_overlap = F, size = 3)

## Warning in st_point_on_surface.sfc(data$geometry): st_point_on_surface may not
## give correct results for longitude/latitude data

#Adding Title/ sub-title.
#Describing Distance and Scale.
#Direction symbols.
#Adding axis labels
#Acknowledging source

gbr2 %>% ggplot()+
  geom_sf()+
  labs(title = "United Kingdom map",
       subtitle = "For PHT 2022 (AMCHSS) module")+
  ggspatial::annotation_scale(location = "bl")+
  ggspatial::annotation_north_arrow(location = "tr")+
  xlab("Longitude")+
  ylab("Latitude")+
  geom_text(aes(x = 0.0, y = 58, label = "Credits: DIVAGIS"),
             color = "orange", check_overlap = F, size = 3)

## Scale on map varies by more than 10%, scale bar may be inaccurate

wm %>% ggplot()+
  geom_sf()+
  labs(title = "Westminster map",
       subtitle = "For PHT 2022 (AMCHSS) module")+
  ggspatial::annotation_scale(location = "bl")+
  ggspatial::annotation_north_arrow(location = "tr")+
  xlab("Longitude")+
  ylab("Latitude")+
  geom_text(aes(x = 0.00, y = 51.50, label = "Credits: DIVAGIS"),
             color = "red", check_overlap = F, size = 3)

#Visualising geometry of sf object
deaths_sf %>% 
  ggplot()+
  geom_sf()+
  labs(title = "Deaths due to cholera")+
  xlab("Longitude")+
  ylab("Latitude")

pumps_sf1 %>% 
  ggplot()+
  geom_sf()+
  labs(title = "Location of pumps")+
  xlab("Longitude")+
  ylab("Latitude")

streets_sf1 %>% 
  ggplot()+
  geom_sf()+
  labs(title = "Streets around the area")+
  xlab("Longitude")+
  ylab("Latitude")

# Trying to incorporate westminster map to leaflet

library(leaflet)
basemap <- leaflet() %>% addTiles()
basemap

basemap %>% addPolygons(data = gbr,
                        color = "orange",
                        fillOpacity = 0.1,
                        opacity = 0.5,
                        weight = 1)

basemap %>% addPolygons(data = wm,
                        color = "blue",
                        fillOpacity = 0.1,
                        opacity = 0.5,
                        weight = 1)

#trying another combination

basemap %>% addPolygons(data = wm,
                        color = "green",
                        fillOpacity = 0.1,
                        opacity = 0.5,
                        weight = 1,
                        popup = wm$block)

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:

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.