Read the Tobias Depth Elevation Model

Setup

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.3     ✓ purrr   0.3.4
## ✓ tibble  3.1.2     ✓ dplyr   1.0.6
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   1.4.0     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sf)

## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1

library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

Setup

library(tidyverse)
library(sf)
library(raster)

Read the Data

Since this is a raster dataset, we use raster().

DEM = raster("DEM_SF.tif")
crs(DEM)

## CRS arguments: +proj=longlat +datum=NAD83 +no_defs

plot(DEM)

Convert the CRS

Since this data is not using the same CRS as the rest of the data, we need to convert it. To do this, we need an sf object with the CRS we want. We can use the shoreline sf object.

shoreline = st_read("Shoreline.shp")

## Reading layer `Shoreline' from data source `/Users/haroldnelson/Dropbox/QGIS/Tobias/Shoreline.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -123.1086 ymin: 37.69325 xmax: -122.3277 ymax: 37.92975
## Geodetic CRS:  WGS84(DD)

DEM_tr = projectRaster(DEM,crs = crs(shoreline))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in Proj4
## definition

## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum WGS84 in Proj4 definition

crs(DEM_tr)

## CRS arguments: +proj=longlat +ellps=WGS84 +no_defs

Convert to Dataframe

DEM_tr_df = as.data.frame(DEM_tr, xy = TRUE)
glimpse(DEM_tr_df)

## Rows: 10,331,100
## Columns: 3
## $ x      <dbl> -122.6111, -122.6110, -122.6109, -122.6108, -122.6107, -122.610…
## $ y      <dbl> 37.89356, 37.89356, 37.89356, 37.89356, 37.89356, 37.89356, 37.…
## $ DEM_SF <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…

Plot

DEM_tr_df %>% 
ggplot() +
     geom_raster(aes(x = x, y = y, 
                  fill = DEM_SF))