Bathymetry data

[pennybeaver]

Hi Martijn, I am trying to work out how to add bathymetry lines to a tmap, I have tried converting a marmot bathymetry map of the area I specifically would like to use, as.raster but I am unable to use so far in tmap. The maps are for a species of shearwater I am researching showing their habitat use. Can I add bathymetry data to tmap?

or is there a better way to do it in tmap by not using marmap, I am getting a lot of errors and wandering if tmap and marmap are not compatible?

Kind regards
Penny
tmap wedgies.pdf
marmap1-bathy.pdf

#Create the bathymetry map using marmap

library(marmap)
papoue <- getNOAA.bathy(lon1 = 100, lon2 = 180, lat1 = -60, lat2 = 40, resolution = 15)
plot(papoue, land = TRUE, deep = c(-9000, -3000, 0),
     shallow = c(-3000, -10, 0),
     step = c(1000, 1000, 0),
     lwd = c(0.8, 0.8, 1), lty = c(1, 1, 1),
     col = c("lightgrey", "darkgrey","black"),
     drawlabel = c(FALSE, FALSE, FALSE))

bathy<-marmap::as.raster(papoue)
class(bathy)

#Create a map using data from wedgetailed shearwaters from Muttonbid Island

library(tmap)
library(tmaptools)
data(World)
tmap_mode("plot")
library(raster)

#add it to tmap maps?
tmap::tm_shape(mtnr)+tm_raster(palette = "YlGnBu", legend.show = FALSE)+tm_shape(World)+tm_borders()+
  tm_grid(col ="gray80", alpha = 0.3)+tm_fill(col = "black")+
  tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+
  tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom'))+
  tm_shape(siteMT)+tm_dots(size = 0.2, shape = 8, col = "red")+
  tm_shape(siteMG)+tm_dots(size = 0.2, shape = 8, col = "green")+
  tm_shape(siteGI)+tm_dots(size = 0.2, shape = 8, col = "green")

[mtennekes]

Hi Penny,

What kind errors do you get?

bathy is a raster object, so you can plot with tm_raster (as you already did). Make sure the palette is a diverging palette (since there are positive and negative values). E.g.:

tmap::tm_shape(bathy)+tm_raster("layer", palette = "-RdBu", legend.show = T, style = "cont")+
tm_shape(World)+tm_borders()+
	tm_grid(col ="gray80", alpha = 0.3)+tm_fill(col = "black")+
	tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+
	tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom'))

The plot function from marmap creates contour lines with the base function contour. The package stars offers a good alternative, st_contour (which is I guess a wrapper around a function from GDAL):

library(stars)
b = st_as_stars(bathy)
x = st_contour(b, breaks = seq(-10000, 5000, by = 1000), contour_lines = TRUE)

tmap::tm_shape(x)+
	tm_lines("layer", palette = pals::kovesi.rainbow(15), lwd = 2, legend.show = FALSE, 
        title.col = "Elevation")+tm_shape(World)+tm_borders()+
	tm_grid(col ="gray80", alpha = 0.3)+
	tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+
	tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,
        title.position = c('left','bottom'), legend.outside = TRUE)

(also pinging @Nowosad since it is proper educational material)

[pennybeaver]

Hi Martijn Thank you for responding and providing some awesome code, big thanks. Based on what you have provided clearly marmap and tmap have similar functions hence the errors, I all need to add was tmap:: so thank you. The crowd goes wild as a map is produced 🙂 woohooo. Sorry to ask but is there a function to add country and ocean names in tmap or is it another layer? [cid:a708ac76-d870-452b-9956-a47d89f28358]

…

________________________________ From: mtennekes ***@***.***> Sent: Monday, 12 April 2021 6:11 PM To: mtennekes/tmap-workshop ***@***.***> Cc: Penny Beaver ***@***.***>; Author ***@***.***> Subject: Re: [mtennekes/tmap-workshop] Bathymetry data (#1) Hi Penny, What kind errors do you get? bathy is a raster object, so you can plot with tm_raster (as you already did). Make sure the palette is a diverging palette (since there are positive and negative values). E.g.: tmap::tm_shape(bathy)+tm_raster("layer", palette = "-RdBu", legend.show<http://legend.show> = T, style = "cont")+ tm_shape(World)+tm_borders()+ tm_grid(col ="gray80", alpha = 0.3)+tm_fill(col = "black")+ tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+ tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom')) [image]<https://user-images.githubusercontent.com/2444081/114357421-01943800-9b72-11eb-85a8-9d05399c4d84.png> The plot function from marmap creates contour lines with the base function contour. The package stars offers a good alternative, st_contour (which is I guess a wrapper around a function from GDAL): library(stars) b = st_as_stars(bathy) x = st_contour(b, breaks = seq(-10000, 5000, by = 1000), contour_lines = TRUE) tmap::tm_shape(x)+ tm_lines("layer", palette = pals::kovesi.rainbow(15), lwd = 2, legend.show<http://legend.show> = FALSE, title.col = "Elevation")+tm_shape(World)+tm_borders()+ tm_grid(col ="gray80", alpha = 0.3)+ tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+ tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8, title.position = c('left','bottom'), legend.outside = TRUE) [image]<https://user-images.githubusercontent.com/2444081/114362014-250db180-9b77-11eb-8aa0-bd9d4d09cafe.png> (also pining @Nowosad<https://github.com/Nowosad> since it is proper educational material) — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFIRKDYNVWCF7NJDCJSKELTTIKTKLANCNFSM42UEL6EA>. This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.

[mtennekes]

Glad I could help:-)

