Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Plot outline around raster cells

Tags:

r

ggplot2

Using ggplot2 I would like to plot a polygon (convex hull?) that honours the cell boundaries of a raster plot (created with geom_raster). I can use convex hull methods but they would pass through the x-y coordinates rather than around the cells of the raster plot (sorry my terminology is affected by my ignorance). As an example, given the following x-y coordinate data with categorical attribute d:

df <- data.frame(x = c(295, 300, 305, 310, 295, 300, 305, 310),
                 y = c(310, 310, 310, 310, 315, 315, 315, 315),
                 d = c(2, 2, 2, 1, 2, 1, 1, 1))

I can use ggplot2's geom_raster to plot the attribute as a fill:

ggplot(df, aes(x, y)) +
  geom_raster(aes(fill = d)) +
  coord_fixed()

To give:

enter image description here

But what I want is an outline hull of the categories defined by d. Something like this:

enter image description here

Can I do this with base R and ggplot2 or is there a way using the raster package or some other grid/raster/GIS package? I prefer to keep it as simple as possible.

like image 646
Alex Trueman Avatar asked Mar 13 '23 15:03

Alex Trueman


2 Answers

You can certainly get those polygons using R's spatial facilities. Here's one way:

library(raster)

## Convert your data.frame to a raster object
r <- rasterFromXYZ(df)

## Extract polygons
pp <- rasterToPolygons(r, dissolve=TRUE)

## Convert SpatialPolygons to a format usable by ggplot2
outline <- fortify(pp)

## Put it all together:
ggplot(df, aes(x, y)) +
    geom_raster(aes(fill = d)) +
    coord_fixed() +
    geom_path(aes(x = long, y = lat, group = group), data = outline, 
              size=1.5, col="gold")

enter image description here

like image 111
Josh O'Brien Avatar answered Mar 17 '23 11:03

Josh O'Brien


Unfortunately, fortify() from {ggplot2} is no longer recommended by the package authors. The suggested replacement, tidy() from the {broom} package is also no longer recommend by those package authors!

Fortunately, we can recreate Josh O'Brien's answer using the {sf} package. We still rely on some functions from the {raster} package, but these might be replaceable with the {stars} package once a presumed bug is fixed (see https://github.com/r-spatial/sf/issues/1389). Using {stars} might be advantageous if your raster is very large, because converting the raster (or dataframe) to contiguous polygons can be very slow (see https://gis.stackexchange.com/a/313550/114075).

Here's the code (borrowing from Josh O'Brien's setup):

library(raster)
library(rgeos)
library(ggplot2)
library(sf)

df <- data.frame(x = c(295, 300, 305, 310, 295, 300, 305, 310),
                 y = c(310, 310, 310, 310, 315, 315, 315, 315),
                 d = c(2, 2, 2, 1, 2, 1, 1, 1))

## Convert your data.frame to a raster object
r <- raster::rasterFromXYZ(df)

## Extract polygons
pp <- raster::rasterToPolygons(r, dissolve = TRUE)

## Convert SpatialPolygons to a format usable by ggplot2
outline <- sf::st_as_sf(pp)

## Put it all together:
## The polygon "outline" is filled by default so will cover up the
## raster values. Use fill = NA for the geom_sf() to just use the
## gold border color. 
ggplot(df) +
  geom_raster(aes(x = x, y = y, fill = d)) +
  geom_sf(data = outline, size = 1.5, col = "gold", fill = NA)

Alternatively, cast the sf object to be a LINESTRING, then there is no need to use the fill = NA argument in the call to geom_sf():

library(raster)
library(rgeos)
library(ggplot2)
library(sf)

df <- data.frame(x = c(295, 300, 305, 310, 295, 300, 305, 310),
                 y = c(310, 310, 310, 310, 315, 315, 315, 315),
                 d = c(2, 2, 2, 1, 2, 1, 1, 1))

## Convert your data.frame to a raster object
r <- raster::rasterFromXYZ(df)

## Extract polygons
pp <- raster::rasterToPolygons(r, dissolve = TRUE)

## Or you can cast to a linestring instead of using a polygon
## and then having to use the fill = NA argument in the 
## geom_sf() call
outline <- sf::st_as_sf(pp) %>% st_cast("LINESTRING")

ggplot(df) +
  geom_raster(aes(x = x, y = y, fill = d)) +
  geom_sf(data = outline, size = 1.5, col = "gold")

enter image description here

like image 31
mikoontz Avatar answered Mar 17 '23 12:03

mikoontz