Update HW1

[ryoliver]

• [x] create qmd for assignment 1 on website
• [x] update assignment with Easter Island example
• [x] create template repository

[julietcohen]

Sub-tasks for "update assignment with Easter Island example":

• [x] use NOAA's marmap API to download bathymetry data around Easter Island, see if it looks good and is easy enough for students to work with, consider its CRS and resolution
• [ ] find 2 other vector datasets that would work with Easter Island example, allowing students to choose between data layers as they improve Easter Island plot as part of their homework
  • [X] seamounts
  • one more?

[julietcohen]

Bathymetry Data

This data can be imported using the marmap API easily without a key. After some exploration of the different marmap functions, the NOAA data is easiest to work with sf in raster format. Therefore I suggest that it is downloaded as a raster, then if we also want to provide the same data to the students in vector format, the raster should be converted to a dataframe using terra::as.data.frame

Raster bathymetry data added to Easter Island map

Vector bathymetry data added to Easter Island map

Some ideas:

• To produce these maps, I converted the bathymetry data to match the CRS of the Easter Island elevation raster, which is EPSG 32712 for UTM zone 12S. We could provide the data in a different CRS for the students to force them to convert it first
• We could provide a bathymetry raster with a much larger extent than the one I downloaded, meaning the students would have to find a way to make the extent smaller to retain an appropriate zoom level for Easter Island and the other data layers
• the resolution of the imported grid is in "minutes" (see resolution parameter in the documentation), we could ask the students to investigate or explain
• the units of the cell values is hard to find. I looked through the documentation and a publication about the package but ended up including meters in the legend units as an assumption
• I can add my qmd that contains code for the marmap API and these maps to the shared data dir in the drive so other can change/replace data files in the future as well (and I can convert the qmd to an R script if that's preferred)

[julietcohen]

Another consideration for mapping the bathymetry data:

The way I mapped both the raster and vector bathymetry data with tmap, including the parameter midpoint = NA within tm_dots() is important for visualizing all the values of the blue palette, because this overrides the default for this parameter, which is 0 when negative and positive values are present in the data. Setting this to NA instead allows for the full range of my palette to be visualized. If this parameter is not included, the bathymetry data may look like just 1- 2 colors. I can see how students may overlook this if they 1) do not have an idea of what the bathymetry data "should" look like with a full spectrum, or 2) they are unfamiliar with this parameter from lecture or reading.

[ryoliver]

This is great!! I'll start working on writing HW1. Could you send the code that you used to produce these maps so I can figure out which/how many hints to include?

[julietcohen]

I converted my qmd for the API and plotting to an R script and uploaded it to the shared data folder in a new subdir week1/bathymetry. The R script saves the raster and vector bathymetry files as a tif and gpkg, and outputs the plots to a single PDF. I uploaded these 3 files to the subdir as well so you don't have to run the script unless you want to re-do the data download or plots.

If we want to make any of the suggested changes to the files, like change the CRS or extent, I can take care of it!

[julietcohen]

Here's the R script here for convenience as well:

ei_bathy_tmap.R

