Introduction
The objective of this project is to introduce the concept of Voronoi polygons and demonstrate their practical application in designing an air quality monitoring network for California. This study shows how to create Voronoi polygons from ozone measurement data, identify suitable monitoring station locations, and explore the significance of such a network in air quality management. You will also compare this with spatial demographic data from the CDC/ATSDR Social Vulnerability Index.
Libraries
Feel free to add any additional libraries you need
library(tmap)
library(raster)
library(gstat)
library(sf)
library(tidyverse)
library(tigris)
library(plotly)
library(maps)
library(terra)
library(spdep)The project uses the Tigris package to read in the census tract and county borders for California and SVI data to run code.
Census tracts are a spatial unit used to collect and analyse census data as a part ACS(american community survey). Census tract are best spatial unit to analyse community level data.These geographies will give you the flexibility to take a deeper look at any small area of a city or county1
AEA conic projection is used for this study. AEA conic projection preserve the spatial area where as UTM preserves the shape of the spatial entity. Therefore AEA is used for thematic maps where accurate area relationships is crucial, such as in demographic or statistical representations. Where as UTM is concerned with shape and direction and it is used for GIS applications, way finding maps etc.,The data for this lab is downloaded from CDC(center for disease control and prevention)2
SOVI is combined metric which is used by the US gov to measure the social vulnerability of communities in the US. SOVI is used for the following purposes:
• Assess community need during emergency preparedness planning.
• Estimate the type and amount of needed supplies such as food, water,
medicine, and bedding.
• Decide how many emergency personnel are required to assist
people.
• Identify areas in need of emergency shelters.
• Create a plan to evacuate people, accounting for those who have
special needs, such as those without vehicles, the elderly, or people
who do not speak English well.
• Identify communities that will need continued support to recover
following an emergency or natural disaster.
SOVI is a combination 4 major themes, which are explained below.
• Socioeconomic Status - RPL_THEME1
This consists of social economic variables such as poverty,
unemployment, education etc., • Household Characteristics -
RPL_THEME2
This consists of percentage of vulnerable population, household
relationships etc., • Racial & Ethnic Minority Status -
RPL_THEME3
This is about the racial and ethnic information •
Housing Type & Transportation - RPL_THEME4
This is about housing typology, density and vehicle ownership
etc.,
map1 <- qtm(tract.sovi.sf, fill = "RPL_THEME1", fill.palette = "Blues")
map2 <- qtm(tract.sovi.sf, fill = "RPL_THEME2", fill.palette = "YlOrBr")
map3 <- qtm(tract.sovi.sf, fill = "RPL_THEME3", fill.palette = "Greens")
map4 <- qtm(tract.sovi.sf, fill = "RPL_THEME4", fill.palette = "Reds")
tmap_arrange(map1, map2, map3, map4)
For the maps, higher value indicate higher vulnerability. There are
several patterns that is evident in the exploratory mapping. In general
we can identify high density of vulnerability in the central regions of
California(which is near Bakersfield, California). For instance, the
intensity of vulnerable population (Theme 2) can lead to lesser
employment rate (in Theme 1). We can identify visual similarities
between the mapping of Theme 1 and Theme 2.
To understand the struggle, we have to identify the sensitive
population. Which are:
Children (under 18), including teenagers, because their lungs are still developing and they breathe more air per pound of body weight than adults
Older adults (65+)
Other adults (18-64) with lung disease, such as asthma, or cardiovascular disease
People who are active outdoors, including outdoor workers
Majority of this is covered in theme 2 mapping.
map2
Therefore the areas in dark red is likely to struggle primarily with air pollution. However there could be secondary stressors as well. For instance, once you caught asthma, income level plays important role in mitigating the struggle. Person requires financial stability to fight asthma and lack of it can exacerbate the condition.
MAUP
The size of the spatial unit can significantly affect our inferences. For instance, lets look at county level mapping and census tract level mapping side by side:
map7 <- qtm(tract.sovi.sf, fill = "RPL_THEME1", fill.palette = "Blues")
map8 <- qtm(county.sovi.sf, fill = "RPL_THEME1", fill.palette = "Blues")
tmap_arrange(map7, map8)
As we can observe in the county-level mapping, the entire lower right portion of California is depicted as vulnerable. However, upon closer examination at the tract level mapping, it becomes evident that this is not the case. This discrepancy leads to the risk of committing an environmental fallacy, wherein conclusions drawn from group-level analysis are inaccurately assumed to hold true for individuals within that group. Hence, it is crucial to utilize higher definition data for increased accuracy.
Ozone data
We will now read in some data on Ozone. You can read about it in the lab book. Summarise what the data is showing.
california_map <- map("state", regions="california", exact=TRUE, plot=TRUE)
california_map$range## [1] -124.38342 -114.13319 32.53827 42.02073point.ozone <- read.csv("CA_OzonePopulation.csv")# Read in Ozone Data
point.ozone.sf <- st_as_sf(point.ozone, coords=c("LONGITUDE","LATITUDE"),crs=4326)
# Transform to projection 3310 - hint look at my set-up code chunk.
# I'm simply overwriting my sf data as my answer
point.ozone.sf <- transform(point.ozone.sf,3310)Qu. 3C. Ozone summary
Use R to summarise the ozone data, then answer these Qu.s in the text
- How many monitoring stations are there in the dataset?
summary(point.ozone)## LOCATION SITE_NAME SHORT_NAME LATITUDE
## Min. :2001 Length:451 Length:451 Min. :32.35
## 1st Qu.:2402 Class :character Class :character 1st Qu.:34.15
## Median :2844 Mode :character Mode :character Median :36.01
## Mean :2802 Mean :36.14
## 3rd Qu.:3134 3rd Qu.:37.94
## Max. :3759 Max. :41.85
## LONGITUDE OZONE_1000PPB POPULATION_DENSITY
## Min. :-124.2 Min. : 3.457 Min. : 0.0000
## 1st Qu.:-121.5 1st Qu.:23.617 1st Qu.: 0.5499
## Median :-120.0 Median :28.304 Median : 14.6029
## Mean :-119.7 Mean :30.347 Mean : 34.9754
## 3rd Qu.:-118.0 3rd Qu.:35.254 3rd Qu.: 53.2731
## Max. :-114.6 Max. :84.655 Max. :406.6252There are 451 monitoring stations
- What is the Ozone at site name: Beverly Hills-Franklin? Remember units!
Ozone_BH <- point.ozone$OZONE_1000PPB[point.ozone$SITE_NAME == "Beverly Hills-Franklin Canyon"]
print(Ozone_BH)## [1] 33.50973The ozone level os 33.5X10^3 parts per billion
- What is the maximum population density?
max(point.ozone$POPULATION_DENSITY)## [1] 406.6252QOzone vs Population density.
library(plotly)
fig <- plot_ly(data = point.ozone, x = point.ozone$OZONE_1000PPB, y = point.ozone$POPULATION_DENSITY)
figI expected to see ozone increase with increase in population, as the VOC from human activity is crucial in the creation of Ozone at this atmospheric level. However it could be due to the fact that the number of ozone measuring devices are less in highly dense areas.Most ozone reading are seen at 0 population areas.
Ozone maps
#population density
tmap_mode('view')
qtm(point.ozone.sf,
dots.col='POPULATION_DENSITY')#population density
tmap_mode('view')
qtm(point.ozone.sf,
dots.col='OZONE_1000PPB')As we can see here higher levels of ozone is correlated with higher elevation. This is because ozone is mainly detected at 20-25 km high from the sea level. The ozone values are really influenced by elevation rather than population density.
# DO NOT TOUCH!
rm(list=ls())
# THIS WON'T READ IN UNLESS YOU HAVE A PROJECT RUNNING AND YOUR DATA IN YOUR PROJECT FOLDER.
# SEE LAB 2 IF YOU ARE CONFUSED
ozone.sf <- st_read("point.ozone.AEA.geojson")## Reading layer `point.ozone.AEA' from data source
## `C:\Users\arjun\OneDrive - The Pennsylvania State University\Desktop\point.ozone.AEA.geojson'
## using driver `GeoJSON'
## Simple feature collection with 451 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -353684.6 ymin: -624323.6 xmax: 501010.6 ymax: 432274.9
## Projected CRS: NAD83 / California Alberscounty.data.sf <- st_read("county.AllData.AEA.geojson")## Reading layer `county.AllData.AEA' from data source
## `C:\Users\arjun\OneDrive - The Pennsylvania State University\Desktop\county.AllData.AEA.geojson'
## using driver `GeoJSON'
## Simple feature collection with 58 features and 38 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -373976.1 ymin: -604478.4 xmax: 540037.5 ymax: 450022.5
## Projected CRS: NAD83 / California Albersstate.border.sf <- states(cb = TRUE) %>%
select(NAME, geometry) %>%
filter(NAME == "California") %>%
st_transform(3310)## | | | 0% | |= | 1% | |== | 2% | |=== | 4% | |==== | 6% | |====== | 9% | |======== | 12% | |========== | 14% | |========== | 15% | |============ | 17% | |============== | 20% | |=============== | 21% | |================ | 23% | |================ | 24% | |================== | 26% | |=================== | 27% | |==================== | 29% | |====================== | 31% | |======================== | 34% | |========================== | 37% | |=========================== | 39% | |============================== | 42% | |=============================== | 45% | |================================ | 45% | |================================= | 48% | |=================================== | 51% | |===================================== | 53% | |======================================= | 56% | |========================================= | 59% | |=========================================== | 61% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |=================================================== | 72% | |===================================================== | 75% | |====================================================== | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================== | 83% | |============================================================ | 85% | |============================================================== | 88% | |================================================================ | 91% | |================================================================= | 93% | |=================================================================== | 96% | |===================================================================== | 99% | |======================================================================| 100%county.border.sf <- counties("CA", cb = TRUE) %>%
select(GEOID, geometry) %>%
st_transform(3310)## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%state.border.terra <- vect(state.border.sf)# you should be able to simply run this code
ozone.terra <- vect(ozone.sf)
ozone.voronoi.sf <- voronoi(ozone.terra)
ozone.voronoi.sf <- st_as_sf(crop(ozone.voronoi.sf,state.border.terra))
tm_shape(ozone.voronoi.sf) +
tm_polygons("OZONE_1000PPB",palette="YlGnBu")+
tm_legend(position = c("right", "top"))+
tm_layout(main.title="Voroni Tesselation of Ozone in CA (1000 PPB)",
main.title.size=.8,main.title.fontface=2)tm_shape(county.data.sf) +
tm_polygons("OZONE_1000PPB",palette="YlGnBu",breaks=seq(0,100,by=20))+
tm_legend(position = c("right", "top"))+
tm_layout(main.title="County Averages of Ozone in CA (1000 PPB)",
main.title.size=.8,main.title.fontface=2)tm_shape(county.data.sf) +
tm_borders()+
tm_shape(ozone.sf) +
tm_dots(col="OZONE_1000PPB",palette="YlGnBu",size=.6,alpha=.6)+
tm_legend(position = c("right", "top"))+
tm_layout(main.title="Point values of Ozone in CA (1000 PPB)",
main.title.size=.8,main.title.fontface=2)Conclusion
Map1-This map is a representation of the tessellation operation that is
done based on the ozone weather station location. This helps to identify
what is the zone of influence for each point. By dividing the space it
makes sure there are no two points in the same polygon. After the
process, each polygon is color-coded based on its attribute, that is
ozone level. An advantage is that it produces fairly accurate results as
divides the areas based on the availability of information. Since we are
dividing based on the “equal access” idea, other biases in area division
is reduced. A advantage could be disadvantage
Map2-In this map area is divided based on county boundary line. Each polygon is cropped by the county and average value is calculated to form the color-code of the census tract. An advantage of doing this is it help to compare with other social data since socio-ecomic information is universally divided based on census boundaries. This can help inform policy decisions. Disadvantages include, loss of detail, by averaging we are smoothing the data, this leads to loss of information; it may lead to ecological fallacy(generalizations).
Maps3-In this map area is divided based circle of influence. The circles are color coded based on the attribute of the points. An advantage of the this method (assuming the radius of the circle is theoretically defined) is that it give accurate information inside the circle. It also shows the areas that is not covered by any stations and this inform the deployment of new station here. This could also be a disadvantage, because certain area are missed in this method.
The division of area and choice of the map will depend on the purpose of the mapping. If its used for socio-economic research with census counties as spatial unit can help, where as it is for research regarding natural systems, such weather analysis, it is better to go with voronoi polygons. Finally the point mapping can be best if the scale of study is small, for instance, you want to study the ozone level in a urban junction or plaza.
https://www.census.gov/data/academy/data-gems/2018/tract.html↩︎
Centers for Disease Control and Prevention/ Agency for Toxic Substances and Disease Registry/ Geospatial Research, Analysis, and Services Program. CDC/ATSDR Social Vulnerability Index [2020] Database [California]. https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html. Accessed on [11/7/23].↩︎