library(tidycensus)
library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(tmaptools)
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
library(ggplot2)
library(rasterVis)
## Warning: package 'rasterVis' was built under R version 4.4.2
## Loading required package: lattice
# 1. Make a map to show the spatial distribution of percentage of Hispanic population of the county you analyzed, and include compass and scale bar
census_api_key("80cd72e3cddeaf00e930b56289358c534d1ecfe6", install="TRUE", overwrite=TRUE)
## Your original .Renviron will be backed up and stored in your R HOME directory if needed.
## Your API key has been stored in your .Renviron and can be accessed by Sys.getenv("CENSUS_API_KEY").
## To use now, restart R or run `readRenviron("~/.Renviron")`
## [1] "80cd72e3cddeaf00e930b56289358c534d1ecfe6"
var <- c(hispanic='B03002_012E',
pop_total='B03002_012E')
st <- "TX"
ct <- "201"
harris_tracts <- get_acs(
geography = "tract",
variables = "B03002_012",
year = 2021,
state = "TX",
county = "201",
geometry = TRUE)
## Getting data from the 2017-2021 5-year ACS
## Downloading feature geometry from the Census website. To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
harris_tracts <- harris_tracts %>%
mutate(percent_hispanic = estimate / sum(estimate) * 100)
tm_shape(harris_tracts) +
tm_polygons(
col = "percent_hispanic",
palette = "YlOrRd",
title = "Hispanic %") +
tm_compass(
type = "8star",
position = c("right", "top"),
size = 1) +
tm_scale_bar(
position = c("left", "bottom")) +
tm_layout(
legend.position = c("center", "bottom"),
legend.outside = TRUE,
legend.outside.position = "right")
# 2. Find any point shapefile data within the county and make a map with census tract as the background
points_data <- st_read("C:/Users/cruzs/Downloads/College.shp")
## Reading layer `College' from data source `C:\Users\cruzs\Downloads\College.shp' using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2988500 ymin: 13767310 xmax: 3271050 ymax: 13995650
## Projected CRS: NAD83 / Texas South Central (ftUS)
# attempting to fix code error "The shape points_data is invalid (after reprojection). See sf::st_is_valid "
tmap_options(check.and.fix = TRUE)
points_data <- st_make_valid(points_data)
st_crs(harris_tracts)
## Coordinate Reference System:
## User input: NAD83
## wkt:
## GEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]]
st_crs(points_data)
## Coordinate Reference System:
## User input: NAD83 / Texas South Central (ftUS)
## wkt:
## PROJCRS["NAD83 / Texas South Central (ftUS)",
## BASEGEOGCRS["NAD83",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4269]],
## CONVERSION["SPCS83 Texas South Central zone (US survey foot)",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",27.8333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-99,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",30.2833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",28.3833333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",1968500,
## LENGTHUNIT["US survey foot",0.304800609601219],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",13123333.333,
## LENGTHUNIT["US survey foot",0.304800609601219],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["US survey foot",0.304800609601219]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["US survey foot",0.304800609601219]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["United States (USA) - Texas - counties of Aransas; Atascosa; Austin; Bandera; Bee; Bexar; Brazoria; Brewster; Caldwell; Calhoun; Chambers; Colorado; Comal; De Witt; Dimmit; Edwards; Fayette; Fort Bend; Frio; Galveston; Goliad; Gonzales; Guadalupe; Harris; Hays; Jackson; Jefferson; Karnes; Kendall; Kerr; Kinney; La Salle; Lavaca; Live Oak; Matagorda; Maverick; McMullen; Medina; Presidio; Real; Refugio; Terrell; Uvalde; Val Verde; Victoria; Waller; Wharton; Wilson; Zavala."],
## BBOX[27.78,-105,30.67,-93.76]],
## ID["EPSG",2278]]
points_data <- st_transform(points_data, st_crs(harris_tracts))
tm_shape(harris_tracts) +
tm_borders() +
tm_shape(points_data) +
tm_dots(
col = "royalblue",
size = 0.3,
title = "Colleges Point Data") +
tm_compass(
type = "8star",
position = c("right", "top"),
size = 1) +
tm_scale_bar(
position = c("left", "bottom")) +
tm_layout(
legend.position = c("center", "bottom"),
legend.outside = TRUE,
legend.outside.position = "right")
# 3. Find any line shapefile data within the county and make a map with census tract as the background
lines_data <- st_read("C:/Users/cruzs/Downloads/hwy.shp")
## Reading layer `HWY' from data source `C:\Users\cruzs\Downloads\HWY.shp' using driver `ESRI Shapefile'
## Simple feature collection with 84 features and 3 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 2933487 ymin: 13747360 xmax: 3273564 ymax: 13992790
## Projected CRS: NAD83 / Texas South Central (ftUS)
tm_shape(harris_tracts) +
tm_borders() +
tm_shape(lines_data) +
tm_lines(
col = "blue",
lwd = 1) +
tm_add_legend(
type = "line",
col = "blue",
title = "Highway Line Data") +
tm_compass(
type = "8star",
position = c("right", "top"),
size = 1) +
tm_scale_bar(
position = c("left", "bottom")) +
tm_layout(
legend.position = c("center", "bottom"),
legend.outside = TRUE,
legend.outside.position = "right")
# 4. Find any polygon shapefile data within the county and make a map with census tract as the background
polygon_data <- st_read("C:/Users/cruzs/Downloads/City.shp")
## Reading layer `City' from data source `C:\Users\cruzs\Downloads\City.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 2939729 ymin: 13748390 xmax: 3269290 ymax: 13976040
## Projected CRS: NAD83 / Texas South Central (ftUS)
polygon_data <- st_simplify(polygon_data, dTolerance = 0.01)
tm_shape(harris_tracts) +
tm_borders() +
tm_shape(polygon_data) +
tm_polygons(
col = "navy",
alpha = 0.5,
border.col = "black",
title = "City Polygon Data") +
tm_compass(
type = "8star",
position = c("right", "top"),
size = 1) +
tm_scale_bar(
position = c("left", "bottom")) +
tm_layout(
legend.position = c("center", "bottom"),
legend.outside = TRUE,
legend.outside.position = "right")
# 5. Have overlay analysis for any two shape file data layers in your study area
tm_shape(harris_tracts) +
tm_borders() +
tm_shape(polygon_data) +
tm_polygons(
col = "orange",
alpha = 0.5,
border.col = "black",
title = "Polygons") +
tm_dots(
col = "blue",
size = 0.2,
title = "Intersected Points") +
tm_compass(
type = "8star",
position = c("right", "top"),
size = 1) +
tm_scale_bar(
position = c("left", "bottom")) +
tm_layout(
legend.position = c("center", "bottom"),
legend.outside = TRUE,
legend.outside.position = "right")
# 6. Get Digital Elevation Model (DEM) data of a selected area of the county
dem <- raster("C:/Users/cruzs/Downloads/USGS_1_n30w096_20240229.tif")
gplot(dem) +
geom_tile(aes(fill = value)) +
scale_fill_viridis_c(option = "C", name = "Elevation (m)") +
labs(
title = "Digital Elevation Model",
x = "Longitude",
y = "Latitude"
) +
theme_minimal()
# 7. Read the DEM data and get the min and max of elevations
DEM_mm <- raster('C:/Users/cruzs/Downloads/USGS_1_n30w096_20240229.tif')
DEM_mm
## class : RasterLayer
## dimensions : 3612, 3612, 13046544 (nrow, ncol, ncell)
## resolution : 0.0002777778, 0.0002777778 (x, y)
## extent : -96.00167, -94.99833, 28.99833, 30.00167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=NAD83 +no_defs
## source : USGS_1_n30w096_20240229.tif
## names : Layer_1
DEM_mm <- setMinMax(DEM_mm)
DEM_mm
## class : RasterLayer
## dimensions : 3612, 3612, 13046544 (nrow, ncol, ncell)
## resolution : 0.0002777778, 0.0002777778 (x, y)
## extent : -96.00167, -94.99833, 28.99833, 30.00167 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=NAD83 +no_defs
## source : USGS_1_n30w096_20240229.tif
## names : Layer_1
## values : -11.40814, 80.06804 (min, max)
cellStats(DEM_mm, min)
## [1] -11.40814
cellStats(DEM_mm, max)
## [1] 80.06804
cellStats(DEM_mm, range)
## [1] -11.40814 80.06804
DEM_mm@crs
## Coordinate Reference System:
## Deprecated Proj.4 representation: NA
DEM_mm@extent
## class : Extent
## xmin : -96.00167
## xmax : -94.99833
## ymin : 28.99833
## ymax : 30.00167
hist(DEM_mm, main = "Distribution of Elevation Values",
col = "blue",
maxpixels = 22000000)
# 8. Make a map that shows areas with elevation lower than 350m
DEM_mm_below_350 <- calc(DEM_mm, fun = function(x) { ifelse(x < 350, x, NA) })
gplot(DEM_mm_below_350) +
geom_tile(aes(fill = value)) +
scale_fill_viridis_c() +
labs(title = "Areas with Elevation Below 350m", fill = "Elevation") +
theme_minimal()
# 9. Make a map that shows the slop degree of topography map is less than 3
slope_DEM <- terrain(DEM_mm, opt = "slope", unit = "degrees")
slope_below_3 <- calc(slope_DEM, fun = function(x) { ifelse(x < 3, x, NA) })
gplot(slope_below_3) +
geom_tile(aes(fill = value)) +
scale_fill_viridis_c() +
labs(title = "Areas with Slope Below 3 Degrees", fill = "Slope (°)") +
theme_minimal()
# 10. Make a map that shows the aspect criteria of topography map is less than between 22.5 and 157.5
aspect_DEM <- terrain(DEM_mm, opt = "aspect", unit = "degrees")
aspect_criteria <- calc(aspect_DEM, fun = function(x) {
ifelse(x >= 22.5 & x <= 157.5, x, NA)
})
aspect_criteria@crs
## Coordinate Reference System:
## Deprecated Proj.4 representation: NA
aspect_criteria@extent
## class : Extent
## xmin : -96.00167
## xmax : -94.99833
## ymin : 28.99833
## ymax : 30.00167
hist(aspect_criteria, main = "Distribution of Aspect Values (22.5° to 157.5°)",
col = "lightblue",
maxpixels = 22000000)
gplot(aspect_criteria) +
geom_tile(aes(fill = value)) +
scale_fill_viridis_c(option = "C", name = "Aspect (°)") +
labs(title = "Areas with Aspect Between 22.5° and 157.5°",) +
theme_minimal()
# 11. Have write-up to summarize your findings
# The findings from this analysis give insight into the demographic distribution, infrastructure and terrain characteristics of Harris County. The Hispanic population concentrations and access to facilities give opportunity for urban planning intervention. Also, the flat terrain and low elevations demonstrate the need for flood prevention improvements. Integrating these changes can allow for improve quality of life across the county.