1

I looked around and I have not found a nice solution for my goal. I want to plot some data on a longitude/ latitude plot using ggplot2 and the coastline plus bathymetry with marmap, everything in one single plot.

This script is to plot mydata

ggplot(data = ctd, aes(x = Longitude, y = Latitude)) +
  geom_raster(aes(fill = Temp)) +
  scale_fill_gradientn(colours = rev(my_colours)) +
  geom_contour(aes(z = Temp), binwidth = 2, colour = "black", alpha = 0.2) +

  #plot stations locations
  geom_point(data = ctd, aes(x = Longitude, y = Latitude),
             colour = 'black', size = 3, alpha = 1, shape = 15) +

  #plot legends
      labs(y = "Latitude", x = "Longitude", fill = "Temp (°C)") +
      coord_cartesian(expand = 0)+
      ggtitle("Temperature distribution") 

Using marmap I download the bathymetry

library(marmap)
Bathy <- getNOAA.bathy(lon1 = 37, lon2 = 38.7,
                       lat1 = -45.5, lat2 = -47.3, resolution = 1)

The result I would like to obtain is the distribution of mydata on Lon/Lat with the land colored in black plus grey lines for the bathymetry.

Strobila
  • 317
  • 3
  • 15

2 Answers2

2

Well, there's a marmap function for that. It's called autoplot.bathy(). Have you checked it's help file?

library(marmap) ; library(ggplot2)

library(marmap)
Bathy <- getNOAA.bathy(lon1 = 37, lon2 = 38.7,
                       lat1 = -45.5, lat2 = -47.3, resolution = 1)

ctd <- data.frame(Longitude = c(37.5, 38, 38.5), Latitude = c(-47, -46.5, -46))

autoplot.bathy(Bathy, geom=c("tile","contour")) +
    scale_fill_gradient2(low="dodgerblue4", mid="gainsboro", high="darkgreen") +
    geom_point(data = ctd, aes(x = Longitude, y = Latitude),
               colour = 'black', size = 3, alpha = 1, shape = 15) +
    labs(y = "Latitude", x = "Longitude", fill = "Elevation") +
    coord_cartesian(expand = 0)+
    ggtitle("A marmap map with ggplot2") 

enter image description here

Or, with base graphics (and the correct aspect ratio):

# Creating color palettes
blues <- c("lightsteelblue4", "lightsteelblue3", "lightsteelblue2", "lightsteelblue1")
greys <- c(grey(0.6), grey(0.93), grey(0.99))

# Plot
plot(Bathy, image = TRUE, land = TRUE, n=30, lwd = 0.1, bpal = list(c(0, max(Bathy), greys), c(min(Bathy), 0, blues)), drawlabels = TRUE)

# Add coastline
plot(Bathy, deep = 0, shallow = 0, step = 0, lwd=2, add = TRUE)

# Add stations
points(ctd, pch=15, cex=1.5)

enter image description here

Benoit
  • 1,154
  • 8
  • 11
  • I really like the `ggplot2` approach with the `autoplot.bathy()` function. However, I think you lose control of the isobath. How would you be able to highlight the (let's say) -300 isobath using this approach? Many thanks in advance! – Antonio Canepa Aug 31 '23 at 08:30
  • Sure, as always with R, you can do it all manually. Take a look [here](https://stackoverflow.com/questions/73582299/creating-bathymetric-lines-from-getnoaa-bathy-and-add-them-to-ggplot-how-do-i-h/73585982?noredirect=1#comment135766501_73585982) for some example code. It is much more complex than the `autoplot.bathy()` approach, but you'll get all the control you need. – Benoit Sep 01 '23 at 09:41
1

Here is an approach:

get bathy data:

library(marmap)
Bathy <- getNOAA.bathy(lon1 = 37, lon2 = 38.7,
                       lat1 = -45.5, lat2 = -47.3, resolution = 1)

convert it to matrix:

Bathy <- as.matrix(Bathy)
class(Bathy) <- "matrix"

now reshape it to long format and plot

library(tidyverse)

Bathy %>%
  as.data.frame() %>%
  rownames_to_column(var = "lon") %>%
  gather(lat, value, -1) %>%
  mutate_all(funs(as.numeric)) %>%
  ggplot()+
  geom_contour(aes(x = lon, y = lat, z = value), bins = 10, colour = "black") +
  coord_map()

enter image description here

missuse
  • 19,056
  • 3
  • 25
  • 47
  • This way is really good. I edited my question with more info. – Strobila May 03 '18 at 10:01
  • @Luca Stirnimann I suggest asking another question instead of editing an existing, already answered question. – missuse May 03 '18 at 18:31
  • @missuse I have seen some tutorials online that simply use region names as arguments, such as "hawaii" or "aleutians". Do you know where I can find a complete list of regions that are recognised? I'm looking for eastern Mediterranean around Israel. – JJGabe Nov 22 '22 at 11:28