jamesgrecian · June 12, 2019 17:35
diff --git a/Generate map wedge in ggplot.r b/Generate map wedge in ggplot.r
 ###########################################################################################
 ### Example script to generate a map using ggplot2 with similar look to Ocean Data View ###
 ###########################################################################################

 # The example map I am trying to recreate is here:
 # https://twitter.com/JamesGrecian/status/1137018574504153089?s=20

 # Load in libraries
 require(marmap) # for bathymetry data
 require(tidyverse)
 require(sf)
 # devtools::install_github("jamesgrecian/mapr", force = T)
 require(mapr) # my wrapper for shapefile generation

 # Some helper functions for colour ramps from Claus Wilke
 # https://github.com/tidyverse/ggplot2/issues/2673
 discrete_gradient_pal <- function(colours, bins = 11) {
  ramp <- scales::colour_ramp(colours)
  
  function(x) {
    if (length(x) == 0) return(character())
    
    i <- floor(x * bins)
    i <- ifelse(i > bins-1, bins-1, i)
    ramp(i/(bins-1))
  }
 }

 scale_colour_discrete_gradient <- function(..., colours, bins = 11, na.value = "grey50", guide = "colourbar", aesthetics = c("colour", "fill"), colors)  {
  colours <- if (missing(colours)) 
    colors
  else colours
  continuous_scale(
    aesthetics,
    "discrete_gradient",
    discrete_gradient_pal(colours, bins),
    na.value = na.value,
    guide = guide,
    ...
  )
 }

 # Extract depth data from NOAA
 # This may take some time...
 batdat <- marmap::getNOAA.bathy(lon1 = -80,
                                lon2 = 20,
                                lat1 = -40,
                                lat2 = -90,
                                resolution = 2) 
 bat = marmap::as.raster(batdat)
 raster::projection(bat) = sp::CRS("+proj=longlat +datum=WGS84")

 # Define a projection and transform bathymetry data
 prj <- '+proj=laea +lat_0=-60 +lon_0=-40 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
 r = raster::projectRaster(bat, crs = sp::CRS(prj), method = 'ngb', res = 50000)
 bat_df = as.data.frame(raster::rasterToPoints(r))
 names(bat_df) = c("x", "y", "z")

 # Example plot ggplot
 p1 <- ggplot() +
  theme_bw() +
  geom_point(aes(x = x, y = y, colour = z), data = bat_df) +
  geom_raster(aes(x = x, y = y, fill = z), data = bat_df) +
  scale_colour_discrete_gradient("Depth (m)",
                                 colours = rev(RColorBrewer::brewer.pal(11, "RdYlBu")),
                                 limits = c(-8000, 800),
                                 breaks = seq(-8000, 800, 800),
                                 guide = guide_colourbar(nbin = 500,
                                                         raster = T,
                                                         frame.colour = "black",
                                                         ticks.colour = "black",
                                                         frame.linewidth = 1,
                                                         barwidth = 1,
                                                         barheight = 25,
                                                         direction = "vertical",
                                                         title.position = "right",
                                                         title.theme = element_text(angle = 90,
                                                                                    hjust = 0.5))) +
  geom_sf(aes(), data = mapr::mapr(ellie, prj, buff = 1e7)) +
  coord_sf(xlim = c(-3200000, 4200000), ylim = c(-2800000, 2100000), crs = prj, expand = T) +
  xlab("") +ylab("")

 print(p1)

 # To make a wedge shaped map create a mask that matches map extent
 # Help again from Claus Wilke
 # https://htmlpreview.github.io/?https://github.com/clauswilke/practical_ggplot2/blob/master/goode.html
  
  
 # Define a clipping polygon based on bathymetry raster
 pol <- sf::st_polygon(list(rbind(c(-80, -90),
                                 c( 20, -90),
                                 c( 20, -40),
                                 c(-80, -40),
                                 c(-80, -90))))
 # Segmentise to ensure projection becomes a 'wedge' and not a trapezoid
 pol <- sf::st_segmentize(pol, 1)
 pol <- pol %>% sf::st_sfc() 
 pol <- pol %>% st_set_crs(4326)
 pol <- pol %>% sf::st_transform(prj)

 # get the bounding box in transformed coordinates and expand by 10%
 xlim <- st_bbox(pol)[c("xmin", "xmax")]*1.1
 ylim <- st_bbox(pol)[c("ymin", "ymax")]*1.1

 # turn into enclosing rectangle
 encl_rect <- 
  list(
    cbind(
      c(xlim[1], xlim[2], xlim[2], xlim[1], xlim[1]), 
      c(ylim[1], ylim[1], ylim[2], ylim[2], ylim[1])
    )
  ) %>%
  st_polygon() %>%
  st_sfc(crs = prj)

 # calculate the area outside the earth outline as the difference
 # between the enclosing rectangle and the earth outline
 cookie <- st_difference(encl_rect, pol)

 # Add this to the plot to mask out the areas not of interest
 # More effort required to move axis labels...

 p2 <- ggplot() +
  theme_minimal() +
  geom_point(aes(x = x, y = y, colour = z), data = bat_df) +
  geom_raster(aes(x = x, y = y, fill = z), data = bat_df) +
  scale_colour_discrete_gradient("Depth (m)",
                                 colours = rev(RColorBrewer::brewer.pal(11, "RdYlBu")),
                                 limits = c(-8000, 800),
                                 breaks = seq(-8000, 800, 800),
                                 guide = guide_colourbar(nbin = 500,
                                                         raster = T,
                                                         frame.colour = "black",
                                                         ticks.colour = "black",
                                                         frame.linewidth = 1,
                                                         barwidth = 1,
                                                         barheight = 25,
                                                         direction = "vertical",
                                                         title.position = "right",
                                                         title.theme = element_text(angle = 90,
                                                                                    hjust = 0.5))) +
  geom_sf(aes(), data = mapr::mapr(ellie, prj, buff = 1e7)) +
  geom_sf(fill = "white", color = "grey", data = cookie) +
  coord_sf(xlim = c(-3200000, 4200000), ylim = c(-2800000, 2100000), crs = prj, expand = T) +
  xlab("") +ylab("")

 gridExtra::grid.arrange(p1, p2, ncol = 2)
	###########################################################################################
	### Example script to generate a map using ggplot2 with similar look to Ocean Data View ###
	###########################################################################################

	# The example map I am trying to recreate is here:
	# https://twitter.com/JamesGrecian/status/1137018574504153089?s=20

	# Load in libraries
	require(marmap) # for bathymetry data
	require(tidyverse)
	require(sf)
	# devtools::install_github("jamesgrecian/mapr", force = T)
	require(mapr) # my wrapper for shapefile generation

	# Some helper functions for colour ramps from Claus Wilke
	# https://github.com/tidyverse/ggplot2/issues/2673
	discrete_gradient_pal <- function(colours, bins = 11) {
	ramp <- scales::colour_ramp(colours)

	function(x) {
	if (length(x) == 0) return(character())

	i <- floor(x * bins)
	i <- ifelse(i > bins-1, bins-1, i)
	ramp(i/(bins-1))
	}
	}

	scale_colour_discrete_gradient <- function(..., colours, bins = 11, na.value = "grey50", guide = "colourbar", aesthetics = c("colour", "fill"), colors) {
	colours <- if (missing(colours))
	colors
	else colours
	continuous_scale(
	aesthetics,
	"discrete_gradient",
	discrete_gradient_pal(colours, bins),
	na.value = na.value,
	guide = guide,
	...
	)
	}

	# Extract depth data from NOAA
	# This may take some time...
	batdat <- marmap::getNOAA.bathy(lon1 = -80,
	lon2 = 20,
	lat1 = -40,
	lat2 = -90,
	resolution = 2)
	bat = marmap::as.raster(batdat)
	raster::projection(bat) = sp::CRS("+proj=longlat +datum=WGS84")

	# Define a projection and transform bathymetry data
	prj <- '+proj=laea +lat_0=-60 +lon_0=-40 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
	r = raster::projectRaster(bat, crs = sp::CRS(prj), method = 'ngb', res = 50000)
	bat_df = as.data.frame(raster::rasterToPoints(r))
	names(bat_df) = c("x", "y", "z")

	# Example plot ggplot
	p1 <- ggplot() +
	theme_bw() +
	geom_point(aes(x = x, y = y, colour = z), data = bat_df) +
	geom_raster(aes(x = x, y = y, fill = z), data = bat_df) +
	scale_colour_discrete_gradient("Depth (m)",
	colours = rev(RColorBrewer::brewer.pal(11, "RdYlBu")),
	limits = c(-8000, 800),
	breaks = seq(-8000, 800, 800),
	guide = guide_colourbar(nbin = 500,
	raster = T,
	frame.colour = "black",
	ticks.colour = "black",
	frame.linewidth = 1,
	barwidth = 1,
	barheight = 25,
	direction = "vertical",
	title.position = "right",
	title.theme = element_text(angle = 90,
	hjust = 0.5))) +
	geom_sf(aes(), data = mapr::mapr(ellie, prj, buff = 1e7)) +
	coord_sf(xlim = c(-3200000, 4200000), ylim = c(-2800000, 2100000), crs = prj, expand = T) +
	xlab("") +ylab("")

	print(p1)

	# To make a wedge shaped map create a mask that matches map extent
	# Help again from Claus Wilke
	# https://htmlpreview.github.io/?https://github.com/clauswilke/practical_ggplot2/blob/master/goode.html


	# Define a clipping polygon based on bathymetry raster
	pol <- sf::st_polygon(list(rbind(c(-80, -90),
	c( 20, -90),
	c( 20, -40),
	c(-80, -40),
	c(-80, -90))))
	# Segmentise to ensure projection becomes a 'wedge' and not a trapezoid
	pol <- sf::st_segmentize(pol, 1)
	pol <- pol %>% sf::st_sfc()
	pol <- pol %>% st_set_crs(4326)
	pol <- pol %>% sf::st_transform(prj)

	# get the bounding box in transformed coordinates and expand by 10%
	xlim <- st_bbox(pol)[c("xmin", "xmax")]*1.1
	ylim <- st_bbox(pol)[c("ymin", "ymax")]*1.1

	# turn into enclosing rectangle
	encl_rect <-
	list(
	cbind(
	c(xlim[1], xlim[2], xlim[2], xlim[1], xlim[1]),
	c(ylim[1], ylim[1], ylim[2], ylim[2], ylim[1])
	)
	) %>%
	st_polygon() %>%
	st_sfc(crs = prj)

	# calculate the area outside the earth outline as the difference
	# between the enclosing rectangle and the earth outline
	cookie <- st_difference(encl_rect, pol)

	# Add this to the plot to mask out the areas not of interest
	# More effort required to move axis labels...

	p2 <- ggplot() +
	theme_minimal() +
	geom_point(aes(x = x, y = y, colour = z), data = bat_df) +
	geom_raster(aes(x = x, y = y, fill = z), data = bat_df) +
	scale_colour_discrete_gradient("Depth (m)",
	colours = rev(RColorBrewer::brewer.pal(11, "RdYlBu")),
	limits = c(-8000, 800),
	breaks = seq(-8000, 800, 800),
	guide = guide_colourbar(nbin = 500,
	raster = T,
	frame.colour = "black",
	ticks.colour = "black",
	frame.linewidth = 1,
	barwidth = 1,
	barheight = 25,
	direction = "vertical",
	title.position = "right",
	title.theme = element_text(angle = 90,
	hjust = 0.5))) +
	geom_sf(aes(), data = mapr::mapr(ellie, prj, buff = 1e7)) +
	geom_sf(fill = "white", color = "grey", data = cookie) +
	coord_sf(xlim = c(-3200000, 4200000), ylim = c(-2800000, 2100000), crs = prj, expand = T) +
	xlab("") +ylab("")

	gridExtra::grid.arrange(p1, p2, ncol = 2)
No results found