- I decided to project my data from WGS 84 to Alaska Albers because all of the data is in Alaska. I only show they crs description code for the first one because they are all transformed to the same projection.
## Reading layer `Qfaults_US_Database' from data source `C:\USU\Spring_2021\Earthquake_Finals\Earthquake_Finals\Qfaults_GIS\SHP\Qfaults_US_Database.shp' using driver `ESRI Shapefile'
## Simple feature collection with 112809 features and 22 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: -180 ymin: 18.76623 xmax: 180 ymax: 70.7517
## geographic CRS: WGS 84
## Coordinate Reference System:
## User input: EPSG:3338
## wkt:
## PROJCRS["NAD83 / Alaska Albers",
## 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["Alaska Albers (meters)",
## METHOD["Albers Equal Area",
## ID["EPSG",9822]],
## PARAMETER["Latitude of false origin",50,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",-154,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",55,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",65,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["USA - Alaska"],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]
## Reading layer `Volcano_Clip' from data source `C:\USU\Spring_2021\Earthquake_Finals\Earthquake_Finals\Volcano_Clip.shp' using driver `ESRI Shapefile'
## Simple feature collection with 89 features and 10 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -178.7926 ymin: 51.78852 xmax: 179.5773 ymax: 65.59928
## geographic CRS: WGS 84
## Reading layer `Alaska_Earthquakes_2020_Layer' from data source `C:\USU\Spring_2021\Earthquake_Finals\Alaska_Earthquakes_2020_Layer.shp' using driver `ESRI Shapefile'
## Simple feature collection with 6439 features and 22 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -179.998 ymin: 49.6097 xmax: 179.9806 ymax: 70.1144
## geographic CRS: WGS 84
## Reading layer `Volcano_Buffer' from data source `C:\USU\Spring_2021\Earthquake_Finals\Earthquake_Finals\Volcano_Buffer.shp' using driver `ESRI Shapefile'
## Simple feature collection with 89 features and 12 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -180 ymin: 50.4403 xmax: 180 ymax: 66.94444
## geographic CRS: WGS 84
## Reading layer `GU_StateOrTerritory' from data source `C:\USU\Spring_2021\Earthquake_Finals\Earthquake_Finals\GOVTUNIT_Alaska_State_Shape\Shape\GU_StateOrTerritory.shp' using driver `ESRI Shapefile'
## Simple feature collection with 2 features and 17 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -179.2311 ymin: 51.17509 xmax: 179.8597 ymax: 71.43979
## geographic CRS: WGS 84
- I took a layer that has earthquake points in Alaska and another layer that has volcano points with a 150km buffer around them and wanted to connect them to see how many earthquakes were within that buffer distance. I wanted to see if there were a lot of earthquakes connected to volcano location, because that could be a cause potentially. I commented the last line because they output was very long and ugly, but it shows which points intersect the buffers.
# Using intersects with and without sparse to find out where earthquake points fall in the buffer area, which is made up of polygons.
Earthquake_Volcano_Buffer_Join = st_intersects(Earthquakes_Project, Volcano_Buffer_Project)
mat = st_intersects(Earthquakes_Project, Volcano_Buffer_Project, sparse = F)
#apply(mat, 1, any)
- I first chose to look at earthquakes with only a magnitude greater than 4 because the majority of them were less than that, and I wanted to see if there was a trend in geographical location for the larger scale earthquakes. The other subset is looking at depth of the epicenters and I chose to look at only those deeper than 150km because they were also the end of more extreme in my dataset. Most of the outputs were the same earthquakes or in similar locations, making me consider a trend that I want to look into more.
Magnitude = Earthquakes_Project %>%
filter(mag > 4)
Depth = Earthquakes_Project %>%
filter(depth > 150)
ggplot()+
geom_sf(data = Magnitude)+
geom_sf(data = State_Boundary_Project)
ggplot()+
geom_sf(data = Depth)+
geom_sf(data = State_Boundary_Project)
- In this new spatial dataframe subset I wanted to look at only volcanoes that were labeled as stratovolcano.
Strato = Volcanoes_Project %>%
filter(TYPE == "Stratovolcano")
class(Strato)
## [1] "sf" "data.frame"
#mapview(Strato) Wouldn't knit with mapview, but is a better visualization of the new dataframe.
- This map is showing volcanoes and earthquakes as of the year 2020. The volcanoes have a buffer of 150km shown. With this visualization I was hoping to see any connections between the earthquakes of magnitude higher than 4 and volcano location. While some volcanoes do not have any earthquakes surrounding them, the majority of earthquakes are near a volcano. There is a cluster along the island chain of volcanoes and earthquakes, showing a relationship in geography.
ggplot() +
geom_sf(data = State_Boundary_Project, color = "black", aes(fill = "State Boundary", State_Boundary_Project = "title"), show.legend = "polygon") +
geom_sf(data = Volcano_Buffer_Project, color = "purple", aes(fill = "Volcano 150km Buffer", Volcano_Buffer_Project = "title"), show.legend = "polygon") +
geom_sf(data = Volcanoes_Project, color = "red", aes(fill = "Volcanoes", Volcanoes_Project = "title"), show.legend = "polygon") +
geom_sf(data = Magnitude, color = "blue", aes(fill = "Earthquakes with Magnitude 4+", Magnitude = "title"), show.legend = "polygon") +
annotation_north_arrow(location = "br", which_north = "true",
pad_x = unit(0.05, "in"), pad_y = unit(0.05, "in"),
height = unit(0.4, "in"), width = unit(0.25, "in"),
style = north_arrow_fancy_orienteering) +
labs(title = "Alaska Earthquakes and Volcanic Activity 2020") +
theme(panel.grid.major = element_line(color = grey(.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "light grey"), axis.text.x=element_blank(), axis.ticks.x=element_blank(), axis.text.y=element_blank(), axis.ticks.y=element_blank()) +
scale_fill_manual(values = c("State Boundary" = "black", "Earthquakes with Magnitude 4+" = "blue", "Volcano 150km Buffer" = "purple", "Volcanoes" = "red"),
guide = guide_legend(override.aes =
list(linetype = "blank", shape = NA)),
name = "Legend")
- This visualization shows fault lines and the same earthquakes (magnitude higher than 4). There are a few smaller fault lines that do not have many large earthquakes surrounding them, but there is a clear large fault line that is mirrored by earthquakes along the island chain. This shows a relationship between faults and earthquakes.
tm_shape(Magnitude) +
tm_dots(col = "blue",
legend.show = TRUE) +
tm_shape(Faults_Project) +
tm_lines(col = "green",
legend.col.show = TRUE) +
tm_shape(State_Boundary_Project) +
tm_polygons(col = "black",
legend.show = TRUE)
