# Load necessary libraries
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(sf)
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(tigris)
## Warning: package 'tigris' was built under R version 4.4.3
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.4.3
##
## 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(readxl)
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.4.3
library(htmlwidgets)
## Warning: package 'htmlwidgets' was built under R version 4.4.3
library(webshot)
## Warning: package 'webshot' was built under R version 4.4.3
# Load New Jersey counties shapefile
nj_counties <- tigris::tracts(state = "NJ", year = 2020, class = "sf", cb = TRUE)
## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 6% | |====== | 8% | |======= | 10% | |======== | 11% | |========= | 12% | |========== | 14% | |=========== | 16% | |============ | 17% | |============= | 18% | |============== | 20% | |=============== | 22% | |================ | 22% | |================== | 25% | |=================== | 27% | |==================== | 29% | |===================== | 30% | |======================= | 32% | |======================== | 34% | |========================= | 36% | |========================== | 37% | |=========================== | 39% | |============================= | 41% | |============================== | 43% | |=============================== | 44% | |================================ | 46% | |================================== | 48% | |=================================== | 50% | |==================================== | 51% | |===================================== | 53% | |======================================= | 55% | |======================================== | 57% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 69% | |================================================== | 71% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 76% | |======================================================= | 78% | |======================================================== | 80% | |========================================================= | 81% | |========================================================== | 83% | |=========================================================== | 85% | |============================================================= | 87% | |============================================================== | 88% | |=============================================================== | 90% | |================================================================ | 92% | |================================================================== | 94% | |=================================================================== | 95% | |==================================================================== | 97% | |===================================================================== | 99% | |======================================================================| 100%
# Load the CCC map data
CCC_map <- read_excel("C:/Users/crook/OneDrive - go.Stockton.edu/American Studies articles/1930s/assignments/CCC map.xlsx")
#View(CCC_map)
str(CCC_map)
## tibble [69 × 12] (S3: tbl_df/tbl/data.frame)
## $ PROJECT : chr [1:69] "BS-2" "SP-9" "S-58" "S-51" ...
## $ CO. # : num [1:69] 1269 1275 1278 218 1264 ...
## $ Designation : chr [1:69] "-C" "-C" NA NA ...
## $ DATE : POSIXct[1:69], format: "1939-10-31" "1933-10-30" ...
## $ RAILROAD : chr [1:69] "Absecon" "Berlin" "Blairstown" "Branchville" ...
## $ POST OFFICE : chr [1:69] "Absecon" "Berlin" "Blairstown" "Branchville" ...
## $ State : chr [1:69] "New Jersey" "New Jersey" "New Jersey" "New Jersey" ...
## $ LOCATION : chr [1:69] "Nacote Creek 7 mi N" "Camp 1 mi W" NA "Kittatinny Lake 5 mi NW" ...
## $ Lat : num [1:69] 39.4 NA NA NA NA ...
## $ Long : num [1:69] -74.5 NA NA NA NA ...
## $ Code Definition: chr [1:69] "Biological Study" "State Park" "State Forest" "State Forest" ...
## $ tidbit : chr [1:69] "Is this a bridge study?" NA NA NA ...
# Ensure Lat and Long are numeric and filter out NAs
CCC_map <- CCC_map %>%
mutate(Lat = as.numeric(Lat), Long = as.numeric(Long)) %>%
filter(!is.na(Lat) & !is.na(Long)) # Remove rows with NA in Lat or Long
# Sample list of post offices to plot
selected_post_offices <- c("Absecon", "New Gretna")
# Filter the CCC_map data
filtered_data <- CCC_map %>%
filter(`POST OFFICE` %in% selected_post_offices) %>%
mutate(Lat = as.numeric(Lat), Long = as.numeric(Long)) %>%
filter(!is.na(Lat) & !is.na(Long)) # Ensure no NAs
Create a spatial dataframe from your tibble
data_sf <- st_as_sf(CCC_map, coords = c("Long", "Lat"), crs = 4326, agr = "constant")
Check if data_sf is created correctly
print(st_geometry(data_sf)) # This should show your geometry
## Geometry set for 49 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -75.215 ymin: 39.0842 xmax: -73.9707 ymax: 41.1964
## Geodetic CRS: WGS 84
## First 5 geometries:
## POINT (-74.4855 39.4268)
## POINT (-74.346 40.9923)
## POINT (-74.346 40.9923)
## POINT (-74.346 40.9923)
## POINT (-74.5408 40.0155)
# Create a spatial dataframe from the filtered tibble
filtered_data_sf <- st_as_sf(filtered_data, coords = c("Long", "Lat"), crs = 4326, agr = "constant")
# Ensure Lat and Long are numeric and filter out NAs
CCC_map <- CCC_map %>%
mutate(Lat = as.numeric(Lat), Long = as.numeric(Long)) %>%
filter(!is.na(Lat) & !is.na(Long)) # Remove rows with NA in Lat or Long
# Create a spatial dataframe from the entire tibble
data_sf <- st_as_sf(CCC_map, coords = c("Long", "Lat"), crs = 4326, agr = "constant")
st_crs(nj_counties)
## 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(data_sf)
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
data_sf <- st_transform(data_sf, st_crs(nj_counties))
# Ensure your data is in WGS84
data_sf <- st_transform(data_sf, crs = 4326)
# Create an interactive map
map <- leaflet(data_sf) %>%
addTiles() %>% # Base map
addCircleMarkers(color = "blue", radius = 5,
popup = ~paste("Project:", PROJECT, "<br>",
"Post Office:", `POST OFFICE`, "<br>",
"Tidbit:", `tidbit`),
fillOpacity = 0.7) %>%
setView(lng = mean(st_coordinates(data_sf)[, 1]),
lat = mean(st_coordinates(data_sf)[, 2]),
zoom = 8) %>%
addMiniMap() %>% # Adds a mini map for easier navigation
addProviderTiles(providers$CartoDB.Positron) # Optional: Different basemap
# Print the map
map
# Save the map as an HTML file
#saveWidget(map, "map.html", selfcontained = TRUE)
# Use webshot to capture the HTML map as a PNG
#webshot("map.html", file = "map.png", vwidth = 800, vheight = 600)
# Group by the POST OFFICE column and summarize
grouped_data <- CCC_map %>%
group_by(`Designation`) %>% # Use backticks for column names with spaces
summarise(
count = n(), # Count of entries for each post office
average_lat = mean(Lat, na.rm = TRUE), # Average latitude
average_long = mean(Long, na.rm = TRUE) # Average longitude
)
# View the grouped data
print(grouped_data)
## # A tibble: 3 × 4
## Designation count average_lat average_long
## <chr> <int> <dbl> <dbl>
## 1 -C 10 39.8 -74.6
## 2 -V 5 40.0 -74.7
## 3 <NA> 34 40.4 -74.5
# Group by the POST OFFICE column and summarize
grouped_data <- CCC_map %>%
group_by(`Code Definition`) %>% # Use backticks for column names with spaces
summarise(
count = n()
) %>%
arrange(desc(count)) # Count of entries for each post office
#average_lat = mean(Lat, na.rm = TRUE), # Average latitude
#average_long = mean(Long, na.rm = TRUE) # Average longitude
# View the grouped data
print(grouped_data)
## # A tibble: 11 × 2
## `Code Definition` count
## <chr> <int>
## 1 State Park 13
## 2 State Forest 12
## 3 MC 7
## 4 Private Forest 4
## 5 Soil Conservation 4
## 6 Army 3
## 7 Air Force 2
## 8 ACH 1
## 9 Biological Study 1
## 10 County Park 1
## 11 Forest 1
