Mt. Baker ⛰️

Homework #1

Markdown Author: Jessie Bell

Download this Rmd: Top right corner → Code → Download Rmd

My Baker Map

All code reasoning is documented in the code below.

#load the things

mtbDEM <- rast("Mt_Baker/data/mtbDEM.tif") #you need data in raster form here to run through spatraster code below. 
mtb_hill <- rast("Mt_Baker/data/mtbHill.tif")
lines <- st_read("Mt_Baker/data/mtbChairLines.shp")
points <- st_read("Mt_Baker/data/mtbLodges.csv")

#check to see what is in the mtb_hill data
baker_hill_df <- as.data.frame(mtb_hill, xy=T)

#do it again with DEM
baker_DEM_df <- as.data.frame(mtbDEM, xy=T)

#create dataframe for points
points <- as.data.frame(points, xy=T)

#create sf for points (use this for ggplot)
lodges <- st_as_sf(points, coords = c("X", "Y"), crs = st_crs(32610))

#make pretty colors and ramp for continuous
cols <- c("#00007f", "#ff00ff", "#01ffff")
elevcols <- colorRampPalette(cols)(100)


#changing these for legend titles
lodges$Lodge <- as.factor(lodges$id)
lines$Lifts <- as.character(c("Chair 1", "Chair 2", "Chair 3", "Chair 4", "Chair 5", "Chair 6", "Chair 7", "Chair 8"))

#create ggplot

#1 call in hillshade and use gray colors with alpha for transparency

p1 <- ggplot()+
  geom_spatraster(data=mtb_hill, mapping = aes(fill=value))+
  scale_fill_gradientn(colors = gray.colors(100, start = 0.0, end=0.988), guide="none") 

#2 Call in DEM and use color ramp that you created above 
p2 <- p1 + new_scale_fill()+
  geom_spatraster(data=mtbDEM, mapping=aes(fill=elev), show.legend = F)+
  scale_fill_viridis_b(values=cols, alpha=0.3)

#3 Add in line data 
p3 <- p2 + geom_sf(data=lines, mapping=aes(color=Lifts), linewidth=1)

#4 Add in point data
p4 <- p3 + geom_sf(data=lodges,show.legend = F)+
  geom_sf_label(mapping = aes(label = Lodge), data=lodges)

#add in contours 
p5 <- p4 + geom_contour(data=baker_DEM_df,aes(x=x,y=y, z=as.integer(elev)),color="white", alpha=0.6, linewidth=0.3) + #addded contours
  guides(col= guide_legend(title= "", override.aes = list(fill=NA)))+
  theme(legend.position = "bottom")+
  labs(x="", y="", caption = "Elevation contours in 100 m intervals")

p5+ theme_minimal()

Mt. Baker 3D Map

Following Milos Makes Maps video

#trying to make a 3d map (using coords from Google Maps)
baker_lat <- 48.782079

baker_long <- -121.803322


#baker_3d <- plot_3d_vista(lat=baker_lat,long=baker_long,radius=5000, overlay_detail = 14, elevation_detail=13, zscale=5, img_provider = 'OpenStreetMap', cache_dir = 'testing',theta=25, phi=25, zoom=0.5,
 #            windowsize =1200, solid=T, background='grey10')

Andy’s Tutorial

Plot

mtbDEM <- rast("Mt_Baker/data/mtbDEM.tif")

mtbDEM

## class       : SpatRaster 
## dimensions  : 409, 628, 1  (nrow, ncol, nlyr)
## resolution  : 5, 5  (x, y)
## extent      : 596210.7, 599350.7, 5411132, 5413177  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source      : mtbDEM.tif 
## name        :      elev 
## min value   :  826.4308 
## max value   : 1689.3694

#quickplot of the data

mtbDEM <- as.data.frame(mtbDEM, xy=T)

