library(sf) # For spatial data manipulation
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(tmap) # For visualization
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(raster) # For working with raster data
## Loading required package: sp
library(dplyr) # For data manipulation
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidycensus) # Load tidycensus for accessing U.S. Census Bureau data directly
storm_channels <- st_read("C:/Users/dagob/OneDrive/Desktop/Project part II/StormChannels/StormChannels.shp")
## Reading layer `StormChannels' from data source
## `C:\Users\dagob\OneDrive\Desktop\Project part II\StormChannels\StormChannels.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 12047 features and 17 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 2024999 ymin: 13600410 xmax: 2249487 ymax: 13824430
## Projected CRS: NAD83 / Texas South Central (ftUS)
corridor_plans <- st_read("C:/Users/dagob/OneDrive/Desktop/Project part II/CorridorPlans/CorridorPlans.shp")
## Reading layer `CorridorPlans' from data source
## `C:\Users\dagob\OneDrive\Desktop\Project part II\CorridorPlans\CorridorPlans.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 12 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 2058688 ymin: 13665110 xmax: 2178301 ymax: 13793820
## Projected CRS: NAD83 / Texas South Central (ftUS)
### Qustion 1 Transform the coordinate reference systems
# Check the CRS of both datasets
### Qustion 1 Transform the coordinate reference systems
# Check the CRS of both datasets
st_crs(storm_channels) # Check the CRS of Storm Channels
## 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]]
st_crs(corridor_plans) # Check the CRS of Corridor Plans
## 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]]
# Transform CRS to EPSG: 4326 if needed
storm_channels <- st_transform(storm_channels, crs = 4326)
corridor_plans <- st_transform(corridor_plans, crs = 4326)
# Check the CRS of both datasets again
st_crs(storm_channels) # Check the CRS of Storm Channels
## 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)"],
## 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]]
st_crs(corridor_plans) # Check the CRS of Corridor Plans
## 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)"],
## 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]]
# Question 2 Find what parts of storm channels are within the corridor plan area
# Perform spatial intersection to find storm channels within corridor plans
intersected_channels <- st_intersection(storm_channels, corridor_plans)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# Question 3 Visualize the intersected storm channels
# Load Bexar County census tracts data
bexar_tracts <- get_acs(
geography = "tract",
variables = "B01003_001", # Total population variable
state = "TX",
county = "Bexar",
geometry = TRUE,
year = 2021
)
## 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% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |================ | 22% | |================ | 23% | |================= | 24% | |================== | 26% | |==================== | 28% | |===================== | 30% | |====================== | 32% | |============================================= | 64% | |====================================================== | 77% | |====================================================== | 78% | |======================================================== | 80% | |=============================================================== | 91% | |================================================================= | 94% | |=================================================================== | 95% | |======================================================================| 100%
# Visualize the intersected storm channels with Bexar County census tracts as the background
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(bexar_tracts) +
tm_borders() + # Show census tract borders
tm_shape(intersected_channels) +
tm_lines(col = "blue", lwd = 2) + # Show intersected storm channels in blue
tm_layout(title = "Intersected Storm Channels in Bexar County")
# Question 4 Read the raster file
# Load the raster file (land surface temperature)
LST <- raster("C:/Users/dagob/OneDrive/Desktop/Project part II/LC09_028040_20230826.tif")
# Convert the Kelvin into Celsius
LST_Celsius = LST - 273.15
plot(LST_Celsius, main="Land Surface Temperature of Dallas")
# Question 5 Only visualize the areas with temperatures higher than 40 Degree Celsius
# Mask areas with temperatures >40°C for use with plot()
LST_above_40 <- LST_Celsius
LST_above_40[LST_above_40 < 40] <- NA # Mask areas with temperatures <40°C
# Plot using the filtered raster
plot(LST_above_40,
main = "Land Surface Temperature > 40°C",
col = heat.colors(50))
