This vignette describes a number of issues that did not come up in the previous vignettes, and that may or may not be categorized as “frequently asked questions”. Readers are encouraged to provide entries for this vignette (as for the others).

EPSG stands for a maintained, well-understood registry of spatial reference systems, maintained by the International Association of Oil & Gass Producers (IOGP). “EPSG” stands for the authority, e.g. “EPSG:4326” stands for spatial reference system with ID 4326 as it is maintained by the EPSG authority. The website for the EPSG registry is found at the epsg.org domain.

`sf`

deal with secondary geometry columns?
`sf`

objects can have more than one geometry list-column,
but always only one geometry column is considered *active*, and
returned by `st_geometry`

. When there are multiple geometry
columns, the default `print`

methods reports which one is
active:

```
library(sf)
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
demo(nc, ask = FALSE, echo = FALSE)
nc$geom2 = st_centroid(st_geometry(nc))
print(nc, n = 2)
## Simple feature collection with 100 features and 14 fields
## Attribute-geometry relationship: 0 constant, 8 aggregate, 6 identity, 1 NA's
## Active geometry column: geom
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS: NAD27
## First 2 features:
## AREA PERIMETER CNTY_ CNTY_ID NAME FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114 1.442 1825 1825 Ashe 37009 37009 5 1091 1
## 2 0.061 1.231 1827 1827 Alleghany 37005 37005 3 487 0
## NWBIR74 BIR79 SID79 NWBIR79 geom
## 1 10 1364 0 19 MULTIPOLYGON (((-81.47276 3...
## 2 10 542 3 12 MULTIPOLYGON (((-81.23989 3...
## geom2
## 1 POINT (-81.49823 36.4314)
## 2 POINT (-81.12513 36.49111)
```

We can switch the active geometry by using
`st_geometry<-`

or `st_set_geometry`

, as in

`plot(st_geometry(nc))`

```
st_geometry(nc) <- "geom2"
plot(st_geometry(nc))
```

`st_simplify`

preserve topology?
`st_simplify`

is a topology-preserving function, but does
this on the level of individual feature geometries. That means, simply
said, that after applying it, a polygon will still be a polygon. However
when two features have a longer shared boundary, applying
`st_simplify`

to the object does not guarantee that in the
resulting object these two polygons still have the same boundary in
common, since the simplification is done independently, *for each
feature geometry*.

`sf`

objects?
They do! However, many developers like to write scripts that never
load packages but address all functions by the `sf::`

prefix,
as in

`i = sf::st_intersects(sf1, sf2)`

This works up to the moment that a `dplyr`

generic like
`select`

for an `sf`

object is needed: should one
call `dplyr::select`

(won’t know it should search in package
`sf`

) or `sf::select`

(which doesn’t exist)?
Neither works. One should in this case simply load `sf`

,
e.g. by

Most (but not all) of the geometry calculating routines used by
`sf`

come from the GEOS
library. This library considers coordinates in a two-dimensional, flat,
Euclidean space. For longitude latitude data, this is not the case. A
simple example is a polygon enclosing the North pole, which should
include the pole:

```
polygon = st_sfc(st_polygon(list(rbind(c(0,80), c(120,80), c(240,80), c(0,80)))),
crs = 4326)
pole = st_sfc(st_point(c(0,90)), crs = 4326)
st_intersects(polygon, pole)
## Sparse geometry binary predicate list of length 1, where the predicate
## was `intersects'
## 1: 1
```

which gives a wrong result (no intersection).

Similar to the above, centroids are computed assuming flat, 2D space:

```
st_centroid(polygon)[[1]]
## POINT (0 90)
```

where the centroid should have been the pole.