Coastline of California
Project summary: How can we classify bathymetry of coasts? Create depth categories to reclassify area depth into categories as defined by Marine Biology: A Very Short Introduction.
1. Imports the libraries and datasets.
options(stringsAsFactors = FALSE)
library(raster)
library(rgdal)Please note that rgdal will be retired during 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-6, (SVN revision 1201)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/gdal
GDAL binary built with GEOS: FALSE
Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-0
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.2. Read in the dataset.
This will take any GeoTIFF file on the NOAA coast website.
coast <- raster("c.tif")Examine their contents
coastclass : RasterLayer
dimensions : 5849, 1697, 9925753 (nrow, ncol, ncell)
resolution : 3.28, 3.28 (x, y)
extent : 7318267, 7323833, 818041.8, 837226.6 (xmin, xmax, ymin, ymax)
crs : +proj=lcc +lat_0=43.6666666666667 +lon_0=-120.5 +lat_1=46 +lat_2=44.3333333333333 +x_0=2500000.0001424 +y_0=0 +ellps=GRS80 +units=ft +no_defs
source : c.tif
names : c summary(coast)Warning: summary is an estimate based on a sample of 1e+05 cells (1.01% of all cells) c
Min. -22.61115
1st Qu. 4.63188
Median 31.52154
3rd Qu. 90.47370
Max. 406.44623
NA's 2704371.00000This prints an image of the rasterized area. Some metadata that was included about the dataset: - “Digital Elevation Model (DEM) files contain rasterized topobathy lidar elevations at a 1 m grid size” - “Horizontal positions, provided in decimal degrees of latitude and longitude, are referenced to the North American Datum of 1983 National Adjustment of 2011 (NAD83 NA11). Vertical positions are provided orthometric heights referenced to the North American Vertical Datum of 1988 (NAVD88)and provided in meters.” - The white refers to NoData values. LiDAR’s capability to scan around large rocks in the water is limited.
library(tidyverse)── Attaching core tidyverse packages ─────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.2 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.0
✔ lubridate 1.9.2 ✔ tibble 3.2.1
✔ purrr 1.0.1 ✔ tidyr 1.3.0── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
✖ dplyr::select() masks raster::select()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorscoast %>%
as.list %>%
lapply(spplot)[[1]]
This code describes how many pixels fit into different depth categories.
hist(coast,
main = "Depth distribution",
xlab = "Depth (m)", ylab = "Number of Pixels")Warning: 1% of the raster cells were used. 100000 values used.
We can estimate that range 0 to -10m has the most pixels.
Some more detailed information about the raster:
histinfo <- hist(coast)Warning: 1% of the raster cells were used. 100000 values used.
histinfo$counts [1] 17 20198 22657 13418 11165 5124 3906 2991 3153 2698 2363 2389 2385 1929 1957 1257 892 581 463 249 115 77 16histinfo$breaks [1] -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 4203. Reclassification
According to “Marine Biology: A Very Short Introduction”, the continental shelf ends at 200m, the continental slope and seamounts at 1000m, the continental rise and abyssal plains at 5000m, and trenches at 11,000. The histogram below shows how many pixels in the raster are classified under each category. I did not include trenches.
hist(coast,
breaks = c(-6000, -3000, -1000, -200, 0, 500),
main = "Histogram with custom breaks",
xlab = "Depth (m)" , ylab = "Number of Pixels")Warning: 1% of the raster cells were used. 100000 values used.
This code defines the breaks of the reclassification matrix.
reclassdf <- c(-6000, -3000, -5, -3000, -1000, -4, -1000, -200, -3, -200, 0, -2, 0, 500, -1)
reclassdf [1] -6000 -3000 -5 -3000 -1000 -4 -1000 -200 -3 -200 0 -2 0 500 -1This code “reshapes the object into a matrix with columns and rows”.
reclass_m <- matrix(reclassdf,
ncol = 3,
byrow = TRUE)
reclass_m [,1] [,2] [,3]
[1,] -6000 -3000 -5
[2,] -3000 -1000 -4
[3,] -1000 -200 -3
[4,] -200 0 -2
[5,] 0 500 -1This code reclassifies the image according to the classes set by the book.
coast_classified <- reclassify(coast,
reclass_m)Shows the pixel distribution.
# view reclassified data
barplot(coast_classified,
main = "Number of pixels in each class")Warning: a sample of 10.1% of the raster cells were used to estimate frequencies
-2 represents continental shelf pixels. -1 represents terrestrial pixels.
plot(coast_classified)
4. Making the map into a neater ggplot map.
library(stars)Loading required package: abind
Loading required package: sf
Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE#--- converting to a data.frame ---#
coast_df <- as.data.frame(coast_classified, xy = TRUE) %>%
na.omit()
#--- take a look ---#
head(coast_df)