+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+GGIS 224 Final Project Data Wrangling and +Spatial Operations
+Warren Jodjana
+December 14th 2023
+ +
+
+Environment Setup
+library(raster)
+library(dplyr)
+library(sf)
+library(terra)
+library(tmap)
+library(readxl)
+library(data.table)
+library(spatstat)
+
+Clean & Wrangle Data
+
+
+Air Quality Index (AQI)
+
+Load AQI data using
+
+Load AQI data using read_excel
+airquality <- read_excel("air-quality-index.xlsx")
+glimpse(airquality)## Rows: 801
+## Columns: 4
+## $ geoid <dbl> 17031291600, 17031580200, 17031590500, 17031590700, 17031610…
+## $ CMAQ_NO2 <dbl> 19.05938, 20.31526, 19.64178, 17.46373, 18.57812, 16.14407, …
+## $ CMAQ_O3 <dbl> 28.53719, 27.68554, 28.14529, 29.99657, 29.17645, 30.71135, …
+## $ CMAQ_PM25 <dbl> 11.103176, 11.435322, 11.288680, 10.704017, 10.921030, 10.60…
+Load tract data using
+
+Load tract data using st_read
+airquality_tracts <- st_read("tracts-aqi.geojson")## Reading layer `tracts-aqi' from data source
+## `/Users/warr/Documents/GitHub/ggis224-finalproject/tracts-aqi.geojson'
+## using driver `GeoJSON'
+## Simple feature collection with 801 features and 9 fields
+## Geometry type: MULTIPOLYGON
+## Dimension: XY
+## Bounding box: xmin: -87.94025 ymin: 41.64429 xmax: -87.52366 ymax: 42.02392
+## Geodetic CRS: WGS 84data.table(airquality_tracts)## statefp10 name10 commarea_n namelsad10 commarea geoid10 notes
+## 1: 17 8424 44 Census Tract 8424 44 17031842400
+## 2: 17 8403 59 Census Tract 8403 59 17031840300
+## 3: 17 8411 34 Census Tract 8411 34 17031841100
+## 4: 17 8412 31 Census Tract 8412 31 17031841200
+## 5: 17 8390 32 Census Tract 8390 32 17031839000
+## ---
+## 797: 17 704 7 Census Tract 704 7 17031070400
+## 798: 17 705 7 Census Tract 705 7 17031070500
+## 799: 17 1303 13 Census Tract 1303 13 17031130300
+## 800: 17 2922 29 Census Tract 2922 29 17031292200
+## 801: 17 6309 63 Census Tract 6309 63 17031630900
+## tractce10 countyfp10 geometry
+## 1: 842400 031 <XY[1]>
+## 2: 840300 031 <XY[1]>
+## 3: 841100 031 <XY[1]>
+## 4: 841200 031 <XY[1]>
+## 5: 839000 031 <XY[1]>
+## ---
+## 797: 070400 031 <XY[1]>
+## 798: 070500 031 <XY[1]>
+## 799: 130300 031 <XY[1]>
+## 800: 292200 031 <XY[1]>
+## 801: 630900 031 <XY[1]>
+
+Merge data using geoid
+airquality <- merge(airquality_tracts, airquality, by.x="geoid10", by.y="geoid", all.x=TRUE)
+
+Subset data
+airquality <- airquality %>%
+ select(-c("statefp10", "commarea_n", "commarea", "notes", "countyfp10"))
+Inspect data using
+
+Inspect data using data.table
+data.table(airquality)## geoid10 name10 namelsad10 tractce10 CMAQ_NO2 CMAQ_O3
+## 1: 17031010100 101 Census Tract 101 010100 11.07882 35.63879
+## 2: 17031010201 102.01 Census Tract 102.01 010201 11.62339 35.07472
+## 3: 17031010202 102.02 Census Tract 102.02 010202 11.23748 35.50884
+## 4: 17031010300 103 Census Tract 103 010300 10.96102 36.05563
+## 5: 17031010400 104 Census Tract 104 010400 11.43238 35.70695
+## ---
+## 797: 17031843700 8437 Census Tract 8437 843700 19.47959 28.76974
+## 798: 17031843800 8438 Census Tract 8438 843800 16.14407 30.71135
+## 799: 17031843900 8439 Census Tract 8439 843900 11.91338 34.78665
+## 800: 17031980000 9800 Census Tract 9800 980000 19.21362 28.57602
+## 801: 17031980100 9801 Census Tract 9801 980100 16.04343 30.61503
+## CMAQ_PM25 geometry
+## 1: 8.636906 <XY[1]>
+## 2: 8.948587 <XY[1]>
+## 3: 8.700630 <XY[1]>
+## 4: 8.518667 <XY[1]>
+## 5: 8.615976 <XY[1]>
+## ---
+## 797: 10.971199 <XY[1]>
+## 798: 10.605793 <XY[1]>
+## 799: 8.794460 <XY[1]>
+## 800: 10.323887 <XY[1]>
+## 801: 11.711032 <XY[1]>
+
+Save cleaned data
+write.csv(airquality, "airquality.csv")
+
+Traffic & Vehicular Emissions
+
+Load traffic data using
+
+Load traffic data using fread
+traffic <- fread("traffic-old.csv")
+glimpse(traffic)## Rows: 1,279
+## Columns: 9
+## $ ID <int> 414, 176, 1367, 316, 129…
+## $ `Traffic Volume Count Location Address` <chr> "5838 West", "320 East",…
+## $ Street <chr> "Lake St", "76th St", "5…
+## $ `Date of Count` <chr> "11/14/2006", "03/28/200…
+## $ `Total Passing Vehicle Volume` <int> 7100, 8600, 53500, 700, …
+## $ `Vehicle Volume By Each Direction of Traffic` <chr> "East Bound: 3600 / West…
+## $ Latitude <dbl> 41.88790, 41.75654, 41.7…
+## $ Longitude <dbl> -87.77106, -87.61733, -8…
+## $ Location <chr> "(41.887904, -87.771064)…
+
+Convert to spatial format
+traffic <- st_as_sf(traffic, coords = c("Longitude", "Latitude"), crs = 4326)
+Inspect data using
+
+Inspect data using data.table
+data.table(traffic)## ID Traffic Volume Count Location Address Street Date of Count
+## 1: 414 5838 West Lake St 11/14/2006
+## 2: 176 320 East 76th St 03/28/2006
+## 3: 1367 1730 East 57th Dr 08/24/2006
+## 4: 316 125 East 24th St 03/30/2006
+## 5: 1294 2924 East 130th St 08/29/2006
+## ---
+## 1275: 257 256 East 35th St 03/29/2006
+## 1276: 950 2032 North Pulaski Rd 10/05/2006
+## 1277: 118 6020 South Halsted St 03/21/2006
+## 1278: 1091 3454 West Franklin St 10/19/2006
+## 1279: 1045 316 West 51st St 10/19/2006
+## Total Passing Vehicle Volume Vehicle Volume By Each Direction of Traffic
+## 1: 7100 East Bound: 3600 / West Bound: 3500
+## 2: 8600 East Bound: 3800 / West Bound: 4800
+## 3: 53500 East Bound: 27800 / West Bound: 25700
+## 4: 700 East Bound: 400 / West Bound: 300
+## 5: 4200 East Bound: 2300 / West Bound: 1900
+## ---
+## 1275: 10200 East Bound: 4800 / West Bound: 5400
+## 1276: 21300 North Bound: 10800 / South Bound: 10500
+## 1277: 13200 North Bound: 6900 / South Bound: 6300
+## 1278: 5000 East Bound: 2200 / West Bound: 2800
+## 1279: 14400 East Bound: 8100 / West Bound: 6300
+## Location geometry
+## 1: (41.887904, -87.771064) POINT (-87.77106 41.8879)
+## 2: (41.756542, -87.617335) POINT (-87.61733 41.75654)
+## 3: (41.792663, -87.582231) POINT (-87.58223 41.79266)
+## 4: (41.849302, -87.622658) POINT (-87.62266 41.8493)
+## 5: (41.659177, -87.552112) POINT (-87.55211 41.65918)
+## ---
+## 1275: (41.831081, -87.619334) POINT (-87.61933 41.83108)
+## 1276: (41.918169, -87.72659) POINT (-87.72659 41.91817)
+## 1277: (41.78453, -87.644905) POINT (-87.6449 41.78453)
+## 1278: (41.890003, -87.713444) POINT (-87.71344 41.89)
+## 1279: (41.801796, -87.633948) POINT (-87.63395 41.8018)
+
+Save cleaned data
+write.csv(traffic, "traffic.csv")
+
+Electric Vehicles
+
+Load Alternative Fuel data using
+
+Load Alternative Fuel data using fread
+alternativefuel <- fread("alternativefuel.csv")
+
+Convert to spatial format
+alternativefuel <- st_as_sf(alternativefuel, coords = c("Longitude", "Latitude"), crs=4326)
+data.table(alternativefuel)## ID Fuel Type Code Station Name
+## 1: 33089 E85 City of Chicago FS#4
+## 2: 180625 ELEC PEOPLES GAS CENTRAL EV2
+## 3: 172397 ELEC IRVINE CO OFC 71SW B2 02-03
+## 4: 205433 ELEC Cityfront Place Parking Garage - Floor 2 Left
+## 5: 201403 ELEC DRAKE TOWER DRAKE 1
+## ---
+## 333: 163425 ELEC BUCK MANAGEMENT LEVEL 2 327
+## 334: 61218 ELEC PARKING GARAGE 111SW STATION 2
+## 335: 151753 ELEC Jewel Osco 4250 N Lincoln Chicago IL
+## 336: 114160 ELEC Millennium Lakeside Garage - Tesla Destination
+## 337: 62312 CNG Ozinga Energy - Lumber St
+## Street Address
+## 1: 3746 S Iron
+## 2: 4207 W 35th Pl Central Shop
+## 3: 71 S Wacker Dr
+## 4: 400 N. McClurg Court
+## 5: 172 E Walton Pl
+## ---
+## 333: 845 W Madison
+## 334: 111 S Wacker Dr
+## 335: 4250 N Lincoln
+## 336: 5 S Columbus Dr
+## 337: 2222 S Lumber St
+## Intersection Directions
+## 1: Iron and 37th
+## 2:
+## 3:
+## 4: Located on the 2nd floor of the parking garage on the south side of the elevators. This is the left charging station.
+## 5:
+## ---
+## 333:
+## 334:
+## 335:
+## 336:
+## 337: 2 Blocks East of Cermak and Halsted.
+## City State ZIP Plus4 Station Phone Status Code
+## 1: Chicago IL 60609 NA E
+## 2: Chicago IL 60646 NA 888-758-4389 T
+## 3: Chicago IL 60606 NA 888-758-4389 E
+## 4: Chicago IL 60611 NA 800-663-5633 T
+## 5: Chicago IL 60611 NA 888-758-4389 E
+## ---
+## 333: Chicago IL 60607 NA 888-758-4389 E
+## 334: Chicago IL 60606 NA 888-758-4389 E
+## 335: Chicago IL 60618 NA 888-264-2208 E
+## 336: Chicago IL 60603 NA 312-616-0600 877-798-3752 E
+## 337: Chicago IL 60616 NA 888-656-0141 E
+## Expected Date Groups With Access Code
+## 1: Private - Government only
+## 2: 12/09/2023 TEMPORARILY UNAVAILABLE (Public)
+## 3: Public
+## 4: 12/20/2022 TEMPORARILY UNAVAILABLE (Public)
+## 5: Public
+## ---
+## 333: Public
+## 334: Public
+## 335: Public
+## 336: Public
+## 337: Private
+## Access Days Time
+## 1:
+## 2: 24 hours daily
+## 3: 24 hours daily
+## 4: 24 hours daily
+## 5: 24 hours daily
+## ---
+## 333: 24 hours daily
+## 334: 24 hours daily
+## 335: Mon 6:00am - 11:59pm; Tue 6:00am - 11:59pm; Wed 6:00am - 11:59pm; Thu 6:00am - 11:59pm; Fri 6:00am - 11:59pm; Sat 6:00am - 11:59pm; Sun 6:00am - 11:59pm
+## 336:
+## 337:
+## Cards Accepted BD Blends NG Fill Type Code NG PSI EV Level1 EVSE Num
+## 1: NA NA
+## 2: NA NA
+## 3: NA NA
+## 4: NA NA
+## 5: NA NA
+## ---
+## 333: NA NA
+## 334: NA NA
+## 335: NA NA
+## 336: NA NA
+## 337: Q 3600 NA
+## EV Level2 EVSE Num EV DC Fast Count EV Other Info EV Network
+## 1: NA NA NA
+## 2: 2 NA NA ChargePoint Network
+## 3: 2 NA NA ChargePoint Network
+## 4: 1 NA NA Blink Network
+## 5: 2 NA NA ChargePoint Network
+## ---
+## 333: 2 NA NA ChargePoint Network
+## 334: 2 NA NA ChargePoint Network
+## 335: 2 NA NA Volta
+## 336: 11 NA NA Tesla Destination
+## 337: NA NA NA
+## EV Network Web Geocode Status
+## 1: 200-8
+## 2: http://www.chargepoint.com/ GPS
+## 3: http://www.chargepoint.com/ GPS
+## 4: http://www.blinkcharging.com/ GPS
+## 5: http://www.chargepoint.com/ GPS
+## ---
+## 333: http://www.chargepoint.com/ GPS
+## 334: http://www.chargepoint.com/ GPS
+## 335: https://voltacharging.com/ GPS
+## 336: https://www.tesla.com/destination-charging 200-8
+## 337: 200-8
+## Date Last Confirmed Updated At Owner Type Code
+## 1: 09/14/2023 09/14/2023 09:01:49 AM +0000 LG
+## 2: 12/13/2023 12/12/2023 06:38:23 PM +0000
+## 3: 12/13/2023 12/12/2023 06:13:10 PM +0000
+## 4: 12/13/2023 12/12/2023 06:08:15 PM +0000
+## 5: 12/13/2023 12/12/2023 06:54:18 PM +0000
+## ---
+## 333: 12/13/2023 12/12/2023 06:40:53 PM +0000
+## 334: 12/13/2023 12/12/2023 06:26:55 PM +0000
+## 335: 12/13/2023 12/12/2023 07:28:45 PM +0000
+## 336: 10/06/2022 01/18/2023 09:15:20 AM +0000 P
+## 337: 06/14/2022 12/06/2023 10:22:24 AM +0000 P
+## Federal Agency ID Federal Agency Name Open Date Hydrogen Status Link
+## 1: NA 10/15/2008 NA
+## 2: NA 01/27/2021 NA
+## 3: NA 01/27/2021 NA
+## 4: NA 12/22/2021 NA
+## 5: NA 11/19/2021 NA
+## ---
+## 333: NA 06/17/2020 NA
+## 334: NA 04/12/2014 NA
+## 335: NA 08/30/2019 NA
+## 336: NA 01/02/2019 NA
+## 337: NA 09/01/2012 NA
+## NG Vehicle Class LPG Primary E85 Blender Pump Location
+## 1: NA FALSE (41.826171, -87.659383)
+## 2: NA NA (41.989265, -87.730213)
+## 3: NA NA (41.880925, -87.636038)
+## 4: NA NA (41.8891899, -87.6181012)
+## 5: NA NA (41.900228, -87.622669)
+## ---
+## 333: NA NA (41.880655, -87.648861)
+## 334: NA NA (41.880164, -87.636483)
+## 335: NA NA (41.9587395, -87.682646)
+## 336: NA NA (41.882094, -87.620496)
+## 337: HD NA NA (41.852403, -87.642494)
+## geometry
+## 1: POINT (-87.65938 41.82617)
+## 2: POINT (-87.73021 41.98926)
+## 3: POINT (-87.63604 41.88092)
+## 4: POINT (-87.6181 41.88919)
+## 5: POINT (-87.62267 41.90023)
+## ---
+## 333: POINT (-87.64886 41.88065)
+## 334: POINT (-87.63648 41.88016)
+## 335: POINT (-87.68265 41.95874)
+## 336: POINT (-87.6205 41.88209)
+## 337: POINT (-87.64249 41.8524)
+
+Subset data
+electric_vehicles <- alternativefuel %>%
+ select(`ID`, `Fuel Type Code`, `Station Name`, `Street Address`, `Location`, `geometry`) %>%
+ filter(`Fuel Type Code` == "ELEC") %>%
+ mutate(`Street Address` = ifelse(`Station Name` == "Paul Simon Chicago JCC", "3348 S Kedzie Ave", `Street Address`)) %>%
+ select(-c("Fuel Type Code"))
+Inspect data using
+
+Inspect data using data.table
+data.table(electric_vehicles)## ID Station Name
+## 1: 180625 PEOPLES GAS CENTRAL EV2
+## 2: 172397 IRVINE CO OFC 71SW B2 02-03
+## 3: 205433 Cityfront Place Parking Garage - Floor 2 Left
+## 4: 201403 DRAKE TOWER DRAKE 1
+## 5: 196870 Paul Simon Chicago JCC
+## ---
+## 305: 173546 UNITED ORD LOT UNITED 5
+## 306: 163425 BUCK MANAGEMENT LEVEL 2 327
+## 307: 61218 PARKING GARAGE 111SW STATION 2
+## 308: 151753 Jewel Osco 4250 N Lincoln Chicago IL
+## 309: 114160 Millennium Lakeside Garage - Tesla Destination
+## Street Address Location
+## 1: 4207 W 35th Pl Central Shop (41.989265, -87.730213)
+## 2: 71 S Wacker Dr (41.880925, -87.636038)
+## 3: 400 N. McClurg Court (41.8891899, -87.6181012)
+## 4: 172 E Walton Pl (41.900228, -87.622669)
+## 5: 3348 S Kedzie Ave (41.878114, -87.629798)
+## ---
+## 305: 5600 Bessie Coleman Dr (41.98128, -87.888622)
+## 306: 845 W Madison (41.880655, -87.648861)
+## 307: 111 S Wacker Dr (41.880164, -87.636483)
+## 308: 4250 N Lincoln (41.9587395, -87.682646)
+## 309: 5 S Columbus Dr (41.882094, -87.620496)
+## geometry
+## 1: POINT (-87.73021 41.98926)
+## 2: POINT (-87.63604 41.88092)
+## 3: POINT (-87.6181 41.88919)
+## 4: POINT (-87.62267 41.90023)
+## 5: POINT (-87.6298 41.87811)
+## ---
+## 305: POINT (-87.88862 41.98128)
+## 306: POINT (-87.64886 41.88065)
+## 307: POINT (-87.63648 41.88016)
+## 308: POINT (-87.68265 41.95874)
+## 309: POINT (-87.6205 41.88209)
+
+Save cleaned data
+write.csv(electric_vehicles, "electric-vehicles.csv")
+
+
+
+
+Spatial Operations
+
+
+Traffic Kernel Density Estimation
+boundary <- st_read("boundaries.geojson")## Reading layer `boundaries' from data source
+## `/Users/warr/Documents/GitHub/ggis224-finalproject/boundaries.geojson'
+## using driver `GeoJSON'
+## Simple feature collection with 77 features and 9 fields
+## Geometry type: MULTIPOLYGON
+## Dimension: XY
+## Bounding box: xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
+## Geodetic CRS: WGS 84
+
+Ensure both points and boundaries are the same CRS
+boundary <- st_transform(boundary, crs = 4326)
+
+Creating a two-dimensional point-pattern object
+coords <- st_coordinates(traffic) # extract coordinates from sf object
+bbox <- st_bbox(boundary) # define the observation window (bounding box) using the boundary data
+window <- owin(xrange = bbox[c("xmin", "xmax")], yrange = bbox[c("ymin", "ymax")])
+traffic.ppp <- ppp(coords[,1], coords[,2], window = window) # Create a ppp (two-dimensional point-pattern) object
+
+Perform Kernel Density Estimation
+kde <- density(traffic.ppp) # perform kernel density estimation
+kde_raster <- raster(kde) # converting KDE back to raster form
+kde_raster_masked <- mask(kde_raster, boundary) # clip the KDE in raster form to chicago boundaries
+
+Check bandwith and type of kernel
+str(kde) ## List of 10
+## $ v : num [1:128, 1:128] 0.28 0.316 0.356 0.4 0.449 ...
+## $ dim : int [1:2] 128 128
+## $ xrange: num [1:2] -87.9 -87.5
+## $ yrange: num [1:2] 41.6 42
+## $ xstep : num 0.00325
+## $ ystep : num 0.00296
+## $ xcol : num [1:128] -87.9 -87.9 -87.9 -87.9 -87.9 ...
+## $ yrow : num [1:128] 41.6 41.6 41.7 41.7 41.7 ...
+## $ type : chr "real"
+## $ units :List of 3
+## ..$ singular : chr "unit"
+## ..$ plural : chr "units"
+## ..$ multiplier: num 1
+## ..- attr(*, "class")= chr "unitname"
+## - attr(*, "class")= chr "im"
+## - attr(*, "sigma")= num 0.0473
+## - attr(*, "kernel")= chr "gaussian"
+## - attr(*, "kerdata")=List of 5
+## ..$ sigma : num 0.0473
+## ..$ varcov : NULL
+## ..$ cutoff : num 0.378
+## ..$ warnings: NULL
+## ..$ kernel : chr "gaussian"
+
+Kernel Density Estimation on Traffic in Chicago
+tmap_mode("view")## tmap mode set to interactive viewingchikde <- tm_shape(boundary) + tm_borders() +
+ tm_shape(kde_raster_masked) +
+ tm_raster(title="Traffic Density", style="cont", palette="Blues") +
+ tm_scale_bar(position = c("left", "bottom"))
+
+chikde
+
+Proximity Analysis
+nearest_traffic_index <- st_nearest_feature(electric_vehicles, traffic)
+nearest_traffic_distances <- st_distance(electric_vehicles, traffic[nearest_traffic_index, ], by_element = TRUE)
+
+
+tmap_mode("view")## tmap mode set to interactive viewingtm <- tm_shape(traffic) +
+ tm_bubbles(size = "ID", col = "red", border.col = "black", alpha = 0.5) +
+ tm_shape(electric_vehicles) +
+ tm_symbols(col = "blue", size = 0.1) +
+ tm_layout(title = "EV Stations and Nearest Traffic Data Points")
+
+tm## Legend for symbol sizes not available in view mode.
+
+Comparative Analysis
+tmap_mode("plot")## tmap mode set to plottingairquality_summary <- airquality %>%
+ group_by(geoid10) %>%
+ summarize(mean_NO2 = mean(CMAQ_NO2, na.rm = TRUE),
+ mean_PM25 = mean(CMAQ_PM25, na.rm = TRUE),
+ mean_O3 = mean(CMAQ_O3, na.rm = TRUE))
+
+
+
+plot_NO2 <- tm_shape(airquality_summary) +
+ tm_polygons("mean_NO2", title = "Mean NO2") +
+ tm_layout(main.title = "Spatial Distribution of NO2", main.title.size = 1)
+
+plot_PM25 <- tm_shape(airquality_summary) +
+ tm_polygons("mean_PM25", title = "Mean PM2.5") +
+ tm_layout(main.title = "Spatial Distribution of PM2.5", main.title.size = 1)
+
+plot_O3 <- tm_shape(airquality_summary) +
+ tm_polygons("mean_O3", title = "Mean O3") +
+ tm_layout(main.title = "Spatial Distribution of O3", main.title.size = 1)
+
+
+tmap_arrange(plot_NO2, plot_PM25, plot_O3)
+
+
+
+
+
+