library(sp) library(raster) library(leaflet) library(oro.nifti) library(rgdal) library(geosphere) # Manipulating and relating spatial data # Let's get admin 1 data for Burkina Faso BF_Adm_1 <- raster::getData("GADM", country="BFA", level=1) # You can subset a SpatialPolygonsDataFrame just like a data frame #i.e. BF_Adm_1_cropped <- BF_Adm_1[1,] BF_Adm_1_cropped plot(BF_Adm_1) lines(BF_Adm_1_cropped, col="red", lwd=2) # You can also subset by name BF_Adm_1_Cascades <- subset(BF_Adm_1, BF_Adm_1$NAME_1=="Cascades") # BF_Adm_1[BF_Adm_1$NAME_1=="Cascades",] will also work BF_Adm_1_Cascades plot(BF_Adm_1) lines(BF_Adm_1_Cascades, col="blue", lwd=2) # Now read in the BF malaria data # This one has no covariates (we'll get these ourselves) BF_malaria_data <- read.csv("https://www.dropbox.com/s/bfs3pinxe1lvvxr/data_bf2_binomial.csv?dl=1", header=T) # Convert to a SPDF BF_malaria_data_SPDF <- SpatialPointsDataFrame(coords = BF_malaria_data[,c("longitude", "latitude")], data = BF_malaria_data[,c("examined", "positives")], proj4string = CRS("+init=epsg:4326")) # To identify which admin unit each point lies in # you can use the 'over' function over(BF_malaria_data_SPDF, BF_Adm_1_Cascades) # The CRS isn't identical, so over kicks out an error. # To transform projections you can use the spTransform function BF_malaria_data_SPDF <- spTransform(BF_malaria_data_SPDF, crs(BF_Adm_1)) BF_malaria_data_SPDF # Now try again BF_Adm_1_per_point <- over(BF_malaria_data_SPDF, BF_Adm_1) # Now we can use this to calculate admin unit specific statistics # We might be interested in the number of sites per admin unit table(BF_Adm_1_per_point$NAME_1) # or we can use the tapply function for more complex calculations Nex_per_Adm1 <- tapply(BF_malaria_data_SPDF$examined, BF_Adm_1_per_point$NAME_1, sum) Npos_per_Adm1 <- tapply(BF_malaria_data_SPDF$positives, BF_Adm_1_per_point$NAME_1, sum) prev_per_Adm1 <- Npos_per_Adm1 / Nex_per_Adm1 # Now we can use this to make a map of prevalence colorPal <- colorQuantile(tim.colors(), prev_per_Adm1, n = 4) plot(BF_Adm_1, col=colorPal(prev_per_Adm1)) #or leaflet() %>% addTiles() %>% addPolygons(data=BF_Adm_1, col=colorPal(prev_per_Adm1), fillOpacity=0.6) # You can also define your own color bins colorPal <- colorBin(tim.colors(), prev_per_Adm1, bins = c(0, 0.25, 0.5, 0.75, 1)) plot(BF_Adm_1, col=colorPal(prev_per_Adm1)) # Let's get some other layers # Elevation BF_elev <- raster::getData("alt", country="BF") plot(BF_elev) # Land use (# see http://due.esrin.esa.int/files/GLOBCOVER2009_Validation_Report_2.2.pdf) BF_land_use <- raster("https://www.dropbox.com/s/hc9m6ac3kb845ip/BF_land_use.tif?dl=1") # Also on the CLE BF_land_use plot(BF_land_use) table(BF_land_use[]) # For a break down of the classes in BF # Its good practice to resample rasters to the same extent and resolution (i.e. same grid) # This makes it easier to deal with later and to relate rasters to each other # The resample command makes this process easy # The default method is bilinear interpolation, which doesn't make sense for our categorical # variable, so we can use the nearest neighbour function 'ngb BF_land_use_resampled <- resample(BF_land_use, BF_elev, method="ngb") # AH - why might they not intercect?? Hint: check the projections... crs(BF_land_use) # Mercator crs(BF_elev) # WGS84 # To reproject a raster, you can use the projectRaster function BF_land_use <- projectRaster(BF_land_use, BF_elev, method="ngb") # Now try resampling BF_land_use_resampled <- resample(BF_land_use, BF_elev, method="ngb") BF_land_use_resampled;BF_elev # You can manipulate rasters in many ways # you might want to change the resolution for analysis res(BF_elev) # in decimal degrees. 1 dd roughly 111km at the equator BF_elev_low_res <- aggregate(BF_elev, fact = 10) # by default, calculates mean res(BF_elev_low_res) plot(BF_elev_low_res) # You can crop them using another spatial layer BF_elev_Cascades <- crop(BF_elev, BF_Adm_1_Cascades) plot(BF_elev_Cascades) lines(BF_Adm_1_Cascades) # You can change their values plot(BF_elev*3.28) # in feet # to recategorize plot(cut(BF_elev, 4)) #or create a new raster and change its values BF_elev_cat <- BF_elev BF_elev_cat[] <- ifelse(BF_elev[]>250,1,2) plot(BF_elev_cat) # If a raster is the same resolution and extent, you can # perform joint operations on them, e.g. plot(BF_elev - BF_land_use_resampled) # Meaningless! Just for illustrative purposes.. # Now let's extract values of elevation at each survey point # you can use the extract function from the raster package extract(BF_elev, BF_malaria_data_SPDF) # We can catch this as an object, or add it directly to our SPDF BF_malaria_data_SPDF$elev <- extract(BF_elev, BF_malaria_data_SPDF) BF_malaria_data_SPDF# now has 3 variables # We can now have a quick look at the relationship between prevalence # and elevation # First generate a prevalence variable BF_malaria_data_SPDF$prevalence <- BF_malaria_data_SPDF$positives / BF_malaria_data_SPDF$examined plot(BF_malaria_data_SPDF$elev, BF_malaria_data_SPDF$prevalence) library(stats) scatter.smooth(BF_malaria_data_SPDF$elev, BF_malaria_data_SPDF$prevalence) # with lowes # You can also extract values using polygons # e.g to get admin 1 level elevations # You just have to define a function to apply, otherwise you get all the pixel # values per polygon BF_Adm_1$elev <- extract(BF_elev, BF_Adm_1, fun=mean, na.rm=TRUE) # takes a little longer.. # you might also be interested in distances to/from other features (e.g. health facilities, water) # e.g. distance to nearest river # Get a rivers file (obtained via http://www.diva-gis.org/Data) waterbodies <- readOGR("/Users/hughsturrock/Dropbox/SpatialEpiCourse/", "BFA_water_areas_dcw") waterbodies plot(waterbodies) # The goesphere package has some nice functions such as (takes a little while to compute) dist_to_water <- dist2Line(BF_malaria_data_SPDF, waterbodies) # Look at the result dist_to_water # Can add to your data frame BF_malaria_data_SPDF$dist_to_water <- dist_to_water[,1] # Questions # 1. Try making a map of infection prevalence at admin 2 level # 2. Use this file of population density to calculate a table of # numbers of people infected per admin 2 level BF_pop <- raster("https://www.dropbox.com/s/a9glj1is86o0xvz/BF_pop.tif?dl=1") # originally from worlpop.org # 3. use the getData function to get hold of more bioclimatic variables. # Add 2-3 more variables to each survey point and explore relationships with plots bioclim_global <- raster::getData('worldclim', var='bio', res=10) # low res but easy to work with for the moment bioclim_global # is a raster stack. Bascially a cube of rasters plot(bioclim_global) # See http://www.worldclim.org/bioclim for more details on these layers # to get the first raster 'bio1' bioclim_global[["bio1"]] # bioclim_global[[1]] will also work plot(bioclim_global[["bio1"]]) # 4. reclassify the BF_land_use layer according to 'LCCS Entry' column at the bottom of the # doc http://due.esrin.esa.int/files/GLOBCOVER2009_Validation_Report_2.2.pdf # And extract values. plot prevalence by category # 5. Calculate distance to nearest health facility in km for each survey point # Health facility shapefile available here # https://data.humdata.org/dataset/burkina-faso-healthsites # Hint1: the distm function generates a distance matrix between all pairs of points # Hint2: the apply() function can find the minimum value across rows or columns # 6. calculate distance (euclidean) to nearest pixel with >200 people # Hint: coordinates(BF_elev) gives you a matrix of all coordinates of the BF_elev raster # 7. Optional extra # Calc numbers positive per admin 2 - HINT: not all admin 2 have data points in them # Repeat on a different point dataset ## Assignment # Find 2 more interesting layers/variables (not from the above) to relate to the survey data . # You can find a raster and extract, or generate your own variable (i.e. distance to something)