```{r} # EDS 223 HW 1 Supplement: # Add NOAA Bathymetry Data to Easter Island tmap # Recreate the Easter Island plot from this figure: # https://r-tmap.github.io/tmap-book/nutshell.html#fig:rmap1 # and overlay NOAA bathymetry layer library(terra) library(tidyterra) library(tmap) library(sf) library(RColorBrewer) library(marmap) # customize name of output PDF pdf(file = "ei_bathy_maps.pdf") # read in layers for original EI tmap ei_elev = terra::rast("data/easter_island/ei_elev.tif") ei_borders = read_sf("data/easter_island/ei_border.gpkg") ei_roads = read_sf("data/easter_island/ei_roads.gpkg") ei_points = read_sf("data/easter_island/ei_points.gpkg") volcanoes = subset(ei_points, type == "volcano") # read in EI bathy data from marmap API ei_bathy <- marmap::getNOAA.bathy( lon1 = -109.582827, lon2 = -109.142687, lat1 = -27.230938, lat2 = -26.995631, resolution = 0.01 # units = "minutes" ) %>% marmap::as.raster() %>% # convert to spatRaster class terra::rast() # reproject CRS to match the elevation raster ei_bathy_transformed <- terra::project(ei_bathy, ei_elev) # save raster as tif file terra::writeRaster(ei_bathy_transformed, "ei_bathy.tif", overwrite = TRUE) # add bathymetry raster data to EI tmap blues_palette <- colorRampPalette(c("darkblue", "blue", "royalblue", "lightskyblue", "#DDDDDD"))(1000) title <- "Easter Island Topography, Volcanoes, Roads, and Bathymetry" tm_shape(ei_bathy_transformed) + tm_raster(style = "cont",, palette = blues_palette, title = "Depth from\nsea level (meters)", midpoint = NA) + # land elevation raster tm_shape(ei_elev) + tm_raster(style = "cont", palette = "-RdYlGn", title = "Elevation (m asl)") + # Easter Island perimeter vector tm_shape(ei_borders) + tm_borders(col = "black", lwd = 1.5) + # roads vector tm_shape(ei_roads) + tm_lines(col = "grey40", lwd = "strokelwd", legend.lwd.show = FALSE) + # volcano vector tm_shape(volcanoes) + tm_symbols(shape = 24, size = "elevation", col = "red", border.col = "black", title.size = "Volcanoes (m asl)") + tm_layout(main.title = title, main.title.size = 1, bg.color = "grey", legend.outside = TRUE, legend.outside.position = "right", frame = FALSE) # convert to dataframe ei_bathy_df <- terra::as.data.frame(ei_bathy_transformed, cells = FALSE, xy = TRUE, na.rm = TRUE) # combine x and y cols to a single point geom col ei_bathy_points <- st_as_sf(ei_bathy_df, coords = c("x", "y"), crs = st_crs(ei_elev)) %>% rename("depth" = "layer") # save vector as geopackage file sf::st_write(ei_bathy_points, "ei_bathy.gpkg", append = FALSE) # add bathymetry vector data to EI tmap tm_shape(ei_bathy_points) + tm_dots(col = "depth", size = 1, title = "Depth from\nsea level (meters)", palette = blues_palette, midpoint = NA) + # land elevation raster tm_shape(ei_elev) + tm_raster(style = "cont", palette = "-RdYlGn", title = "Elevation (m asl)") + # Easter Island perimeter vector tm_shape(ei_borders) + tm_borders(col = "black", lwd = 1.5) + # roads vector tm_shape(ei_roads) + tm_lines(col = "grey40", lwd = "strokelwd", legend.lwd.show = FALSE) + # volcano vector tm_shape(volcanoes) + tm_symbols(shape = 24, size = "elevation", col = "red", border.col = "black", title.size = "Volcanoes (m asl)") + tm_layout(main.title = title, main.title.size = 1, bg.color = "grey", legend.outside = TRUE, legend.outside.position = "right", frame = FALSE) print("Script complete. Plots saved to file 'ei_bathy_maps.PDF'") ```

[julietcohen]

In pursuit of another vector dataset option for HW1, I downloaded this seamount dataset and subset the points in R to the same bounding box used to subset the bathymetry data. Here they're visualized with mapview:

It's only 4 points and arguably not that interesting, but has a few attributes that could be used for the palette (height, depth, etc) and aligns well with the bathymetry dataset.

Just need tmap to zoom a little bit out to gather that 4th point, I think...

[ryoliver]

thanks so much! this is all great! working on putting together HW1 now.

[julietcohen]

Great, in case we want to use this data, I'll clean up the code for subsetting it/tranforming etc. and share it.

I added the data to the folder in week1/seamounts. global_seamounts contains the original shapefile file downloaded from the source with every seamount documented, and ei_seamounts is just the 4 relevant points for Easter Island.

[ryoliver]

that's perfect, thanks!

[julietcohen]

I wonder if it be confusing for students to be introduced to geopackage format so early in the quarter? It requires no extra steps, it's imported the same as a shapefile, but I don't want to frazzle them.

[julietcohen]

My code for subsetting the seamount data from the global file to just the Easter Island seamounts, as well as the code to plot all the data layers together using tmap, is now in an R script in the data dir here, called ei_sm_tmap.R