Spatial objects (polygons or points) are required to add country and ocean names. For countries, you can use tmap's World, which you already used for the country borders.
I don't know where to get ocean data. But since there are only 2 oceans, I've created a spatial object myself. See below the (cleaned) code.

# Coordinates of the ocean labels (which can be changed if needed)

library(sf)
o = st_sf(name = c("Indian Ocean", "Pacific Ocean"), st_sfc(st_point(c(110, -50)), st_point(c(160, 20))), crs = 4326)

## plot 1 (raster data)
tm_shape(bathy)+tm_raster("layer", palette = "-RdBu", legend.show = FALSE, style = "cont")+
tm_shape(World)+
	tm_polygons(col = "black")+
	tm_text("name", bg.color = "black", bg.alpha = 0.6) +
tm_shape(o) +
	tm_text("name") +
tm_grid(col ="gray80", alpha = 0.3)+
tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+
tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom'))

## plot 2 (contour lines)
tm_shape(x)+
    tm_lines("layer", palette = pals::kovesi.rainbow(15), lwd = 2, legend.show = FALSE, title.col = "Elevation")+
tm_shape(World)+
    tm_borders()+ 
	tm_text("name") +
tm_shape(o) +
	tm_text("name") +
tm_grid(col ="gray80", alpha = 0.3)+
tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+
tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom'), legend.outside = TRUE)

[pennybeaver]

very cool thanks Martijn I think I have landed with the raster layer, it allows the migratory pathway to be much clearer. As you can see the below map shows a lot more detail with regards to the oceans this species uses, ocean troughs if you know where they are located (Mariana and Philippines trough) in relation to land masses for this migratory seabird called a wedge-tailed shearwater :) now back to the write up zzzzzzzzz.......zzzz not as exciting as producing maps! [cid:e8cecdda-bfe9-4af4-8b2f-8d5aa6178627]

…

________________________________ From: mtennekes ***@***.***> Sent: Monday, 12 April 2021 10:17 PM To: mtennekes/tmap-workshop ***@***.***> Cc: Penny Beaver ***@***.***>; Author ***@***.***> Subject: Re: [mtennekes/tmap-workshop] Bathymetry data (#1) Glad I could help:-) Spatial objects (polygons or points) are required to add country and ocean names. For countries, you can use tmap's World, which you already used for the country borders. I don't know where to get ocean data. But since there are only 2 oceans, I've created a spatial object myself. See below the (cleaned) code. # Coordinates of the ocean labels (which can be changed if needed) library(sf) o = st_sf(name = c("Indian Ocean", "Pacific Ocean"), st_sfc(st_point(c(110, -50)), st_point(c(160, 20))), crs = 4326) ## plot 1 (raster data) tm_shape(bathy)+tm_raster("layer", palette = "-RdBu", legend.show<http://legend.show> = FALSE, style = "cont")+ tm_shape(World)+ tm_polygons(col = "black")+ tm_text("name", bg.color = "black", bg.alpha = 0.6) + tm_shape(o) + tm_text("name") + tm_grid(col ="gray80", alpha = 0.3)+ tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+ tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom')) ## plot 2 (contour lines) tm_shape(x)+ tm_lines("layer", palette = pals::kovesi.rainbow(15), lwd = 2, legend.show<http://legend.show> = FALSE, title.col = "Elevation")+ tm_shape(World)+ tm_borders()+ tm_text("name") + tm_shape(o) + tm_text("name") + tm_grid(col ="gray80", alpha = 0.3)+ tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+ tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,title.position = c('left','bottom'), legend.outside = TRUE) [image]<https://user-images.githubusercontent.com/2444081/114392865-a5450e80-9b99-11eb-8909-12ca89269d95.png> [image]<https://user-images.githubusercontent.com/2444081/114392895-aece7680-9b99-11eb-8281-ae0d61323ef0.png> — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub<#1 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFIRKDYIQS7TRSSOSGTO5OTTILQEZANCNFSM42UEL6EA>. This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.

[AnaLauraMachado]

Hi ,

I just found this and it is just what I was looking for to work with my bathymetry data.
When I run the following line I get an error. Can you help me to fix it?

b = st_as_stars(Bathy)
x = st_contour(b, breaks = seq(-1000, 200, by = 1000), contour_lines = TRUE)
Error in Ops.units(max(breaks), max(x[[1]], na.rm = TRUE)) :
both operands of the expression should be "units" objects

Thank you.

[Nowosad]

@AnaLauraMachado -- can you provide a reproducible example (similar to the one below)?

library(marmap)
library(tmap)
library(stars)
data(World)
papoue = getNOAA.bathy(lon1 = 100, lon2 = 180, lat1 = -60, lat2 = 40, 
                        resolution = 15)
bathy = marmap::as.raster(papoue)
b = st_as_stars(bathy)
x = st_contour(b, breaks = seq(-10000, 5000, by = 1000), contour_lines = TRUE)

tm_shape(x) +
  tm_lines("layer", palette = pals::kovesi.rainbow(15), lwd = 2, 
           legend.show = FALSE, title.col = "Elevation") +
  tm_shape(World) +
  tm_borders() +
  tm_grid(col = "gray80", alpha = 0.3)+
  tm_xlab("Longitude", size = 0.8) + 
  tm_ylab("Latitude", size = 0.8) +
  tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8,
            title.position = c("left", "bottom"), legend.outside = TRUE)

