1 Load packages

## Load packages

library(tidyverse) #for visualization (ggplot2) and data wrangling (dplyr)
library(rgdal) # to read shape file (.shp)
library(janitor) #to clean data frames
library(broom) #to convert shp --> data frame

2 Load map (geospatial) data

## Load the shape file using readOGR function of rgdal package


lk_data <- rgdal::readOGR("D:/TBP_Data Analysis/Maps/data/polbnda_lka.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\TBP_Data Analysis\Maps\data\polbnda_lka.shp", layer: "polbnda_lka"
## with 214 features
## It has 9 fields
## Check the limits of the shapefile (i.e. extent of the coordinates)
summary(lk_data)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min       max
## x 79.519148 81.877154
## y  5.919079  9.834897
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=GRS80 +no_defs]
## Data attributes:
##     f_code              coc                nam                laa           
##  Length:214         Length:214         Length:214         Length:214        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##       pop              ypc      adm_code             salb          
##  Min.   :-1e+08   Min.   :0   Length:214         Length:214        
##  1st Qu.:-1e+08   1st Qu.:0   Class :character   Class :character  
##  Median :-1e+08   Median :0   Mode  :character   Mode  :character  
##  Mean   :-1e+08   Mean   :0                                        
##  3rd Qu.:-1e+08   3rd Qu.:0                                        
##  Max.   :-1e+08   Max.   :0                                        
##      soc           
##  Length:214        
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
## Check how many divisions are in the shape file
length(lk_data)
## [1] 214
#Check the shape file meta data
head(lk_data@data)
##   f_code coc           nam          laa       pop ypc adm_code salb soc
## 0  FA001 LKA       Eastern       Ampara -99999999   0      UNK  UNK LKA
## 1  FA001 LKA North Western   Kurunegala -99999999   0      UNK  UNK LKA
## 2  FA001 LKA       Central       Matale -99999999   0      UNK  UNK LKA
## 3  FA001 LKA  Sabaragamuwa   Rathnapura -99999999   0      UNK  UNK LKA
## 4  FA001 LKA North Central Anuradhapura -99999999   0      UNK  UNK LKA
## 5  FA001 LKA       Central        Kandy -99999999   0      UNK  UNK LKA

3 Load and clean Epi data

#load data
tick_prev <- read.csv("D:/TBP_Data Analysis/Maps/tbp_prev_lk.csv")

#check data
head(tick_prev)
##       District Bgibsoni   Bvogeli   Hcanis            zone  X
## 1       Ampara 20.00000 10.000000 20.00000 low_country_dry NA
## 2 Anuradhapura 20.00000 10.000000 20.00000 low_country_dry NA
## 3      Badulla 27.36156  2.931596 17.58958  up_country_wet NA
## 4   Batticaloa 20.00000 10.000000 20.00000 low_country_dry NA
## 5      Colombo 33.33333 10.666667 32.00000 low_country_wet NA
## 6        Galle 33.33333 10.666667 32.00000 low_country_wet NA
#Clean data
tick_prev <- tick_prev%>%
  select(1:5)

#Arrange data
tick_prev_pivot <- tick_prev%>%
    pivot_longer(cols = c("Bgibsoni","Bvogeli","Hcanis")  ,  names_to = "pathogen", values_to = "prevalence")

4 Merge geospatial and epi data

#Convert .shp file into a data frame

lk_data_fortified <- tidy(lk_data, region = "laa")

#Merge with epi data frame

all_data <- left_join(lk_data_fortified, tick_prev_pivot,by = c("id" = "District"))

For the conversion of the shape file note down the column that indicate the boundaries of interest. e.g., districts, countries etc.

For merging the two data frame note down the column names of the common data e.g., districts, countries etc.

5 Plot with ggplot2

plot_map <- ggplot() +
  geom_polygon(data = all_data, aes( x = long, y = lat, group = group, fill= prevalence)) + #Use fill option to plot the data in
  scale_fill_gradient(low ="#B3D6DB", high = "#0D343D" )+ #Set the colour to gradient
  theme_void() + #Use this theme to remove axes
  facet_wrap(~pathogen) # Plot individual maps for each pathogen

print(plot_map) #Visualize the map