I was able to make the plot include all seamount points by creating a new bbox for the tmap, as shown in the script. The plot is output to a PDF in the same subdir called ei_all_layers.pdf

For convenience, the script is here as well:

ei_sm_tmap.R

```{r} # EDS 223 HW 1 Supplement: # Add Seamount Data to Easter Island tmap # Recreate the Easter Island plot from this figure: # https://r-tmap.github.io/tmap-book/nutshell.html#fig:rmap1 # and overlay NOAA bathymetry layer library(terra) library(tmap) library(sf) library(RColorBrewer) # customize name of output PDF pdf(file = "ei_all_layers.pdf") # read in layers for original EI tmap ei_elev = terra::rast("data/week1/easter_island/ei_elev.tif") ei_borders = read_sf("data/week1/easter_island/ei_border.gpkg") ei_roads = read_sf("data/week1/easter_island/ei_roads.gpkg") ei_points = read_sf("data/week1/easter_island/ei_points.gpkg") volcanoes = subset(ei_points, type == "volcano") ei_bathy = terra::rast("data/week1/bathymetry/ei_bathy.tif") # read in all seamount points in world sm <- st_read("data/week1/seamounts/global_seamounts/Seamounts.shp") # transform to CRS for EI sm_transformed <- st_transform(sm, st_crs(ei_elev)) # define bbox with coords from bboxfinder.com # (same coords used for marmap API) bbox <- st_bbox(c(xmin = -109.582827, ymin = -27.230938, xmax = -109.142687, ymax = -26.995631), crs = 4326) # convert bbox to an sfs obj bc it is of # class bbox, then convert to sf obj # note: cannot convert class bbox straight to sf obj # bc no geometry col sf_bbox <- st_as_sfc(bbox) |> st_sf() # project to ei data CRS bbox_transformed <- st_transform(sf_bbox, st_crs(ei_elev)) # subset seamounts to just bbox (EI) ei_sm <- st_join(sm_transformed, bbox_transformed, join = st_within, left = FALSE) # alternatively to all above processing code, just # read in seamount file saved to data dir: # ei_sm <- st_read("data/week1/seamounts/ei_seamounts.gpkg") # For tmap with ALL data layers: # define larger plot bbox that is buffered around both of # the two largest layers to display all 4 seamounts in view bbox_sm <- st_bbox(ei_sm) bbox_bathy <- st_bbox(ei_bathy) bbox_largest <- st_bbox(c(xmin = min(bbox_bathy[1], bbox_sm[1]), ymin = min(bbox_bathy[2], bbox_sm[2]), xmax = max(bbox_bathy[3], bbox_sm[3]), ymax = max(bbox_bathy[4], bbox_sm[4]))) blues <- c("darkblue", "blue", "royalblue", "lightskyblue", "#DDDDDD") blues_palette <- colorRampPalette(blues)(1000) purples <- c("purple", "orchid", "#CF9FFF", "lightpink") purples_palette <- colorRampPalette(purples)(4) title = "Easter Island Topography, Volcanoes,\nRoads, Bathymetry, and Seamounts" tm_shape(ei_bathy, bbox = bbox_largest) + tm_raster(style = "cont", palette = blues_palette, title = "Depth from\nsea level (meters)", midpoint = NA) + tm_shape(ei_sm) + tm_dots(col = "HEIGHT", size = 0.5, title = "Height (meters)", palette = purples_palette, midpoint = NA) + tm_shape(ei_elev) + tm_raster(style = "cont", palette = "-RdYlGn", title = "Elevation (m asl)") + tm_shape(ei_borders) + tm_borders(col = "black", lwd = 1.5) + tm_shape(ei_roads) + tm_lines(col = "#DDDDDD", lwd = "strokelwd", legend.lwd.show = FALSE) + tm_shape(volcanoes) + tm_symbols(shape = 24, size = "elevation", col = "red", border.col = "black", title.size = "Volcanoes (m asl)") + tm_layout(main.title = title, main.title.size = 1, bg.color = "#DDDDDD", legend.outside = TRUE, legend.outside.position = "right", frame = FALSE) ```