inst/doc/sf7.R

In sf: Simple Features for R

## ----echo=FALSE, include=FALSE------------------------------------------------
knitr::opts_chunk$set(fig.height = 4.5)
knitr::opts_chunk$set(fig.width = 6)
knitr::opts_chunk$set(collapse = TRUE)

## ---- eval=FALSE--------------------------------------------------------------
#  install.packages("s2")

## -----------------------------------------------------------------------------
library(sf)
sf_use_s2()

## -----------------------------------------------------------------------------
library(s2)

## -----------------------------------------------------------------------------
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf")) # wrong ring directions
s2_area(st_as_s2(nc, oriented = FALSE)[1:3]) # corrects ring direction, correct area:
s2_area(st_as_s2(nc, oriented = TRUE)[1:3]) # wrong direction: Earth's surface minus area
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"), check_ring_dir = TRUE)
s2_area(st_as_s2(nc, oriented = TRUE)[1:3]) # no second correction needed here:

## -----------------------------------------------------------------------------
g = as_s2_geography(TRUE)
g

## -----------------------------------------------------------------------------
co = s2_data_countries()
oc = s2_difference(g, s2_union_agg(co)) # oceans
b = s2_buffer_cells(as_s2_geography("POINT(-30 52)"), 9800000) # visible half
i = s2_intersection(b, oc) # visible ocean
plot(st_transform(st_as_sfc(i), "+proj=ortho +lat_0=52 +lon_0=-30"), col = 'blue')

## -----------------------------------------------------------------------------
s2_area(oc) / s2_area(g)

## -----------------------------------------------------------------------------
a = as_s2_geography("POINT(0 0)")
b = as_s2_geography("POLYGON((0 0,1 0,1 1,0 1,0 0))")
s2_intersects(a, b, s2_options(model = "open")) 
s2_intersects(a, b, s2_options(model = "closed"))
s2_intersects(a, b, s2_options(model = "semi-open")) # a toss
s2_intersects(a, b) # default: semi-open

## -----------------------------------------------------------------------------
fiji = s2_data_countries("Fiji")
aa = s2_data_countries("Antarctica")
s2_bounds_cap(fiji)
s2_bounds_rect(c(fiji,aa))

## -----------------------------------------------------------------------------
sf_use_s2()

## -----------------------------------------------------------------------------
sf_use_s2(FALSE)

## -----------------------------------------------------------------------------
sf_use_s2(TRUE)

## ----eval=require("lwgeom", quietly = TRUE)-----------------------------------
library(sf)
library(units)
nc = read_sf(system.file("gpkg/nc.gpkg", package="sf"))
sf_use_s2(TRUE)
a1 = st_area(nc)
sf_use_s2(FALSE)
a2 = st_area(nc)
plot(a1, a2)
abline(0, 1)
summary((a1 - a2)/a1)

## -----------------------------------------------------------------------------
nc_ls = st_cast(nc, "MULTILINESTRING")
sf_use_s2(TRUE)
l1 = st_length(nc_ls)
sf_use_s2(FALSE)
l2 = st_length(nc_ls)
plot(l1 , l2)
abline(0, 1)
summary((l1-l2)/l1)

## -----------------------------------------------------------------------------
sf_use_s2(TRUE)
d1 = st_distance(nc, nc[1:10,])
sf_use_s2(FALSE)
d2 = st_distance(nc, nc[1:10,])
plot(as.vector(d1), as.vector(d2))
abline(0, 1)
summary(as.vector(d1)-as.vector(d2))

## -----------------------------------------------------------------------------
sf_use_s2(TRUE)
st_intersects(nc[1:3,], nc[1:3,]) # self-intersections + neighbours
sf_use_s2(TRUE)
st_intersects(nc[1:3,], nc[1:3,], model = "semi-open") # only self-intersections

## ---- fig.show='hold', out.width="50%"----------------------------------------
uk = s2_data_countries("United Kingdom")
class(uk)
uk_sfc = st_as_sfc(uk) 
uk_buffer = s2_buffer_cells(uk, distance = 20000)
uk_buffer2 = s2_buffer_cells(uk, distance = 20000, max_cells = 10000)
uk_buffer3 = s2_buffer_cells(uk, distance = 20000, max_cells = 100)
class(uk_buffer)
plot(uk_sfc)
plot(st_as_sfc(uk_buffer))
plot(st_as_sfc(uk_buffer2))
plot(st_as_sfc(uk_buffer3))
uk_sf = st_as_sf(uk) 

## -----------------------------------------------------------------------------
# the sf way
system.time({
  uk_projected = st_transform(uk_sfc, 27700)
  uk_buffer_sf = st_buffer(uk_projected, dist = 20000)
})
# sf way with few than the 30 segments in the buffer
system.time({
  uk_projected = st_transform(uk_sfc, 27700)
  uk_buffer_sf2 = st_buffer(uk_projected, dist = 20000, nQuadSegs = 4)
})
# s2 with default cell size
system.time({
  uk_buffer = s2_buffer_cells(uk, distance = 20000)
})
# s2 with 10000 cells
system.time({
  uk_buffer2 = s2_buffer_cells(uk, distance = 20000, max_cells = 10000)
})
# s2 with 100 cells
system.time({
  uk_buffer2 = s2_buffer_cells(uk, distance = 20000, max_cells = 100)
})