In today’s lab, we will be looking at some of the more fundamental exploratory analysis tools. Cluster Analysis (using K-Means), and Quadrat Analysis,

## K-Means Clustering

Clustering algorithms are useful for separating data into groups with similar qualities, clustering is often not the end goal in and of itself, but it does provide a means to start looking for patterns, once groups have been formed it is then possible to take a deeper look into what is driving those clusters allowing us to find patterns * in the attributes*.

Below is a demo of the K-Means algorithm created using shiny, much like we made our maps in lab 5. This demo allows you to step through the iterations to see the impacts.

https://connect.gis.unbc.ca/content/8453154e-321e-48eb-9304-7375277ae188/

In this example, we are working in 2D space that is to say 2 attributes, and distance is calculated as the distance between point p, and cluster k is.

`\[d\ =\ \sqrt{\left(p.x-k.x\right)^2\ +\left(p.y\ -\ k.y\right)^2}\]`

This is relatively easy to understand, when moving into n columns this formula will change to

`\[\mathrm{d=}\sqrt{\sum_{\mathrm{i=1}}^{\mathrm{n}}\mathrm{p.}\mathrm{i}^{\mathrm{2}}\mathrm{-k.}\mathrm{i}^{\mathrm{2}}}\]`

As a refresher, we will do the same exercise we did in class last week.

