ant109_EC
Alexander Marquez
2024-06-13
These are the sites from where I obtained the shapefiles: Alaska Boundary and Reserves: https://www.sciencebase.gov/catalog/item/59d5b565e4b05fe04cc53a91 Lakes in Alaska: https://learning.nceas.ucsb.edu/2020-11-RRCourse/session-13-geospatial-analysis-in-r.html
Below is the sites which I used to help develop my code: Cite: https://learning.nceas.ucsb.edu/2020-11-RRCourse/session-13-geospatial-analysis-in-r.html Cite2: https://r-charts.com/spatial/interactive-maps-leaflet/
## In order to read shapefile in R
#install.packages("tidyverse")
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors#install.packages("sf")
library(sf) ## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE#install.packages("maptools")
#library(maptools) # creates maps and work with spatial files
#install.packages("broom")
library(broom) # assists with tidy data
#install.packages("ggplot2")
library(ggplot2) # graphics package
#install.packages("leaflet")
library(leaflet) # interactive graphics (output does not show in RMD files)
#install.packages("dplyr")
library(dplyr) # joining data frames
#install.packages("readr")
library(readr) # quickly reads files into R
library(scales)##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factor#install.packages("ggmap")
library(ggmap)## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
## Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
## OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.## To make the interactive map
#install.packages("leaflet")
library(leaflet)AlaskaBoundary <- st_read("/Users/misterx/GOVTUNIT_Alaska_State_Shape/Shape/GU_StateOrTerritory.shp",
stringsAsFactors = FALSE)## Reading layer `GU_StateOrTerritory' from data source
## `/Users/misterx/GOVTUNIT_Alaska_State_Shape/Shape/GU_StateOrTerritory.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.2311 ymin: 51.17509 xmax: 179.8597 ymax: 71.43979
## Geodetic CRS: NAD83ProtectedAreas <- st_read("/Users/misterx/Alaska_DEC_Sensitive_Areas/Alaska_DEC_Sensitive_Areas.shp",
stringsAsFactors = FALSE)## Reading layer `Alaska_DEC_Sensitive_Areas' from data source
## `/Users/misterx/Alaska_DEC_Sensitive_Areas/Alaska_DEC_Sensitive_Areas.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 19335 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -179.1501 ymin: 51.20992 xmax: 179.7717 ymax: 70.32941
## Geodetic CRS: WGS 84Rivers <- st_read("/Users/misterx/ant109_MapOfAlaska_ShpFolder/NHD_Major_Rivers/NHD_Major_Rivers.shp",
stringsAsFactors = FALSE)## Reading layer `NHD_Major_Rivers' from data source
## `/Users/misterx/ant109_MapOfAlaska_ShpFolder/NHD_Major_Rivers/NHD_Major_Rivers.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 57007 features and 16 fields
## Geometry type: LINESTRING
## Dimension: XYZM
## Bounding box: xmin: -124.409 ymin: 32.50998 xmax: -114.1308 ymax: 42.54999
## z_range: zmin: -2e-04 zmax: 0
## m_range: mmin: 0 mmax: 100
## Geodetic CRS: NAD83 + NAVD88 heightRivers_new <- st_read("/Users/misterx/Downloads/shapefile_demo_data/ak_rivers_simp.shp",
stringsAsFactors = FALSE)## Reading layer `ak_rivers_simp' from data source
## `/Users/misterx/Downloads/shapefile_demo_data/ak_rivers_simp.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 39 features and 3 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -622169.8 ymin: 557375.7 xmax: 1554975 ymax: 2321138
## Projected CRS: Albers#plot(AlaskaBoundary)
st_crs(AlaskaBoundary)## 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]]#Alaska
AlaskaShape <- AlaskaBoundary %>%
st_transform(crs = 3338)
st_crs(AlaskaShape)## 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["Topographic mapping (small scale)."],
## AREA["United States (USA) - Alaska."],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]#Protected Areas
ProtectedAreas_3338 <- ProtectedAreas %>%
st_transform(crs = 3338)
st_crs(ProtectedAreas_3338)## 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["Topographic mapping (small scale)."],
## AREA["United States (USA) - Alaska."],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]#Rivers
Rivers_3338 <- Rivers %>%
st_transform(crs = 3338)
st_crs(Rivers_3338)## 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["Topographic mapping (small scale)."],
## AREA["United States (USA) - Alaska."],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]AlaskaProtected<- st_join(ProtectedAreas_3338, AlaskaShape, join=st_within)
head(AlaskaProtected)## Simple feature collection with 6 features and 30 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -601111.3 ymin: 574720.5 xmax: 568464.7 ymax: 2275120
## Projected CRS: NAD83 / Alaska Albers
## OBJECTID IFWSNo NWRName NWRCode UnitName UnitCode SubunitNam
## 1 1 992 Alaska Peninsula akp Ugashik uga <NA>
## 2 2 992 Alaska Peninsula akp Chignik chi <NA>
## 3 3 992 Alaska Peninsula akp Pavlof pav <NA>
## 4 4 992 Alaska Peninsula akp North Creek ncr <NA>
## 5 5 993 Arctic arc <NA> <NA> <NA>
## 6 6 0 not applicable na <NA> <NA> <NA>
## IslandName Source Updated
## 1 <NA> USGS paper quads, BLM PLSS 2009-02-20
## 2 <NA> USGS paper quads, BLM PLSS 2009-02-20
## 3 <NA> USGS paper quads, BLM PLSS 2009-02-20
## 4 <NA> USGS paper quads, BLM PLSS 2009-02-20
## 5 <NA> USGS paper quads, BLM PLSS 2009-02-25
## 6 <NA> COGO - USS 4234 2009-05-13
## Notes
## 1 Updated 02/20/2009 by Scott McGee to conform to the joint FWS/NPS common boundary between Alaska Peninsula NWR and Aniakchak NP. ANY ALTERATIONS OF THE FWS/NPS COMMON BOUNDARY LINE MUST BE DONE IN COLLABORATION WITH THE NATIONAL PARK SERVICE.
## 2 Updated 02/20/2009 by Scott McGee to conform to the joint FWS/NPS common boundary between Alaska Peninsula NWR and Aniakchak NP. ANY ALTERATIONS OF THE FWS/NPS COMMON BOUNDARY LINE MUST BE DONE IN COLLABORATION WITH THE NATIONAL PARK SERVICE.
## 3 Last update: snapped arcs that follow section lines to the new NPS PLSS
## 4 Last update: snapped arcs that follow section lines to the new NPS PLSS
## 5 Last update: snapped arcs that follow section lines to the new NPS PLSS; snapped Alaska/Yukon boundary arc to the IBC official NAD83 coordinates (see L:\\GeoData\\Boundaries\\141st Meridian\\International Boundary Commission Official NAD83 141st Meridian Coo
## 6 Last update: USS 4234 entered via COGO. Village exclusion (Kaktovik)
## SHAPE_Leng ShapeSTAre ShapeSTLen permanent_
## 1 571373.435 5661559996 571373.435 0c2b464d-3e76-4be9-bc7c-0af54214cc86
## 2 816001.370 8356051605 816001.370 0c2b464d-3e76-4be9-bc7c-0af54214cc86
## 3 817502.741 5852082084 817502.741 0c2b464d-3e76-4be9-bc7c-0af54214cc86
## 4 30734.582 34287529 30734.582 0c2b464d-3e76-4be9-bc7c-0af54214cc86
## 5 1866914.756 80346360344 1866914.756 <NA>
## 6 4455.602 548801 4455.602 0c2b464d-3e76-4be9-bc7c-0af54214cc86
## source_fea source_dat
## 1 02 74a33eda-5543-40af-a9c6-293ab88466f4
## 2 02 74a33eda-5543-40af-a9c6-293ab88466f4
## 3 02 74a33eda-5543-40af-a9c6-293ab88466f4
## 4 02 74a33eda-5543-40af-a9c6-293ab88466f4
## 5 <NA> <NA>
## 6 02 74a33eda-5543-40af-a9c6-293ab88466f4
## source_d_1 source_ori
## 1 2022 TIGER/Line Shapefile, Current State and Equivalent U.S. Census Bureau
## 2 2022 TIGER/Line Shapefile, Current State and Equivalent U.S. Census Bureau
## 3 2022 TIGER/Line Shapefile, Current State and Equivalent U.S. Census Bureau
## 4 2022 TIGER/Line Shapefile, Current State and Equivalent U.S. Census Bureau
## 5 <NA> <NA>
## 6 2022 TIGER/Line Shapefile, Current State and Equivalent U.S. Census Bureau
## loaddate fcode state_fips state_name population areasqkm gnis_id
## 1 2022-11-17 61100 02 Alaska 733391 1724322 1785533
## 2 2022-11-17 61100 02 Alaska 733391 1724322 1785533
## 3 2022-11-17 61100 02 Alaska 733391 1724322 1785533
## 4 2022-11-17 61100 02 Alaska 733391 1724322 1785533
## 5 <NA> NA <NA> <NA> NA NA <NA>
## 6 2022-11-17 61100 02 Alaska 733391 1724322 1785533
## gnis_name shape_Leng shape_Area ObjectID geometry
## 1 State of Alaska NA NA 1 MULTIPOLYGON (((-156900.6 8...
## 2 State of Alaska NA NA 1 MULTIPOLYGON (((-259186.7 7...
## 3 State of Alaska NA NA 1 MULTIPOLYGON (((-423667.9 6...
## 4 State of Alaska NA NA 1 MULTIPOLYGON (((-526429.9 6...
## 5 <NA> NA NA NA MULTIPOLYGON (((412747.9 22...
## 6 State of Alaska NA NA 1 MULTIPOLYGON (((399522.7 22...AlaskaClean <- AlaskaProtected %>%
group_by(NWRName) %>%
summarize(Area = mean(areasqkm), do_union = TRUE)
group_mean<- aggregate(x= AlaskaProtected$areasqkm,
# Specify group indicator
by = list(AlaskaProtected$NWRName),
# Specify function (i.e. mean)
FUN = mean, na.rm=TRUE)
print(group_mean)## Group.1 x
## 1 Alaska Maritime 1724322
## 2 Alaska Peninsula 1724322
## 3 Arctic NaN
## 4 Becharof 1724322
## 5 Innoko 1724322
## 6 Izembek 1724322
## 7 Kanuti 1724322
## 8 Kenai 1724322
## 9 Kodiak 1724322
## 10 Koyukuk 1724322
## 11 not applicable 1724322
## 12 Nowitna 1724322
## 13 Selawik 1724322
## 14 Tetlin 1724322
## 15 Togiak 1724322
## 16 Yukon Delta 1724322
## 17 Yukon Flats 1724322ProtectedSurfaceArea <- AlaskaProtected %>%
as.data.frame() %>%
group_by(NWRName) %>%
summarise(surface_area = mean(areasqkm))
ProtectedSurfaceArea## # A tibble: 17 × 2
## NWRName surface_area
## <chr> <dbl>
## 1 Alaska Maritime NA
## 2 Alaska Peninsula 1724322.
## 3 Arctic NA
## 4 Becharof 1724322.
## 5 Innoko 1724322.
## 6 Izembek 1724322.
## 7 Kanuti 1724322.
## 8 Kenai 1724322.
## 9 Kodiak 1724322.
## 10 Koyukuk 1724322.
## 11 Nowitna 1724322.
## 12 Selawik 1724322.
## 13 Tetlin 1724322.
## 14 Togiak NA
## 15 Yukon Delta NA
## 16 Yukon Flats 1724322.
## 17 not applicable 1724322.AlaskaJoin <- left_join(AlaskaProtected,ProtectedSurfaceArea)## Joining with `by = join_by(NWRName)`names(AlaskaJoin)## [1] "OBJECTID" "IFWSNo" "NWRName" "NWRCode" "UnitName"
## [6] "UnitCode" "SubunitNam" "IslandName" "Source" "Updated"
## [11] "Notes" "SHAPE_Leng" "ShapeSTAre" "ShapeSTLen" "permanent_"
## [16] "source_fea" "source_dat" "source_d_1" "source_ori" "loaddate"
## [21] "fcode" "state_fips" "state_name" "population" "areasqkm"
## [26] "gnis_id" "gnis_name" "shape_Leng" "shape_Area" "ObjectID"
## [31] "surface_area" "geometry"#plot(AlaskaJoin["surface_area"]) pal <- leaflet::colorNumeric(palette = "Reds", domain = AlaskaJoin$surface_area)
AlaskaPro <- ggplot(AlaskaJoin) +
geom_sf(data = Rivers_new, aes(size = StrOrder), color = "blue") +
geom_sf(data = AlaskaProtected, aes(fill = "red")) #we want to add protected areas name and geo area
AlaskaPro
epsg3338 <- leaflet::leafletCRS(
crsClass = "L.Proj.CRS",
code = "EPSG:3338",
proj4def = "+proj=aea +lat_1=55 +lat_2=64 +lat_0=50 +lon_0=-152 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
resolutions = 2^(16:7))
st_crs(AlaskaJoin)## 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["Topographic mapping (small scale)."],
## AREA["United States (USA) - Alaska."],
## BBOX[51.3,172.42,71.4,-129.99]],
## ID["EPSG",3338]]AlaskaJoinNew <- AlaskaJoin %>% st_transform(crs = 4326)
#Alaska_basemap <- leaflet() %>%
#
# setView(lng = -152, lat = 64, zoom =4)
#Alaska_basemap
AlaskaPro +
leaflet(options = leafletOptions(crs = epsg3338)) %>%
addTiles() %>%
setView(lng = -152, lat = 64, zoom =4) %>%
addPolygons(data = AlaskaJoinNew,
fillColor = ~pal(surface_area),
weight = 1,
color = "red",
fillOpacity = 1,
label = ~paste0(AlaskaJoin$NRWName, ", surface area ", AlaskaJoin$surface_area)) %>%
scale_size(range = c(0.01, 0.2), guide = F) +
theme_bw() +
labs(fill = "Alaska Protected Areas")## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.