library(sp) library(raster) library(rgdal) # Working with spatial data in R # Normally dealing with shapefiles (or sometimes geojson) # or tables with coordinates # The simplest data is a table with coordinates (i.e. point data) # here's a link to some malaria prevalence data from Burkina Faso BF_malaria_data <- read.table("https://www.dropbox.com/s/ioiv609g3g3nc5u/data_bf2_binomial.txt?dl=1", header=T) head(BF_malaria_data) # gives you the first few rows names(BF_malaria_data) # gives you the column names BF_malaria_data$examined # gives you the # examined column as a vector/string BF_malaria_data$positives # gives you the # positives column as a vector/string # If you want to create a new variable, such as # prevalence, you can use the $ sign BF_malaria_data$infection_prevalence <- BF_malaria_data$positives / BF_malaria_data$examined hist(BF_malaria_data$infection_prevalence) # Use R's basic plotting function to plot plot(BF_malaria_data$longitude, BF_malaria_data$latitude, ylab = "Latitude", xlab="Longitude") #boring! # Use the cex function to plot circle size as a function of a variable plot(BF_malaria_data$longitude, BF_malaria_data$latitude, cex=BF_malaria_data$infection_prevalence, ylab = "Latitude", xlab="Longitude") # Add some color prev_class <- as.numeric(cut(BF_malaria_data$infection_prevalence, 4)) prev_class color_class <- c("blue", "yellow", "orange", "red") color_class color_class[prev_class] # Use the cex function to plot circle size as a function of a variable plot(BF_malaria_data$longitude, BF_malaria_data$latitude, cex=BF_malaria_data$infection_prevalence, ylab = "Latitude", xlab="Longitude", col=color_class[prev_class]) # With filled circles plot(BF_malaria_data$longitude, BF_malaria_data$latitude, cex=BF_malaria_data$infection_prevalence, ylab = "Latitude", xlab="Longitude", col=color_class[prev_class], pch=16) # With larger filled circles plot(BF_malaria_data$longitude, BF_malaria_data$latitude, cex=BF_malaria_data$infection_prevalence*2, ylab = "Latitude", xlab="Longitude", col=color_class[prev_class], pch=16) # Let's make a SpatialPointsDataFrame object (useful for other operations) BF_malaria_data_SPDF <- SpatialPointsDataFrame(coords = BF_malaria_data[,c("longitude", "latitude")], data = BF_malaria_data[,c("examined", "positives", "infection_prevalence")], proj4string = CRS("+proj=longlat +datum=WGS84")) # WGS 1984 using lat/long. Optional but good to specify # can also use "+init=epsg:4326" # Summary of object BF_malaria_data_SPDF # SPDFs partition data elements, e.g. the coordinates are stored separately from the data BF_malaria_data_SPDF@coords BF_malaria_data_SPDF@data # You can quickly access the data frame as per a # standard data frame, e.g. BF_malaria_data_SPDF$infection_prevalence # You can use the plot or spplot function to get quick plots plot(BF_malaria_data_SPDF) spplot(BF_malaria_data_SPDF, zcol = "infection_prevalence") #### Spatial data types (e.g. shapefiles, geojson, raster) BF_Adm_1 <- readOGR("/Users/abennett1/Dropbox (Lao PDR - MEI)/SpatialEpiCourse", "BF_Adm_1") # Let's put it in context. The raster package allows you # to pull data from GADM # to get country info ccodes() BF_Adm_1 <-raster::getData('GADM', country="BFA", level=1) # let's take a look at that object BF_Adm_1 plot(BF_Adm_1) # plot both country and data points plot(BF_Adm_1) points(BF_malaria_data$longitude, BF_malaria_data$latitude, cex=BF_malaria_data$infection_prevalence*2, ylab = "Latitude", xlab="Longitude", col=color_class[prev_class], pch=16) # Plot over a web map library(leaflet) # leaflet allows you to layer # a basemap leaflet() %>% addTiles() # you can use one of several basemaps leaflet() %>% addTiles("https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}") leaflet() %>% addTiles("https://cartodb-basemaps-{s}.global.ssl.fastly.net/dark_all/{z}/{x}/{y}.png") # Let's choose a simple one basemap <- leaflet() %>% addTiles("https://cartodb-basemaps-{s}.global.ssl.fastly.net/light_all/{z}/{x}/{y}.png") # you can use the 'piping' command %>% to add layers # as our point and polygon data are already SP objects # this is easy basemap %>% addPolygons(data=BF_Adm_1) # Can you play with the way it plots # e.g. to change the colors/line weight basemap %>% addPolygons(data=BF_Adm_1, color = "red", weight = 1, fillOpacity = 0.2) # You can also add popups basemap %>% addPolygons(data=BF_Adm_1, popup = BF_Adm_1$NAME_1) # If you want to add points as well basemap %>% addPolygons(data=BF_Adm_1, weight = 2, popup = BF_Adm_1$NAME_1) %>% addCircleMarkers(data=BF_malaria_data_SPDF, color="red", radius = 2) #If you want to change the color according #to a variable, i.e. prevalence, you can use the # colorNumeric function library(oro.nifti) # for a nice color palette colorPal <- colorQuantile(tim.colors(), BF_malaria_data_SPDF$infection_prevalence, n = 4) # colorPal is now a function you can apply to get the corresponding # color for a value colorPal(0.1) basemap %>% addPolygons(data=BF_Adm_1, weight = 2, fillOpacity=0, popup = BF_Adm_1$NAME_1) %>% addCircleMarkers(data=BF_malaria_data_SPDF, color = colorPal(BF_malaria_data_SPDF$infection_prevalence), radius = 2, popup = as.character(BF_malaria_data_SPDF$infection_prevalence)) # Might want to add a legend. This just goes on as another layer # First calculate the labels. In this case, the quantiles of # prevalence prev_quantiles <- quantile(BF_malaria_data_SPDF$infection_prevalence, prob = seq(0,1,0.25)) prev_quantiles basemap %>% addPolygons(data=BF_Adm_1, weight = 2, fillOpacity=0, popup = BF_Adm_1$NAME_1) %>% addCircleMarkers(data=BF_malaria_data_SPDF, color = colorPal(BF_malaria_data_SPDF$infection_prevalence), radius = 2, popup = as.character(round(BF_malaria_data_SPDF$infection_prevalence,2))) %>% addLegend(colors=tim.colors(4), labels = c("<268%", "26 - 45%", "45 - 67%", ">67%")) # For more complex popups, you can define the html basemap %>% addPolygons(data=BF_Adm_1, weight = 2, fillOpacity=0, popup = BF_Adm_1$NAME_1) %>% addCircleMarkers(data=BF_malaria_data_SPDF, color = colorPal(BF_malaria_data_SPDF$infection_prevalence), radius = 2, popup = paste("