[AnaLauraMachado]

Hi, thanks for the quick reply. library(tmap) library(marmap) library(stars) Data(Site) Bathy <- raster("gebco_2020_n-60.0_s-64.0_w-61.0_e-57.0.tif") Fildes <- st_read ("C:/Users/almac/Desktop/PhD/shp/ADD_Coastline_high_res_polygon_Sliced.shp") st_crs(Fildes) = 3031 #Quantarctica usa proyeccion esterografica Fildes_WGS84 <- st_transform (Fildes, 4326) tmap_mode("plot") tmap_options(check.and.fix = TRUE) bb <-st_bbox(proj_trips) map_fildes <- tm_shape(Fildes,bbox = bb)+ #c(-63, -64, -57, -62)) + tm_borders(lwd=0.3, col = "black")+ tm_fill("grey") map_pen <- tm_shape(Bathy)+ tm_raster(style= "cont",palette="-RdBu", breaks = seq(-2000, 100, by = 200))+ tm_layout(legend.outside = TRUE)+ tm_shape(Site,is.master = TRUE) + tm_raster(alpha = .8,"N_animals", style="cont", palette = "-magma", colorNA=NULL, colorNULL=NULL)+ tm_graticules(col = "grey", lwd = 0.2, labels.size = 1)+ tm_layout(bg.color = "lightblue") + map_fildes map_pen b = st_as_stars(Bathy) x = st_contour(b, breaks = seq(-2000, 100, by = 200), contour_lines = TRUE) Error in Ops.units(max(breaks), max(x[[1]], na.rm = TRUE)) : both operands of the expression should be "units" objects El vie, 3 mar 2023 a las 0:37, Jakub Nowosad ***@***.***>) escribió:

…

@AnaLauraMachado <https://github.com/AnaLauraMachado> -- can you provide a reproducible example (similar to the one below)? library(marmap) library(tmap) library(stars) data(World)papoue = getNOAA.bathy(lon1 = 100, lon2 = 180, lat1 = -60, lat2 = 40, resolution = 15)bathy = marmap::as.raster(papoue)b = st_as_stars(bathy)x = st_contour(b, breaks = seq(-10000, 5000, by = 1000), contour_lines = TRUE) tm_shape(x) + tm_lines("layer", palette = pals::kovesi.rainbow(15), lwd = 2, legend.show = FALSE, title.col = "Elevation") + tm_shape(World) + tm_borders() + tm_grid(col = "gray80", alpha = 0.3)+ tm_xlab("Longitude", size = 0.8) + tm_ylab("Latitude", size = 0.8) + tm_layout(title = "(a) n = 19", scale = 0.7, main.title.size = 0.8, title.position = c("left", "bottom"), legend.outside = TRUE) <https://camo.githubusercontent.com/04ccbb3f0d940ea7c608b36e854d448f9615c2743c5b36066579fa6c4e328c0e/68747470733a2f2f692e696d6775722e636f6d2f525776474661742e706e67> — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANABFLLJUVAGC7L7RFVQ3FDW2GUUZANCNFSM42UEL6EA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[Nowosad]

Could you share the gebco_2020_n-60.0_s-64.0_w-61.0_e-57.0.tif data with us?

[AnaLauraMachado]

El dom, 5 mar 2023 a las 9:11, Jakub Nowosad ***@***.***>) escribió:

…

Could you share the gebco_2020_n-60.0_s-64.0_w-61.0_e-57.0.tif data with us? — Reply to this email directly, view it on GitHub <#1 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ANABFLJK7GJWDNDBVO74CSLW2R7HRANCNFSM42UEL6EA> . You are receiving this because you were mentioned.Message ID: ***@***.***>