Open RStudio Workbench (https://rstudio.gis.unbc.ca) and start by importing your libraries

library(tidyverse) library(maptools) library(sf) library(rgdal) library(leaflet)

Next, open the data from the L drive and plot it

guerry_shp <- readOGR(dsn = "~/L/GEOG413/lab07/guerry", layer = 'guerry') plot(guerry_shp)

The K-Means is simple to calculate in R, however it must first be converted into a dataframe

# Convert to data.frame guerry_df <- as.data.frame(guerry_shp) # Compute K-Means num_clusters <- 5 guerry_ke <- kmeans(guerry_df[, c('Crm_prs', 'Crm_prp', 'Litercy', 'Donatns', 'Infants', 'Suicids')], num_clusters)

Currently, we have a list of clusters, but to see what region belongs to each cluster we must bind the computed clusters to our original data.

# Bind cluster to input data !!! Do Not Sort Clusters !!! guerry_clusters <- cbind(guerry_shp, guerry_ke$cluster) # Convert to sf for leaflet guerry_sf <- st_as_sf(guerry_clusters) # Rename cluster column (index is 1 larger than number of input columns) colnames(guerry_sf)[ncol(guerry_shp)+1] <- "cluster"

Finally look at the results on the map

# Make colour palate cols = rainbow(num_clusters) # Show Leaflet map leaflet(guerry_sf %>% st_transform(4326)) %>% # Leaflet does not understand UTM! addProviderTiles("CartoDB.Positron") %>% # Add basemap addPolygons(color = cols, fillOpacity = 1) # Show different clusters as different colours

One more tip, as the clusters start in random locations, the output may vary slightly based on the starting state, if you need consistent output for testing purposes you can ‘seed’ the random number generator, before calling kmeans()

set.seed(42)

## Quadrat Analysis

Quadrat Analysis allows for exploratory searching of patterns in * locations of points*. In today’s lab we will be looking at the location of crashes on the lower mainland, using the spatStat package. *

*Note that spatStat does not natively work with*

*sf or ggplot, so our code will look a little different than previous labs.*

Like always we start by loading some data, this is a three-step process first we open the data, then convert the coordinates into locations, finally converting the sf data into spatial point patterns (ppp).

# Start of Lab crash_data <- read.csv2("~/L/GEOG413/lab07/Lower_mainland_Crashes_expanded.csv") %>% select(Latitude, Longitude, Location) # Convert spreadsheet to sf objects, note data is in WGS84 CRS crash_locations <- st_as_sf(crash_data, coords = c("Longitude", "Latitude"), crs = 4326) %>% st_transform(32610) # SpatStat requires data as ppp, transform to UTM 10N. Do not # Calculate statistics on unprojected data! crash_points <- as.ppp(crash_locations %>% st_transform(32610))

Examine the data by plotting

# Minimize plot margins to help see points, set default point to red # Strange bug where col in plot() changes background not point colour so moved to par() par(plt=c(0, 1, 0, 1), bg = 'white', col = alpha('red', 0.5)) # Plot points plot(crash_points, pch=20, cex = 0.5, main='q-plot', legend = F, label_value = F)

Next, to perform quadrat analysis the data needs to be split into quadrats, and the points inside counted, note that the size of your quadrats will affect the answer, we will talk more about this; for now to make the plots easy to see we will work with a 5×5 grid of quadrats.

# Count quadrats Q <- quadratcount(crash_points, nx= 5, ny=5) # Fade points to make counts readable rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col=alpha('white', 0.7)) # Add quadrat grid plot(Q, main='q-plot', col = alpha('black', 1), add=T)

Once all the points have been counted the next step is the chi-squared test, done with the quadrat.test() function.

# Chi-squared test qtest <- quadrat.test(Q) # Print results to console qtest # Optionally plot quadrats (count, expected count, deviation) plot(qtest, cex = .7, col = 'Black')

Lets talk about the output

Chi-squared test of CSR using quadrat counts data: X2 = 140895, df = 24, p-value < 2.2e-16 alternative hypothesis: two.sided Quadrats: 5 by 5 grid of tiles

The interpretation of this is that Chi-squared (X2) is larger than the critical value (df), suggesting that this data is clustered; as the p-value is less than 0.05 we can take these results to be statistically significant.

### Data Prep matters

You may have noticed the data we loaded was called Lower_mainland_Crashes_expanded; what does this mean? On the original data each point contained a count of crashes at that location, the expanded dataset contains extra points for each extra crash. If you try the above again changing to the original data how do your results differ?

crash_data <- read.csv2("~/L/GEOG413/lab07/Lower_mainland_Crashes_data.csv", sep = ',')

The original data has a smaller Chi-squared, what does that tell us about the data?

## Looking Through the Window

Another interesting point in this data is that there is a lot of ocean, and cars generally do not crash in the ocean. We can account for this by applying a window to the quadrat analysis.

# Open lowermainland shp as window outline <- readOGR(dsn = "~/L/GEOG413/lab07/lower_mainland", layer = 'lower_mainland') crash_data <- read.csv2("~/L/GEOG413/lab07/Lower_mainland_Crashes_expanded.csv") %>% select(Latitude, Longitude, Location) # Convert spreadsheet to sf objects, note data is in WGS84 CRS crash_locations <- st_as_sf(crash_data, coords = c("Longitude", "Latitude"), crs = 4326) %>% st_transform(32610) # SpatStat requires data as ppp, transform to UTM 10N. crash_points <- as.ppp(crash_locations %>% st_transform(32610)) # Set window on crash points window <- as.owin(outline) crash_points_window <- crash_points[window] # Count quadrats Q_window <- quadratcount(crash_points_window, nx=5, ny=5) # Minimize plot margins to help see points, set default point to red par(plt=c(0, 1, 0, 1), bg = 'white', col = alpha('red', 0.5)) # Plot points plot(crash_points_window, pch=20, cex = 0.5, main='q_window-plot', legend = F, label_value = F) # Fade points to make counts readable rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col=alpha('white', 0.7)) # Add quadrat grid plot(Q_window, main='q_window-plot', col = alpha('black', 1), add=T) # Chi-squared test q_window_test <- quadrat.test(Q_window) q_window_test plot(q_window_test, cex = 0.7, col = 'black')

In this case, you can see that the data is appearing less clustered and that the expected counts are not the same for all quadrats as they are scaled based on area.

### Evenly Spaced Data

regular_points <- readOGR(dsn = "~/L/GEOG413/lab07/synthetic_data", layer = 'regular_points') regular_locations <- st_as_sf(regular_points) # SpatStat requires data as ppp, transform to UTM 10N. crash_points_regular <- as.ppp(regular_locations %>% st_transform(32610)) # Open lowermainland shp as window outline <- readOGR(dsn = "~/L/GEOG413/lab07/lower_mainland", layer = 'lower_mainland') plot(outline) window <- as.owin(outline) # Set window on crash points crash_points_regular_window <- crash_points_regular[window] # Count quadrats Q_regular_window <- quadratcount(crash_points_regular_window, nx= 5, ny=5) # Minimize plot margins to help see points, set default point to red par(plt=c(0, 1, 0, 1), bg = 'white', col = alpha('red', 0.5)) # Plot points plot(crash_points_regular_window, pch=20, cex = 0.5, main='q_window_regular-plot', legend = F, label_value = F) # Fade points to make counts readable rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col=alpha('white', 0.7)) # Add quadrat grid plot(Q_regular_window, main='q_window_regular-plot', col = alpha('black', 1), add=T) q_regular_window_test <- quadrat.test(Q_regular_window) q_regular_window_test plot(q_regular_window_test, cex = 0.5, col = alpha('black', 1))

Chi-squared test of CSR using quadrat counts data: X2 = 7.8595, df = 21, p-value = 0.008577 alternative hypothesis: two.sided Quadrats: 22 tiles (irregular windows)

### Randomly Distributed Data

random_points <- readOGR(dsn = "~/L/GEOG413/lab07/synthetic_data", layer = 'random_points') random_locations <- st_as_sf(random_points) # SpatStat requires data as ppp, transform to UTM 10N. crash_points_random <- as.ppp(random_locations %>% st_transform(32610)) # Open lowermainland shp as window outline <- readOGR(dsn = "~/L/GEOG413/lab07/lower_mainland", layer = 'lower_mainland') plot(outline) window <- as.owin(outline) # Set window on crash points crash_points_random_window <- crash_points_random[window] # Count quadrats Q_random_window <- quadratcount(crash_points_random_window, nx= 5, ny=5) # Minimize plot margins to help see points, set default point to red par(plt=c(0, 1, 0, 1), bg = 'white', col = alpha('red', 0.5)) # Plot points plot(crash_points_random_window, pch=20, cex = 0.5, main='q_window_random-plot', legend = F, label_value = F) # Fade points to make counts readable rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col=alpha('white', 0.7)) # Add quadrat grid plot(Q_random_window, main='q_window_random-plot', col = alpha('black', 1), add=T) q_random_window_test <- quadrat.test(Q_random_window) q_random_window_test plot(q_random_window_test, cex = 0.7, col = alpha('black', 1))

Chi-squared test of CSR using quadrat counts data: X2 = 23.694, df = 21, p-value = 0.6163 alternative hypothesis: two.sided Quadrats: 22 tiles (irregular windows)

## Assignment

Answer the 3 questions below adding the missing lines to the code.

Submit your completed R script and screenshots of plots produced (or right click on graphs and save as to download to your computer). Make sure to caption the plots in a word document as the labels have been hidden to give more space.

library(tidyverse) library(maptools) library(sf) library(sp) library(rgdal) library(leaflet) library(spatstat) library(dplyr) library(rgeos) ################################################ ### Question 1: Which city in the layer 5CCITIES ### is the most unique based upon Marital Status ### NEVERMARRY, MARRIED, SEPARATED, WIDOWED, DIVORCED ### Use 5 clusters in K-Means ################################################ # Open data five_ccities <- readOGR(dsn = '~/L/GEOG413/lab07/assignment', layer = '5CCITIES') # Convert to DataFrame # K-Means Cluster # Bind clusters, and name column colnames(city_clusters@data)[length(five_ccities@data) + 1] <- 'clusters' # Convert to spatial dataframe for leaflet compatability city_sf <- st_as_sf(city_clusters) # Create colour pallet for leaflet pal <- colorFactor(rainbow(nrow(city_sf)), domain = city_sf$clusters) # Show Leaflet map # Display how many points in each cluster city_sf %>% group_by(clusters) %>% summarize(num_points = n()) ################################################ ### Question 2: Which two cities in the layer 5CCITIES ### is the most unique based upon housing status ### VACANT, OWNER_OCC, RENTER_OCC ### Normalize columns by HSE_UNITS ### Use 5 clusters in K-Means ################################################ # Open data five_ccities <- readOGR(dsn = '~/L/GEOG413/lab07/assignment', layer = '5CCITIES') # Convert to DataFrame # K-Means Cluster # Bind clusters, and name column city_clusters <- cbind(five_ccities, klusters$cluster) colnames(city_clusters@data)[length(five_ccities@data) + 1] <- 'clusters' # Convert to spatial dataframe for leaflet compatability city_sf <- st_as_sf(city_clusters) # Create colour pallet for leaflet pal <- colorFactor(rainbow(nrow(city_sf)), domain = city_sf$clusters) # Show Leaflet map # Display how many points in each cluster city_sf %>% group_by(clusters) %>% summarize(num_points = n()) ################################################ ### 3 Pts ### Question 3: emap, and dipin layers both contain ### information on water quality in lakes; using ### quadrat analysis determine how these layers ### compare in terms of point distribution ################################################ # Open Data dipin <- readOGR(dsn = '~/L/GEOG413/lab07/assignment', layer = 'dipin') emap <- readOGR(dsn = '~/L/GEOG413/lab07/assignment', layer = 'emap') # Shapfiles do not have CRS defined, define to allow transformation proj4string(dipin) <- proj4string(emap) <- # Transform to psuedo-webmercator projection dipin_proj <- emap_proj <- # Create spatial point patterns # Create window of North East States ne_usa <- readOGR(dsn = '~/L/GEOG413/lab07/assignment', layer = 'ne-usa') proj4string(ne_usa) <- CRS("+init=epsg:4326") ne_usa_proj <- spTransform(ne_usa, CRS("+init=epsg:3857")) outline <- gUnaryUnion(ne_usa_proj) window <- as.owin(outline) # Set window on ppp data # Plot points par(plt=c(0, 1, 0, 1), bg = 'white', col = 'red') plot(dipin_points_window, pch=20, cex = 0.5, main='q-plot', legend = F, label_value = F, which.marks=1) par(col='blue') plot(emap_points_window, pch=20, cex = 0.5, main='q-plot', legend = F, label_value = F, which.marks=1, add=T) # Count quadrats # Fade points to make counts readable rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col=alpha('white', 0.7)) # Add quadrat grid plot(dipin_qc, main='q-plot', col = alpha('black', 1), add=T) # Chi-squared test # Print results to console # Optionally plot quadrats (count, expected count, deviation) plot(dipin_chi2, cex = .5, col = 'Black') #### Both dipin and emap have same critical value, both are clustered, and have #### good p-values, however we can say that dipin is more clustered than e-map #### an interesting result given the obvious cluster in emap in North West region