## It is nice practice to set up a working directory where your data is stored/will output ## note: use "/" instead of "\", or "\\" wd <- "CHANGE TO YOUR DIRECTORY [FILE WITH TUTORIAL DATA]" # working directory setwd(wd) # set the working directory getwd() # to get your working directory details #load needed packages library(sf) library(dplyr) library(ggplot2) library(tmap) library(mapview) library(tidyverse) ################################### # OTAGO SUMMER SCHOOL TURTORIAL 3 # ################################### ######## Tasks and goals ######## # # # make buffers # # summarize points in a polygon # # calculate distance # # # ################################## ##Bring in greenspace shp and make buffers greenspace <- st_read("CHANGE HERE TO NAME OF SHP") greenspace <- greenspace %>% st_transform(2193) ##we want to make buffers in meters so make sure it is in metric coordinate system greenspace$unique_ID <- seq.int(nrow(greenspace)) #make a unique id for greenspaces # Make 200 and 500meter buffers around the greenspace greenspace_200 <- st_buffer(greenspace, 200) greenspace_500 <- st_buffer(greenspace, 500) #visualize the greenspace and buffers mapview(greenspace[1,]) + mapview(greenspace_200[1,]) + mapview(greenspace_500[1,]) ## Bring in obesity data and plot obesity <- st_read("CHANGE HERE TO NAME OF SHP") obesity <- obesity %>% st_transform(2193) obesity$unique_ID_obesity <- seq.int(nrow(obesity)) #make a unique id for obesity cases mapview(obesity) #do a spatial join:one-to-many [buffers can overlap] obesity_greenspace_200 <- st_intersection(greenspace_200, obesity) obesity_greenspace_500 <- st_intersection(greenspace_500, obesity) # count by park spjoin_200 <- obesity_greenspace_200 %>% group_by(unique_ID, unique_ID_obesity) %>% count () obesity_per_greenspace_200 <- spjoin_200 %>% group_by(unique_ID) %>% summarise(obesity_per_greenspace = sum(n)) summary(obesity_per_greenspace_200) ### get descriptive stats on number of obesity cases within 200m of greenspace spjoin_500 <- obesity_greenspace_500 %>% group_by(unique_ID, unique_ID_obesity) %>% count () obesity_per_greenspace_500 <- spjoin_500 %>% group_by(unique_ID) %>% summarise(obesity_per_greenspace = sum(n)) summary(obesity_per_greenspace_500) ### get descriptive stats on number of obesity cases within 500m of greenspace ## Bring in food outlet data food_outlet <- st_read("CHANGE HERE TO NAME OF SHP") ## Assess distance from food outlet to obesity case food_outlet_w_distance <- st_distance(obesity, food_outlet) fartherest_food_outlet <- max(food_outlet_w_distance) print(fartherest_food_outlet) ### see the greatest distance between food outlet and obesity case nearest_food_outlet <- min(food_outlet_w_distance) print(nearest_food_outlet) ### see the how close the nearest food outlet is to an obesity case ################################################ # OTAGO SUMMER SCHOOL TURORIAL 3 BONUS (#2) # ################################################ ######## Tasks and goals ######## # # # protect sensitive data # # jitter points # # alter display # # # ################################## ## Make new df of obesity cases to not accidentally jitter original file jitter_obesity <- obesity jittered <- st_jitter(jitter_obesity, amount = 200) #amount is the distance in meters that points will be shifted #compare jittered with original data mapview(obesity, col.regions = "darkolivegreen", col = "darkolivegreen") + mapview(jittered, col.regions = "red", col = "red") mapview(obesity[1,], col.regions = "darkolivegreen", col = "darkolivegreen") + mapview(jittered[1,], col.regions = "red", col = "red") #to just look at one case ################################### # OTAGO SUMMER SCHOOL TURTORIAL 4 # ################################### ######## Tasks and goals ######## # # # Kernel density mapping # # # ################################## install.packages("spatstat") install.packages("terra") install.packages("spData") install.packages("adehabitatHR") if (!require("rspat")) remotes::install_github("rspatial/rspat") library(spatstat) library(terra) library(spData) library(raster) library(adehabitatHR) library(rspat) library(sf) library(mapview) ## Now bring in your data gambling <- st_read("CHANGE HERE TO NAME OF SHP") #since your working directory was set you just need your shp file name gambling <- st_as_sf(gambling, coords = c("X", "Y"), crs = 2193) ## Bring in CAU shp to restrict the map to the auckland area auckland_cau <- st_read("CHANGE HERE TO NAME OF SHP") #since your working directory was set you just need your shp file mapview(auckland_cau, col.regions = "grey", col = "black") + mapview(gambling,col.regions = "hotpink", col = "hotpink") ## Assess distance between gambling locations gambling_w_distance <- st_distance(gambling) fartherest_gambling <- max(gambling_w_distance) print(fartherest_gambling) #greatest distance between gambling outlets nearest_gambling <- min(gambling_w_distance) print(nearest_gambling) #smallest distance between gambling outlets hist(gambling_w_distance) #histogram of distance between gambling outlets #make kernel density maps ##project to meters coordinate system auckland_cau_2193 <- st_transform(auckland_cau, 2193) gambling_2193 <- st_transform(gambling, 2193) ## need to first convert to ppp [point pattern object] and set the analysis window to Auckland auckland_union <- st_union(auckland_cau_2193) auck_window <- as.owin(auckland_union) gambling_xy <- st_coordinates(gambling_2193) gambling_ppp <- ppp(x = gambling_xy [,1], y = gambling_xy [,2], window = auck_window) plot(gambling_ppp) ## Calculate different kernel density bandwidths. sigma is meters because of the projected coordinate system set above kde_50m<- density.ppp(gambling_ppp, sigma = 1000, eps = 50, edge = TRUE) plot(kde_50m) plot(auckland_union, add = TRUE) # Make the spatstat object a terra raster, then transform to mapview to make a prettier the map kde_r_50m <- rast(kde_50m) crs(kde_r_50m) <- "EPSG:2193" mapview(auckland_cau_2193, col.regions = "transparent", alpha.regions = 0, col = "black") + #this will allow you to click on the region and see the details mapview(kde_r_50m, layer.name = "KDE Gambling Places (50 m cells, sigma = 1000 m)", alpha.regions = 0.9, na.color = "transparent") #the legend is all zero becuase the values are in scientific notation kde_km2 <- kde_r_50m * 1e6 #convert them mapview(auckland_cau_2193, col.regions = "transparent", alpha.regions = 0, col = "black") + mapview(kde_km2, layer.name = "KDE Gambling Places (per km2; 50m cells)", alpha.regions = 0.9, na.color = "transparent") ############################## ##Now try a different cell size kde_737m<- density.ppp(gambling_ppp, sigma = 1000, eps = 737, edge = TRUE) plot(kde_737m) kde_r_737m <- rast(kde_737m) crs(kde_r_737m) <- "EPSG:2193" mapview(auckland_cau_2193, col.regions = "transparent", alpha.regions = 0, col = "black") + mapview(kde_r_737m, layer.name = "KDE Gambling Places (737 m cells, sigma = 1000 m)", alpha.regions = 0.9, na.color = "transparent") #again the legend is all zero becuase the values are in scientific notation kde_km2_737 <- kde_r_737m * 1e6 #convert them mapview(auckland_cau_2193, col.regions = "transparent", alpha.regions = 0, col = "black") + mapview(kde_km2_737, layer.name = "KDE Gambling Places (per km2; 737m cells)", alpha.regions = 0.9, na.color = "transparent") ################ #Now do KDE for number of machines marks(gambling_ppp) <- gambling_2193$MACHINES ##link the variable to the ppp object marks(gambling_ppp) <- as.numeric(marks(gambling_ppp)) #ensure it is numeric kde_machines <- density.ppp(gambling_ppp, sigma = 1000, eps = 50, weights = marks(gambling_ppp), edge= TRUE) #same thing again, make the ppp into a raster, then into a mapview object kde_w_machines <- rast(kde_machines) crs(kde_w_machines) <- "EPSG:2193" #we know legend will be all zeros so first re-scale kde_w_machines <- kde_w_machines * 1e6 mapview(kde_w_machines, layer.name = "KDE Gambling Machines (per km2)", alpha.regions = 0.9, na.color = "transparent") ################################### # OTAGO SUMMER SCHOOL TURTORIAL 5 # ################################### ######## Tasks and goals ######## # # # Spatial autocorrelation # # # ################################## ## Load libraries and data files library(sp) # for working with spatial data library(spdep) # a big library of spatial functions library(tmap) #help(package=spdep) # lists general help for this library ## Bring in your data alcohol <- st_read("CHANGE HERE TO NAME OF SHP") #since your working directory was set you just need your shp file name crime <- st_read("CHANGE HERE TO NAME OF SHP") #since your working directory was set you just need your shp file ##Global Moran'I process ## Need to first build elements for calculation crime_coords <- st_centroid(st_geometry(crime)) #use polygon centroids for the distance calculation crime_nb_1000 <- dnearneigh(crime_coords, d1 = 0, d2 = 1000) #use a neighborhood threshold of 1000m dist_list <- nbdists(crime_nb_1000, crime_coords) ## calculate Euclidean distances b/w neighbors idw_list <- lapply(dist_list, function(x) 1/x) #create idw weights crime_w_idw <- nb2listw(crime_nb_1000, glist = idw_list, style = "W", zero.policy = TRUE) #W is row standardized option moranI_crime <- moran(crime$AvgAllCrim, crime_w_idw, length(crime_nb_1000), Szero(crime_w_idw))[1]#Calculate Global Moran's I print(moranI_crime) ## BUILD NEIGHBOR MATRIX ## # K-nearest Neighbors ## # using knearneigh() in spdep; k= number of nearest neighbors to calculate nn8 <- knearneigh(crime_coords, k=8) # Creates a matrix of nn indexes using 8 nearest neighbors crime8.nlist <- knn2nb(nn8, row.names = NULL, sym = FALSE) # create neighbor list crime8.sw <- nb2listw(crime8.nlist) # create a spatial weights matrix #Calculate autocorrelation statistics moran.test(crime$AvgAllCrim, crime8.sw) #now add spatial weights crime8.dist <- nbdists(crime8.nlist, crime_coords) # calculates distance # now invert distance to determine weights (closer=higher) crime8.dist1 <- lapply(crime8.dist, function(x) {ifelse(is.finite(1/x), (1/x), (1/0.001))}) crime8.dist2 <- lapply(crime8.dist, function(x) {ifelse(is.finite(1/x^2), (1/x^2), (1/0.001^2))}) summary(unlist(crime8.dist1)) # check out the distance distribution # now create sp. weights matrix weighted on distance crime8.d1sw <- nb2listw(crime8.nlist, glist=crime8.dist1) crime8.d2sw <- nb2listw(crime8.nlist, glist=crime8.dist2) moran.test(crime$AvgAllCrim, crime5.d1sw) moran.test(crime$AvgAllCrim, crime5.d2sw) # Print distances and weights distwts <- cbind(c(sort(crime8.dist[[1]]),sort(crime8.dist[[2]]),sort(crime8.dist[[3]])), c(sort(crime8.sw$weights[[1]]),sort(crime8.sw$weights[[2]]),sort(crime8.sw$weights[[3]])), c(sort(crime8.d1sw$weights[[1]],T),sort(crime8.d1sw$weights[[2]],T),sort(crime8.d1sw$weights[[3]],T)), c(sort(crime8.d2sw$weights[[1]],T),sort(crime8.d2sw$weights[[2]],T),sort(crime8.d2sw$weights[[3]],T))) print(distwts[1:5,]) print(distwts[6:10,]) print(distwts[11:15,]) #calculate Anselin Local Moran's I lmoran <- localmoran(crime$AvgAllCrim, crime8.sw) crime$LocalI <- lmoran[,"Ii"] # the Local Moran's I statistic for each area crime$LocalIz <- lmoran[,"Z.Ii"] # z-scores for each local I crime$LocalIp <- lmoran[,"Pr(z != E(Ii))"] # p-value for the local Moran # Prep the output for plotting x <- as.numeric(scale(crime$AvgAllCrim)) lag_x <- lag.listw(crime8.sw, x) alpha <- 0.05 #set significance for p value sig <- crime$LocalIp <= alpha # create a LISA cluster column # Assign each cluster to a category crime$lisa_cluster <- "Insignificant" crime$lisa_cluster[sig & x > 0 & lag_x > 0] <- "High-High" crime$lisa_cluster[sig & x < 0 & lag_x < 0] <- "Low-Low" crime$lisa_cluster[sig & x > 0 & lag_x < 0] <- "High-Low" crime$lisa_cluster[sig & x < 0 & lag_x > 0] <- "Low-High" #ensure the order of the legend crime$lisa_cluster <- factor( crime$lisa_cluster, levels = c("High-High","Low-Low","High-Low","Low-High","Insignificant") ) #plot the map tm_shape(crime) + tm_polygons("lisa_cluster", palette = c("High-High" = "red", "Low-Low" = "blue", "Low-High" = "lightblue", "High-Low" = "pink", "Insignificant" = "white"), title = "LISA Clusters") + tm_layout(main.title = "LISA for Crime", legend.outside = TRUE)