How to create contour lines from simple features of an sf object
I had a task to draw counter lines from a simple feature of an
sf object of polygons. The question is: Can we directly draw contour lines from polygons?
In what follows, I first convert the polygons to lines and you will see this will not do the trick. Then I provide an alternative way to do so.
The basic idea is that we first convert the polygons to points and then we interpolate the points to a raster. Finally, we draw contour lines from that raster. The workhorse is the
st_contour function in the
Prepare the data
Here I use the Agricultural Statistics Districts (ASD) for the entire U.S. as our spatial data.
An ASD is defined as a contiguous group of counties having relatively similar agricultural characteristics. The raw data (shapefiles) can be downloaded from here.
I have already converted the shapefile to a
library(tidyverse) library(sf) library(stars) library(gstat) library(tmap) library(automap) dataMap <- st_read(dsn = 'contour.gpkg')
## Reading layer `contour' from data source ## `C:\work\PersonalWebsite\mywebsite\content\notes_in_r\2021-08-29-how-to-create-contour-lines-from-simple-features-of-an-sf-object\contour.gpkg' ## using driver `GPKG' ## Simple feature collection with 303 features and 4 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -124.7332 ymin: 24.5447 xmax: -66.9499 ymax: 49.38435 ## Geodetic CRS: NAD83
# Here I use the ESRI:102004 projection dataMap <- st_transform(dataMap, "ESRI:102004")
I created a
contour column which contains arbitrary values.
## Simple feature collection with 6 features and 4 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: 708708.4 ymin: -940047.8 xmax: 1052223 ymax: -389495.9 ## Projected CRS: USA_Contiguous_Lambert_Conformal_Conic ## STATE STASD_N STASD_A contour geom ## 1 01 110 0110 141 MULTIPOLYGON (((733752.1 -4... ## 2 01 120 0120 190 MULTIPOLYGON (((945896.1 -3... ## 3 01 130 0130 76 MULTIPOLYGON (((745661.9 -5... ## 4 01 140 0140 70 MULTIPOLYGON (((760374.9 -6... ## 5 01 150 0150 112 MULTIPOLYGON (((765076.4 -9... ## 6 01 160 0160 167 MULTIPOLYGON (((1032953 -65...
tm_shape(dataMap) + tm_polygons(col = 'contour')
First we try converting polygons to lines
Let’s see what it looks like if we convert the polygons to lines. Note that the
tmap package provides a
tm_iso which, according to the documentation, should be used to draw contour lines. The function is essentially a combination of
tm_text. In the code below, we just use
dataLines <- st_cast(dataMap, 'MULTILINESTRING') tm_shape(dataLines) + tm_lines(col = 'contour')
We end up getting a map of ‘border’ with the color indicating the value.
An alternative way to do the task
It turns out that contour lines need to be drawn from raster objects to make the lines smooth. In cases of polygons, we can directly rasterize the polygons and draw contour lines. But you will see that this is not we want. Another way is to first convert the polygons to points and then interpolate the points to a raster object and we draw contour lines from there. The interpolation acts as a smoothing process.
After the raster object is created, we will use the
st_contour function from the
stars package to draw contour lines. We will also need to define the breaks where contour lines are made.
Let us first examine how it looks like if we draw the contour lines from rasterized polygons.
# rasterize the polygons dataRaster <- st_rasterize(dataMap['contour']) # Define the breaks for drawing contour lines breaks = seq(0, 200, 5) dataRasterLines <- st_contour(dataRaster, breaks = breaks, contour_lines = T) tm_shape(dataRasterLines) + tm_lines(col = 'contour')
This map looks no different from the “border” map above. Next let us try the second approach.
# We first convert the polygon object to points (i.e., centroids). dataPoints <- st_point_on_surface(dataMap)
## Warning in st_point_on_surface.sf(dataMap): st_point_on_surface assumes ## attributes are constant over geometries of x
# We prepare the empty raster that we will interpolate the data to. # Note that the resulted raster is regular. dataRasterBase <- st_bbox(dataMap) %>% st_as_stars(dx = 10000) %>% st_crop(dataMap) # Here we run a simple kriging interpolation v_mod_ok = autofitVariogram(contour ~ 1, as(dataPoints, "Spatial")) gg = gstat::gstat(formula = contour ~ 1, model = v_mod_ok$var_model, data = dataPoints) zz = predict(gg, dataRasterBase)
## [using ordinary kriging]
zz = zz["var1.pred", ] names(zz) = "contour" # Let us take a look at the raster object we just created plot(zz)
# Then we convert to contour lines using st_contour breaks1 <- seq(80, 120, 1) contourline <- st_contour(zz, contour_lines = T, na.rm = T, breaks = breaks1) tm_shape(dataMap) + tm_borders(col = "gray80") + tm_shape(contourline) + tmap::tm_lines(col = "contour", lwd = 3, palette = "viridis") + tm_layout(legend.title.size = 1.5, legend.text.size = 1, legend.position = c("left", "bottom"), legend.format = list(digits = 0))
This map looks much better. Note that the more detailed you define the breaks in
st_contour, the more lines you will get.
Please note that there may exist more efficiency ways to draw contour lines from polygons. Please do not hesitate to let me know if you find some.