","Prevalence:", round(BF_malaria_data_SPDF$infection_prevalence,2), "

")) # You can plot rasters as well # you can use the getData function from the raster # package to get some environmental/climatological data elev <- raster("https://www.dropbox.com/s/kjkz22xc4cko7g4/elev_BFA.tif?dl=1") # Can also use the raster package getData function (very useful!) elev <- raster::getData("alt", country="BFA") elev # you can plot using the standard plot function plot(elev) # and you can use leaflet basemap %>% addRasterImage(elev) # If you want to add a legend, you have to define the color # palette first raster_colorPal <- colorNumeric(topo.colors(64), values(elev), na.color = NA) raster_colorPal(500) basemap %>% addRasterImage(elev, color = raster_colorPal) %>% addLegend(values = values(elev), pal = raster_colorPal) # If you want to export the data, there are several options. # 1. Export button - # 2. save as kml for someone to open in Google Earth library(plotKML) plotKML(BF_malaria_data_SPDF) # see ?plotKML for more options ## Now try creating a map using the Scottish lip cancer data # Hint: its an RData file called "Scotland.RData" in the dropbox folder... # you can use the 'load' function to load an RData object # You will be using this dataset later. The important variables are # 'COUNT' - number of cases, "PY" - person years, "EXP_" - expected number of cases # 'SMR' - Standardized morbidity ratio