ggplot label contours with sf object

user_2995175 0
Viewed 341

I have a contour dataset imported using st_read() from a shapefile. I have plotted this using ggplot() and geom_sf(). It renders nicely. Now I want to label the contours. Using geom_st_label() does not produce great output. Contours are dense in places and labels overlap.

I had a look at the metR package. This has a geom_contour_label() function that controls contour placement nicely in ggplot(). However, the geom_contour() and associated functions do not recognise the geometry contained in sf_objects. I get this error: stat_contour requires the following missing aesthetics: x, y and z.

How can I get geom_contour_label() to work with an sf object? This is what geom_contour_label() can produce:

enter image description here

My contour data is available from https://cloudstor.aarnet.edu.au/sender/?s=download&token=241de91b-2015-4a19-a18f-c2125a12f2a7.

library(sf)
library(tidyverse)

isobath <- read_sf("1misobath.shp")

ggplot() +
  geom_sf(data = isobath, color = "blue", lwd = 0.25) +
  geom_sf_label(data = isobath, aes(label = DEPTH), size = 2) +
  coord_sf(xlim = c(18.42, 18.5), ylim = c(-34.20, -34.16), expand = T) +
  theme_bw()
sf ggplot2 r
1 Answers

After examining the geom_contour() function in depth. Here is my understanding of the problem.

The function geom_contour() does not expect an object of type sf but a dataframe or a data.table which contains rows of points (not linestrings as in your case) with three columns: longitude (i.e. X), latitude (i.e. Y) and a Z column (i.e., in your case, this would be the variable DEPTH). The geom_contour() function will then compute the isolines and plot them on the map.

From your sf object, I was able to create a dataframe of points with X, Y and Z coordinates as the function expects, but then I was faced with an unsolvable problem. The grid of your data is rectilinear, but the function expects a grid of regular type. So, to convert the grid from rectilinear to regular, I tried to interpolate the points with the interp() function of the akima package but it was impossible because of the duplicated values in Z (i.e. variable DEPTH).

So I changed my mind! And finally, I think I found a very simple solution to your mapping problem by using the tm_iso() function of the very good tmap package.

Please find the reprex below showing the code and the resulting map result.

Reprex

library(sf) 
library(s2)
library(tmap)

# 1. Your data
isobath <- st_read("D:/test/Projet_essai4/Isobath/1misobath.shp")
#> Reading layer `1misobath' from data source 
#>   `D:\test\Projet_essai4\Isobath\1misobath.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 162 features and 9 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 18.43944 ymin: -34.19143 xmax: 18.49996 ymax: -34.16805
#> Geodetic CRS:  Cape


# 2. Visualization of your data with tm_iso() of the tmap package
tm_shape(isobath) +
  tm_grid(col ="lightgray", line = TRUE) +
  tm_iso(col = "blue", text = "DEPTH",  fontface = "bold") +
  tm_layout(frame = TRUE, inner.margins = 0)

Created on 2021-10-31 by the reprex package (v0.3.0)

  • Another version, perhaps closer to what you are looking for:
tm_shape(isobath) +
  tm_grid(col ="white", line = TRUE) +
  tm_iso(col = "blue", text = "DEPTH",  fontface = "bold", bg.color = "lightgray", shadow = TRUE, size = 0.6) +
  tm_layout(bg.color = "lightgray", frame = TRUE, inner.margins = 0)

Created on 2021-10-31 by the reprex package (v0.3.0)

user_136610610
Related