GIS 4091 Project 4: R

#importing necessary libraries
library(ggplot2)
library(raster)
library(arcgisbinding)
library(sf)
#Validating connection to ArcGIS Pro.
arc.check_product()

product: ArcGIS Pro (13.1.0.41833)
license: Advanced
version: 1.0.1.300 

#Creating variable with path to the online shapefile
lava.dir <- "https://services.arcgis.com/nzS0F0zdNLvs7nc8/arcgis/rest/services/Mauna_Loa_Lava_Flows_200yrs_Wolfe_Morris_1996/FeatureServer/26"
#Retrieving the meta data from the shapefile
lava.meta <- arc.open(lava.dir)

#Creating variable with path to local raster file.
ncnp.dir <- "C:/Users/rmmay/GIS4091/FinalProject/nocavegmap/noca_3m.tif"
#Retrieving the meta data from the raster file.
ncnp.meta <- arc.open(humm.dir)
#Displaying raster meta data
ncnp.meta

dataset_type    : RasterDataset
path            : C:/Users/rmmay/GIS4091/FinalProject/nocavegmap/noca_3m.tif 
format          : TIFF
pixel_type      : U8 (8bit)
compression_type: LZW
nrow            : 30184
ncol            : 30068
extent          : xmin=595920, ymin=5339270, xmax=686124, ymax=5429822
WKT             : PROJCS["NAD_1983_UTM_Zone_10N",GEOGCS["GCS_North_American_19...
WKID            : 26910
bands           : 1

#Displaying raster meta data
lava.meta

dataset_type    : FeatureClass
path            : https://services.arcgis.com/nzS0F0zdNLvs7nc8/arcgis/rest/services/Mauna_Loa_Lava_Flows_200yrs_Wolfe_Morris_1996/FeatureServer/26 
fields          : OBJECTID, ID, ISLAND, VOLCANO, STRAT_CODE, SYMBOL, 
fields          : AGE_GROUP, AGE_RANGE, NAME, NAME0, UNIT, ROCK_TYPE, 
fields          : LITHOLOGY, VOLC_STAGE, COMPOSITIO, SOURCE, Shape__Area, 
fields          : Shape__Length, Shape
extent          : xmin=-17357396, ymin=2148749, xmax=-17264126, ymax=2260571
geometry type   : Polygon
WKT             : PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_...
WKID            : 3857 

#Creating an ArcGIS raster with raster meta data.
arcRaster <- arc.raster(ncnp.meta, nrow=1500, ncol=1000)
#Converting that ArcGIS raster into a spatial data frame
ncnp.sp <- arc.data2sp(arcRaster)
#Displaying spatial data frame
head(ncnp.sp)

#Creating a regualr data frame in order to utilize ggplot. 
#Removing all NA/Null values from data frame.
#Pulling out x and y coordinates from raster data.
ncnp.df <- as.data.frame(ncnp.sp, xy = TRUE, na.rm=TRUE)
#Displaying data frame.
head(ncnp.df)

#Plotting the raster data frame using ggplot. X and y coordinate to point values xmin and ymin. 
ggplot(data = ncnp.df, aes(x=xmin, y=ymin)) + 
  #Using geom_tile to accuratley disply the tiles of the raster from their centers
  #Defining the fill as raster values to accurately visualize the data
  geom_tile(aes(fill = ncnp.df$Band_1)) +
  #Setting title
  ggtitle("Vegetation Classification of North Cascades National Park")+
  #Setting colorramp and title for the color ramp
  scale_fill_viridis_c(name = "Vegetation Inventory", option = "magma")

  
#Selecting the shapefile meta data and loading it into a data frame
l.select <- arc.select(lava.meta)
#Turning the data frame into an sf object in order to utilize ggplot
lava.sf <- arc.data2sf(l.select)

#Plotting the raster utilizing ggplot
ggplot(data = lava.sf)+
  #Utilizing geom_sf to correctly display our sf object shapefile
  #Setting the fill to the AGE_RANGE column, so the polygons are all colored
  geom_sf(aes(fill= AGE_RANGE))+
  #Setting title
  ggtitle("Mauna Loa Lava Flows 1843-1984")

#Displaying the heads of our sf object and data frame
head(ncnp.df)

head(lava.sf)

Simple feature collection with 6 features and 18 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -17318710 ymin: 2202651 xmax: -17264130 ymax: 2240873
Projected CRS: WGS 84 / Pseudo-Mercator
  OBJECTID   ID ISLAND VOLCANO STRAT_CODE SYMBOL AGE_GROUP      AGE_RANGE
1        1 5522 Hawaii    mloa        209    Qk5         1      A.D. 1843
2        2 5715 Hawaii    mloa        209    Qk5         1 A.D. 1880-1881
3        3 5733 Hawaii    mloa        209    Qk5         1      A.D. 1855
4        4 5787 Hawaii    mloa        209    Qk5         1      A.D. 1935
5        5 5986 Hawaii    mloa        209    Qk5         1 A.D. 1880-1881
6        6 6058 Hawaii    mloa        209    Qk5         1      A.D. 1852
        NAME NAME0 UNIT  ROCK_TYPE       LITHOLOGY VOLC_STAGE
1 Kau Basalt            Lava flows Pahoehoe and aa     shield
2 Kau Basalt            Lava flows Pahoehoe and aa     shield
3 Kau Basalt            Lava flows Pahoehoe and aa     shield
4 Kau Basalt            Lava flows Pahoehoe and aa     shield
5 Kau Basalt            Lava flows Pahoehoe and aa     shield
6 Kau Basalt            Lava flows Pahoehoe and aa     shield
         COMPOSITIO                  SOURCE Shape__Area Shape__Length
1 Tholeiitic basalt Wolfe and Morris, 1996a    40143421    138870.740
2 Tholeiitic basalt Wolfe and Morris, 1996a    49217245    239640.430
3 Tholeiitic basalt Wolfe and Morris, 1996a    53277692    203196.521
4 Tholeiitic basalt Wolfe and Morris, 1996a    21233708    108935.882
5 Tholeiitic basalt Wolfe and Morris, 1996a    41195348    199844.653
6 Tholeiitic basalt Wolfe and Morris, 1996a     1226408      8535.583
                            geom
1 MULTIPOLYGON (((-17314368 2...
2 MULTIPOLYGON (((-17265278 2...
3 MULTIPOLYGON (((-17275771 2...
4 MULTIPOLYGON (((-17308705 2...
5 MULTIPOLYGON (((-17295692 2...
6 POLYGON ((-17281177 2232958...

#Plotting the raster data frame cell value count utilizing ggplot and geom_histogram
#Setting x equal to the band values.
ggplot(data = ncnp.df, aes(x = ncnp.df$Band_1))+
  #Utilizing geom_histogram to make a histogram and defining 10 bins.
  geom_histogram(bins=10)+
  #Setting Title
  ggtitle("Pixel Counts for North Cascades National Park Vegetation Inventory")+
  xlab("Pixel Values")

#Plotting sp object shapefile data; lava flow count per year.
ggplot(data=lava.sf)+
  #Utilizing geom_bar to create a bar graph showing lava flow count per age class.
  #Removing legend because it warps barplot.
  geom_bar(aes(x=AGE_RANGE, fill=AGE_RANGE), show.legend = FALSE)+
  #Adding title
  ggtitle("Count of Mauna Loa Lava Flows between 1843-1984")+
  #Preventing x axis labels from overlapping.
  scale_x_discrete(guide = guide_axis(n.dodge=3))

#Importing necessary library for st_extract.
library(stars)

#Reading in the North cascades tif in with stars.
st.ncnp <- read_stars(humm.dir)
#Displaying the data.
st.ncnp

stars_proxy object with 1 attribute in 1 file(s):
$noca_3m.tif
[1] "[...]/noca_3m.tif"

dimension(s):

#Creating sample points from our raster data within the same dimensions of the raster.
ncnp.pnt <- st_sample(st_as_sfc(st_bbox(st.humm)), 10)
#Displaying the points and their values.
humm.sfc

Geometry set for 10 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 606403.7 ymin: 5339547 xmax: 675429 ymax: 5428699
Projected CRS: NAD83 / UTM zone 10N
First 5 geometries:

POINT (647165.7 5368074)
POINT (668002.6 5373423)
POINT (619313.8 5404168)
POINT (606403.7 5339547)
POINT (669409.7 5344235)

#Utilzing st_extract to extract cell values from the North Cascades raster from the random points within the rasters dimensions.
st_extract(st.ncnp, ncnp.pnt)

Simple feature collection with 10 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 602773.3 ymin: 5340046 xmax: 675885.2 ymax: 5422087
Projected CRS: NAD83 / UTM zone 10N
   noca_3m.tif                 geometry
1            7 POINT (663837.5 5367559)
2           NA POINT (630757.5 5340046)
3           NA POINT (671979.9 5386416)
4           NA POINT (661000.9 5410456)
5           NA POINT (602773.3 5379630)
6           44 POINT (630597.1 5377925)
7           27 POINT (636821.9 5422087)
8           NA POINT (675885.2 5420633)
9           NA   POINT (675590 5397761)
10          23   POINT (634078 5408920)

#Writing the new points and their values to a ESRI Geodatabase.
arc.write("C:/Users/rmmay/GIS4091/FinalProject/nocavegmap/points_val.gdb/fc", ncnp.pnt)