ggplot()+
  geom_raster(data=mtbDEM, aes(x, y, fill=elev))+
  scale_fill_gradientn(colors=terrain.colors(100))+
  labs(x="Easting [m]", y="Northing [m]")+
                coord_equal()

Flow

[1] Take .tif and run it through rast() function

[2] Turn that into data frame by running it through as.data.frame function

[3] Use ggplot and ggraster function to create a nice map

Tidyterra

Using tidyterra package, grab geom_spatraster function and make it easier!

For fun colors you can use ?hypso.colors

#install.packages("tidyterra")

mtbDEM <- rast("Mt_Baker/data/mtbDEM.tif") #you need data in raster form here to run through spatraster code below. 

ggplot()+
  geom_spatraster(data=mtbDEM)+
  scale_fill_hypso_c(palette="usgs-gswa2")

|

The difference between this and geom_rast is that geom_spatraster uses coordinates to plot the data, while geom_rast uses meters. You can use spatrast in exactly the same way if you add coord_sf() into your arguments.

ggplot()+
  geom_spatraster(data=mtbDEM)+
  scale_fill_hypso_c(palette="usgs-gswa2")+
  coord_sf(datum = 32610)

Hillshade

mtb_hill <- rast("Mt_Baker/data/mtbHill.tif")

mtb_hill

## class       : SpatRaster 
## dimensions  : 409, 628, 1  (nrow, ncol, nlyr)
## resolution  : 5, 5  (x, y)
## extent      : 596210.7, 599350.7, 5411132, 5413177  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source      : mtbHill.tif 
## name        :    value 
## min value   : 0.000000 
## max value   : 0.999856

#make a quickplot

ggplot()+
  geom_spatraster(data=mtb_hill)

## oooo pretty!

Combined

p1 <- ggplot()+
  geom_spatraster(data=mtb_hill)+
  scale_fill_gradientn(colors=gray.colors(100, start=0.1, end=0.9), guide="none")+
  new_scale_fill()+
  geom_spatraster(data=mtbDEM)+
  scale_fill_hypso_c(name="Elevation [m]", 
                     palette="usgs-gswa2", alpha = 0.6)+
  theme_minimal()

p1

##ooooo v nice

Important

[1] Add hillshade in gray first

[2] Add elevation second with transparency

[3] Add multiple color scales (like gray and USGS) with new_scale_fill function

Shapefile

lines <- st_read("Mt_Baker/data/mtbChairLines.shp")

## Reading layer `mtbChairLines' from data source 
##   `C:\Users\jessi\Desktop\Rdata\ESCI_505\Mt_Baker\data\mtbChairLines.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: 596713.7 ymin: 5411451 xmax: 599131.5 ymax: 5412970
## Projected CRS: UTM_Zone_10_Northern_Hemisphere

p2 <- p1+geom_sf(data=lines)
p2

.csv

(See datacarpentry .csv to sf object section)

points <- st_read("Mt_Baker/data/mtbLodges.csv")

## Reading layer `mtbLodges' from data source 
##   `C:\Users\jessi\Desktop\Rdata\ESCI_505\Mt_Baker\data\mtbLodges.csv' 
##   using driver `CSV'

## Warning: no simple feature geometries present: returning a data.frame or tbl_df

points <- as.data.frame(points, xy=T)

names(points)

## [1] "X"  "Y"  "id"

lodges <- st_as_sf(points, coords = c("X", "Y"), crs = st_crs(32610))


st_crs(lodges)

## Coordinate Reference System:
##   User input: EPSG:32610 
##   wkt:
## PROJCRS["WGS 84 / UTM zone 10N",
##     BASEGEOGCRS["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]],
##         ID["EPSG",4326]],
##     CONVERSION["UTM zone 10N",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-123,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Navigation and medium accuracy spatial referencing."],
##         AREA["Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)."],
##         BBOX[0,-126,84,-120]],
##     ID["EPSG",32610]]

Multi-vector

p2+
  geom_sf(data=lodges)+
  ggtitle("Hi")