sf tutorial

1. `sf`中的空间对象

1.1 `sf`与`simple features`

1.1.1 `sf`包简介

sf包集成了GDAL、PROJ、GEOS、S2的功能，提供了一个R语言的接口，它支持简单要素类simple features，并提供了丰富的函数进行数据转换、分析和可视化。与R语言中其他空间分析包(如terra、stars、qgisprocess)兼容性良好，并支持很多tidyverse中的数据操作函数。同时sf包也具备高效的性能和强大的空间分析功能，如缓冲区分析、叠加分析、创建渔网等。sf类似于Python中的geopandas但比geopandas的API设计要更简洁和规范(sf参考了很多PostGIS的设计),同时sf默认地理坐标系使用谷歌开发的S2空间索引算法运算，相较于GEOS在planr下的地理坐标系运算有较大的进步。sf包主要由Edzer Pebesma教授维护和开发(edzer教授同时也是sp、gstat、stars、s2、units等包的作者)，关于sf包的bug反馈请前往sf的github repo

1.1.2 `simple features` 简单要素类

simple features是一个由OGC开发和推广的地理空间矢量数据标准，它是一个层次化的数据模型。在该数据标准支持的最高达18种几何类型中，只有7种被大多数地理空间数据广泛使用(如下图所示)，而sf包完全支持这7种常用的几何类型,我们在后面将详细介绍这7种常用的几何类型。

1.2 常用的7种几何类型

sf中常用的7种几何类型分别是point 与 multipoint、linestring 与 multilinestring、polygon 与 multipolygon以及geometry collection，我们需要输入R中的数据对象从而使用sf提供的函数从头创建几何类型，具体来说我们使用st_前缀加上对应的小写几何类型组成的函数结合R中的数据对象从头创建sf中的几何类型(如使用st_point和数字向量创建点，st_linestring和矩阵创建线段)。 R中的数据对象和使用sf包可以创建的几何类型的对应关系如下:

vector: 对应一个单独的点point
matrix: 对应多点对象multipoint、单独的线段linestring
list: 矩阵对象的集合;对应multilinestring、polygon、multipolygon、 geometry collections

library(sf)

# XY point
st_point(c(5, 2))
# XYZ point
st_point(c(5, 2, 3))              
# XYM point
st_point(c(5, 2, 1), dim = "XYM") 
# XYZM point
st_point(c(5, 2, 3, 1)) 
# MULTIPOINT
multipoint_matrix = matrix(c(5, 2, 1, 3, 3, 4, 3, 2),
                           ncol = 2,byrow = T)
st_multipoint(multipoint_matrix)
plot(st_multipoint(multipoint_matrix),axes = T, col = 'red')

# LINESTRING
linestring_matrix = matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
                           ncol = 2,byrow = T)
st_linestring(linestring_matrix)
plot(st_linestring(linestring_matrix),axes = T, col = 'red')

# POLYGON 需要注意从头创建面时起点和终点点坐标要重合以保证面闭合
polygon_list = list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
                           ncol = 2,byrow = T))
st_polygon(polygon_list)
plot(st_polygon(polygon_list),axes = T,col = 'grey', border = 'red')

# POLYGON with a hole 内部有空洞的面
polygon_border = matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
                        ncol = 2,byrow = T)
polygon_hole = matrix(c(2, 4, 3, 4, 3, 3, 2, 3, 2, 4),
                      ncol = 2,byrow = T)
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list)
plot(st_combine(st_sfc(st_polygon(list(polygon_border)),
                       st_polygon(list(polygon_hole)))),
                axes = T, col = "transparent", border = "red")

# MULTILINESTRING
multilinestring_list = list(matrix(c(1, 5, 4, 4, 4, 1, 2, 2, 3, 2),
                                   ncol = 2,byrow = T), 
                            matrix(c(1, 2, 2, 4),
                                   ncol = 2,byrow = T))
st_multilinestring(multilinestring_list)
plot(st_multilinestring(multilinestring_list), col = "red", axes = T)

# MULTIPOLYGON
multipolygon_list = list(list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
                                     ncol = 2,byrow = T)),
                         list(matrix(c(0, 2, 1, 2, 1, 3, 0, 3, 0, 2),
                                     ncol = 2,byrow = T)))
st_multipolygon(multipolygon_list)
plot(st_multipolygon(multipolygon_list),axes = T, border = "red", col = "grey")

# GEOMETRYCOLLECTION
geometrycollection_list = list(st_multipoint(multipoint_matrix),
                               st_linestring(linestring_matrix))
st_geometrycollection(geometrycollection_list)
plot(st_geometrycollection(geometrycollection_list),axes = T, col = 'grey')

这里我们从头创建的几何类型都是sfg类，我们知道实际上空间矢量数据包括属性字段和几何要素，那么在sf包中空间矢量数据是如何组织的呢，请继续往下看：

1.3 `sf`中的空间数据组织

简单要素类由两部分组成: 几何对象和非空间的属性字段，下图展示了一个sf对象是如何被创建的。

在sf包中，创建sf类的几何对象来自sfc类，属性字段来自data.frame或tibble.sfg(simple feature geometries)即为我们在1.2部分提到的7种几何类型,然而一个sfg对象只包括一个简单要素类几何对象.实际上，一个空间矢量数据可能具备很多几何对象(比如西安市有11个区县级行政区划面).一个sfc对象就是一个或多个sfg对象的集合,sfc对象还能够包含关于所使用的坐标系统(坐标系统这里不再赘述,请参考这篇知乎文章)的信息.要将两个sfg对象合成一个具有sfc对象,我们可以使用st_sfc()函数.

point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc

## Geometry set for 2 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1 ymin: 2 xmax: 5 ymax: 3
## CRS:           NA

我们可以通过st_geometry_type()查看sfc对象的几何类型

st_geometry_type(points_sfc)

## [1] POINT POINT
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

不同类型的sfg对象也可以合并为一个sfc对象.

st_sfc(st_polygon(polygon_with_hole_list),
       st_multilinestring(multilinestring_list))

## Geometry set for 2 features 
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: 1 ymin: 1 xmax: 4 ymax: 5
## CRS:           NA

你可能已经注意到输出信息中CRS显示为NA,前面部分也提到sfc对象还可以存储坐标系统信息。默认值是 NA (Not Vacable) ，可以用st_crs()查看sfc对象的坐标系统:

st_crs(points_sfc)

## Coordinate Reference System: NA

我们可以在创建sfc对象时通过指定crs参数以定义坐标系统:

points_sfc = st_sfc(point1, point2, crs = "EPSG:4326")
points_sfc

## Geometry set for 2 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1 ymin: 2 xmax: 5 ymax: 3
## Geodetic CRS:  WGS 84

也可以在创建完sfc对象之后使用st_set_crs()为其设置坐标系统:

points_sfc = st_sfc(point1, point2)
points_sfc = st_set_crs(points_sfc,"EPSG:4326")
points_sfc

## Geometry set for 2 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1 ymin: 2 xmax: 5 ymax: 3
## Geodetic CRS:  WGS 84

sfc对象是sf对象中的几何列.sf对象额外包括了一些属性信息，下面以一个简单的sf对象创建为例演示:

xian_point = st_point(c(108.56, 34.15)) 
xian_geom = st_sfc(xian_point, crs = "EPSG:4326")
xian_attrib = tibble::tibble(                   
  name = "Xi'an",
  temperature = units::set_units(1,'°C'),
  date = as.Date("2024-02-05")
)
xian_sf = st_sf(xian_attrib, geometry = xian_geom) 
xian_sf

## Simple feature collection with 1 feature and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 108.56 ymin: 34.15 xmax: 108.56 ymax: 34.15
## Geodetic CRS:  WGS 84
## # A tibble: 1 × 4
##   name  temperature date             geometry
##   <chr>        [°C] <date>        <POINT [°]>
## 1 Xi'an           1 2024-02-05 (108.56 34.15)

1.4 `sf`包空间数据的读入和写出

常用的矢量数据格式有shapefile、geojson、geopackage、geodatabase等，根据矢量数据的存储格式,使用sf读取数据略有不同:针对shapefile格式的数据,读取时提供.shp后缀的路径即可;针对geojson格式的数据,读取时提供.geojson后缀的路径即可;而geopackage、geodatabase是空间数据库，不仅需要提供.gdb和.gpkg格式的后缀，还需提供layer参数以指定加载的图层，如果是单一图层的空间数据库则无需指定layer参数.sf提供st_read()和read_sf()两个函数以读取矢量数据，两者区别在于后者读取数据时属性列用tibble存储,和tidyverse更搭。下面以sf包自带的数据为例演示数据的读入:

system.file("gpkg/nc.gpkg", package="sf")

## [1] "C:/Users/dell/R/pkgR/sf/gpkg/nc.gpkg"

system.file("gpkg/nc.gpkg", package="sf") |> 
  read_sf() -> nc
nc

## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## # A tibble: 100 × 15
##     AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS  FIPSNO CRESS_ID BIR74 SID74 NWBIR74
##    <dbl>     <dbl> <dbl>   <dbl> <chr> <chr>  <dbl>    <int> <dbl> <dbl>   <dbl>
##  1 0.114      1.44  1825    1825 Ashe  37009  37009        5  1091     1      10
##  2 0.061      1.23  1827    1827 Alle… 37005  37005        3   487     0      10
##  3 0.143      1.63  1828    1828 Surry 37171  37171       86  3188     5     208
##  4 0.07       2.97  1831    1831 Curr… 37053  37053       27   508     1     123
##  5 0.153      2.21  1832    1832 Nort… 37131  37131       66  1421     9    1066
##  6 0.097      1.67  1833    1833 Hert… 37091  37091       46  1452     7     954
##  7 0.062      1.55  1834    1834 Camd… 37029  37029       15   286     0     115
##  8 0.091      1.28  1835    1835 Gates 37073  37073       37   420     0     254
##  9 0.118      1.42  1836    1836 Warr… 37185  37185       93   968     4     748
## 10 0.124      1.43  1837    1837 Stok… 37169  37169       85  1612     1     160
## # ℹ 90 more rows
## # ℹ 4 more variables: BIR79 <dbl>, SID79 <dbl>, NWBIR79 <dbl>,
## #   geom <MULTIPOLYGON [°]>

system.file("gpkg/nc.gpkg", package="sf") |> 
  st_read() -> nc

## Reading layer `nc.gpkg' from data source `C:\Users\dell\R\pkgR\sf\gpkg\nc.gpkg' using driver `GPKG'
## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27

nc

## Simple feature collection with 100 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## Geodetic CRS:  NAD27
## First 10 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
## 3  0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
## 4  0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
## 5  0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
## 6  0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
## 7  0.062     1.547  1834    1834      Camden 37029  37029       15   286     0
## 8  0.091     1.284  1835    1835       Gates 37073  37073       37   420     0
## 9  0.118     1.421  1836    1836      Warren 37185  37185       93   968     4
## 10 0.124     1.428  1837    1837      Stokes 37169  37169       85  1612     1
##    NWBIR74 BIR79 SID79 NWBIR79                           geom
## 1       10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
## 2       10   542     3      12 MULTIPOLYGON (((-81.23989 3...
## 3      208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
## 4      123   830     2     145 MULTIPOLYGON (((-76.00897 3...
## 5     1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
## 6      954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...
## 7      115   350     2     139 MULTIPOLYGON (((-76.00897 3...
## 8      254   594     2     371 MULTIPOLYGON (((-76.56251 3...
## 9      748  1190     2     844 MULTIPOLYGON (((-78.30876 3...
## 10     160  2038     5     176 MULTIPOLYGON (((-80.02567 3...

写出矢量格式的数据使用st_write()或write_sf()即可,首先给出要写入的sf/sfc/sfg对象，接着提供文件路径，如果是geopackage和geodatabase格式的文件还需指定layer参数，重复导出至相同的数据文件请指定overwrite参数为TRUE.(PS:属性列中带中文的导出为shp由于GDAL引擎编码的问题可能会乱码，请参考这篇博客解决)

write_sf(nc,'nc.gpkg',layer = 'nc',overwrite = TRUE)

1.5 关于`sf`中`s2`引擎的说明

sf包加载时显示的信息会提示你相应的GDAL、PROJ、GEOS版本，你可能也注意到了特殊的提示:sf_use_s2() is TRUE。此时sf默认打开s2运算引擎，当数据对象是地理坐标系时，sf将使用s2作为后端计算引擎而不是geos,如果你想使用geos进行地理坐标系下的数据运算，请使用这行命令关闭s2引擎:sf_use_s2(FALSE);关于数据是否是地理坐标系请使用st_is_longlat()查看；无论s2引擎是否关闭，投影坐标系下的运算都是geos引擎~

1.6 常用`sf`构造数据的一些小tips

WGS84经纬度坐标点转矢量点对象:

tibble::tibble(
  name = c('沈阳','长春','哈尔滨','北京','天津','呼和浩特','银川','太原','石家庄','济南','郑州','西安','武汉','南京','合肥','上海','长沙','南昌','杭州','福州','广州','台北','海口','南宁','重庆','昆明','贵阳','成都','兰州','西宁','拉萨','乌鲁木齐','香港','澳门'),
  lat = c(41.796768,43.886841,45.756966,39.904987,39.125595,40.84149,38.48644,37.857014,38.045475,36.675808,34.757977,34.263161,30.584354,32.041546,31.861191,31.231707,28.19409,28.676493,30.287458,26.075302,23.125177,25.030724,20.04422,22.82402,29.533155,25.040609,26.578342,30.659462,36.06138,36.61729,29.64415,43.82663,22.27534,22.19875),
  long = c(123.429092,125.324501,126.642464,116.405289,117.190186,111.75199,106.23248,112.549248,114.502464,117.000923,113.665413,108.948021,114.298569,118.76741,117.283043,121.472641,112.982277,115.892151,120.15358,119.306236,113.28064,121.520076,110.19989,108.320007,106.504959,102.71225,106.713478,104.065735,103.83417,101.77782,91.1145,87.61688,114.16546,113.54913)
) |> 
  st_as_sf(coords = c('long','lat'),
           crs = "EPSG:4326") -> provincep
provincep

## Simple feature collection with 34 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 87.61688 ymin: 20.04422 xmax: 126.6425 ymax: 45.75697
## Geodetic CRS:  WGS 84
## # A tibble: 34 × 2
##    name                geometry
##  * <chr>            <POINT [°]>
##  1 沈阳     (123.4291 41.79677)
##  2 长春     (125.3245 43.88684)
##  3 哈尔滨   (126.6425 45.75697)
##  4 北京     (116.4053 39.90499)
##  5 天津     (117.1902 39.12559)
##  6 呼和浩特  (111.752 40.84149)
##  7 银川     (106.2325 38.48644)
##  8 太原     (112.5492 37.85701)
##  9 石家庄   (114.5025 38.04548)
## 10 济南     (117.0009 36.67581)
## # ℹ 24 more rows

plot(provincep)

2. `sf`中的空间运算

2.1 `predicates` 谓词运算

2.1.1 `unary predicates` 一元谓词运算

一元谓词运算描述了空间几何对象的某种性质

st_is_simple 判断矢量对象是否自相交self-intersecting;自相交为FALSE否则为TRUE

library(sf)
library(tidyverse)

ls = st_linestring(rbind(c(0,0), c(1,1), c(1,0), c(0,1)))
plot(ls)

tibble(
  long = c(110,120),
  lat = c(10,20)
) |> 
  st_as_sf(coords = c('long','lat'),
           crs = 4326) |> 
  st_bbox() |> 
  st_as_sfc() -> p
plot(p)

st_is_simple(ls)

## [1] FALSE

st_is_simple(p)

## [1] TRUE

st_is_valid 判断矢量几何对象拓扑是否正确 错误的可以使用st_make_valid修复 st_is_empty 判断矢量几何对象是否为空 st_is_longlat 判断矢量几何对象是否为地理坐标系

2.1.2 `binary predicates` 二元谓词运算

二元谓词运算涉及到两个空间对象间的拓扑关系,我们在这里不深入讲解，有兴趣的可以搜索 DE-9IM进行了解

二元谓词运算具体如下所示:

predicate	meaning	inverse of
`contains`	None of the points of A are outside B	`within`
`contains_properly`	A contains B and B has no points in common with the boundary of A
`covers`	No points of B lie in the exterior of A	`covered_by`
`covered_by`	Inverse of `covers`
`crosses`	A and B have some but not all interior points in common
`disjoint`	A and B have no points in common	`intersects`
`equals`	A and B are topologically equal: node order or number of nodes may differ; identical to A contains B and A within B
`equals_exact`	A and B are geometrically equal, and have identical node order
`intersects`	A and B are not disjoint	`disjoint`
`is_within_distance`	A is closer to B than a given distance
`within`	None of the points of B are outside A	`contains`
`touches`	A and B have at least one boundary point in common, but no interior points
`overlaps`	A and B have some points in common; the dimension of these is identical to that of A and B
`relate`	Given a mask pattern, return whether A and B adhere to this pattern

其中st_relate可以基于DE-9IM字符串定义特殊的拓扑关系;所有上述的函数前两个参数都是x和y；这里的x和y可以是sfg、sfc、sf对象；如果只对函数提供一个操作对象，则默认对x及其自身做谓词运算。用粉红色表示x,用蓝色表示y,对应的二元谓词运算(二元拓扑关系运算)图解如下图所示:

2.2 `measures` 测量运算

2.2.1 `unary measures` 一元测量运算

一元测量运算返回一个描述几何属性的测度或量

measure	returns
`dimension`	0 for points, 1 for linear, 2 for polygons, possibly `NA` for empty geometries
`area`	the area of a geometry
`length`	the length of a linear geometry

st_dimension判断矢量几何对象形状维度,points为0,lines为1,surfaces为2,空的几何对象为NA(当NA_if_empty参数保持默认值TRUE时)即默认情况下空的几何对象返回NA

x = st_sfc(
  st_point(0:1),
  st_linestring(rbind(c(0,0),c(1,1))),
  st_polygon(list(rbind(c(0,0),c(1,0),c(0,1),c(0,0)))),
  st_multipoint(),
  st_linestring(),
  st_geometrycollection())
st_dimension(x)

## [1]  0  1  2 NA NA NA

st_area 返回矢量几何对象围成的面积 st_length 返回矢量几何对象长度

2.2.2 `binary measures`

st_distance 返回两个矢量几何对象间的距离

beijing = st_point(c(116.4053,39.90499)) |> 
  st_sfc(crs = "EPSG:4326")

xian = st_point(c(108.948021,34.263161)) |> 
  st_sfc(crs = "EPSG:4326")

st_distance(beijing,xian)

## Units: [m]
##          [,1]
## [1,] 911024.2

某些特殊的二元谓词运算也可以归属到二元测量运算，比如rook和queen邻接:

st_queen = function(x, y) st_relate(x, y, pattern = "F***T****")
st_rook = function(x, y) st_relate(x, y, pattern = "F***1****")

grid = st_make_grid(x, n = 3)
grid_sf = st_sf(grid)
grid_sf$queens = lengths(st_queen(grid, grid[5])) > 0
plot(grid, col = grid_sf$queens)

grid_sf$rooks = lengths(st_rook(grid, grid[5])) > 0
plot(grid, col = grid_sf$rooks)

2.3 `transformations` 转换运算

2.3.1 `unary transformations` 一元转换运算

一元转换运算在每个几何对象上作用并为每个几何对象返回一个新的几何对象(如质心、最小外接圆等)

transformer	returns a geometry …
`centroid`	of type `POINT` with the geometry’s centroid
`buffer`	that is larger (or smaller) than the input geometry, depending on the buffer size
`jitter`	that was moved in space a certain amount, using a bivariate uniform distribution
`wrap_dateline`	cut into pieces that no longer cover or cross the dateline
`boundary`	with the boundary of the input geometry
`convex_hull`	that forms the convex hull of the input geometry (@fig-vor)
`line_merge`	after merging connecting `LINESTRING` elements of a `MULTILINESTRING` into longer `LINESTRING`s.
`make_valid`	that is valid
`node`	with added nodes to linear geometries at intersections without a node; only works on individual linear geometries
`point_on_surface`	with a (arbitrary) point on a surface
`polygonize`	of type polygon, created from lines that form a closed ring
`segmentize`	a (linear) geometry with nodes at a given density or minimal distance
`simplify`	simplified by removing vertices/nodes (lines or polygons)
`split`	that has been split with a splitting linestring
`transform`	transformed or convert to a new coordinate reference system (@sec-cs)
`triangulate`	with Delauney triangulated polygon(s) (@fig-vor)
`voronoi`	with the Voronoi tessellation of an input geometry (@fig-vor)
`zm`	with removed or added `Z` and/or `M` coordinates
`collection_extract`	with sub-geometries from a `GEOMETRYCOLLECTION` of a particular type
`cast`	that is converted to another type
`+`	that is shifted over a given vector
`*`	that is multiplied by a scalar or matrix

上述提及的一元转换运算除了最后的+和*运算符涉及到矢量几何对象的仿射变换，其操作对象是sfg或sfc对象，其他的一元转换运算函数sfg、sfc、sf对象均适用。一元转换运算使用较多的是提取质心st_centroid()、缓冲区分析st_buffer()、投影变换st_transform():当提取的质心不在相应几何对象内部时而我们需要生成一个相应几何对象内部点时可以使用st_point_on_surface();在投影时可以通过st_can_transform()确定两个坐标系间是否可以进行投影变换.

2.3.2 `binary transformations` 二元转换运算

二元转换运算对一对几何矢量对象执行几何运算

function	returns	infix operator
`intersection`	the overlapping geometries for pair of geometries	`&`
`union`	the combination of the geometries; removes internal boundaries and duplicate points, nodes or line pieces	`\|`
`difference`	the geometries of the first after removing the overlap with the second geometry	`/`
`sym_difference`	the combinations of the geometries after removing where they intersect; the negation (opposite) of `intersection`	`%/%`

需要注意的是使用&,|,/,%/%四个运算符时需要提供sfc或sfg对象，相应的操作函数(st_前缀函数)则是sf,sfc,sfg对象均可

图解如下(其中左边部分圆为x,右边部分圆为y,灰色区域代表运算结果):

2.3.3 `n-ary transformations` 多元转换运算

多元转换运算是对二元转换运算的拓展，最常用的是多个矢量几何对象的联合操作st_union, 参考二元转换运算理解即可，这里不过多赘述.

3. `sf`中几何类型转化

3.1 问题引入:

在实际空间数据处理过程中，对空间对象进行几何类型转换是一个经常会遇到的需求.比如我们可能需要根据获取的GPS点位数据生成运动轨迹;根据线生成面等。在sf包中这一操作主要通过st_cast()实现,本文就详细讨论sf中几何类型转化问题。

3.2 从`sfg`对象开始的几何类型转化

我们首先创建一个具有3个点的MULTIPOINT对象:

multipoint = st_multipoint(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))
plot(multipoint, col = 'red', cex = .5, axes = T)

把这个多点对象转换为LINESTRING 单一的线段:

linestring = st_cast(multipoint, "LINESTRING")
plot(linestring, col = 'red', axes = T)

把这个多点对象转换为POLYGON 一个面:

polyg = st_cast(multipoint, "POLYGON")
plot(polyg, col = "grey", border = "red", axes = TRUE)

我们同样也可以将POLYGON和LINESTRING转化为MULTIPOINT，他是前面转换操作的逆操作

multipoint_2 = st_cast(linestring, "MULTIPOINT")
multipoint_3 = st_cast(polyg, "MULTIPOINT")
all.equal(multipoint, multipoint_2)

## [1] TRUE

all.equal(multipoint, multipoint_3)

## [1] TRUE

我们同样也可以将单个的点线面sfg对象转换为对应的多点、多线、多面对象

point_single = st_point(c(1,1))
print(point_single)
point_multi = st_cast(point_single,'MULTIPOINT')
print(point_multi)
line_single = st_linestring(matrix(c(1, 3, 5, 1, 3, 1), ncol = 2))
print(line_single)
line_multi = st_cast(line_single,"MULTILINESTRING")
print(line_multi)
polyg_single = st_polygon(list(matrix(c(1, 1, 2, 1, 2, 2, 1, 2, 1, 1),
                                      ncol = 2,byrow = T)))
print(polyg_single)
polyg_multi = st_cast(polyg_single,"MULTIPOLYGON")
print(polyg_multi)

值得一提的是,当multi-类型的sfg对象转换为single的sfg对象时,仅会保留multi- 类型的sfg对象中的第一个

polyg = st_multipolygon(list(list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
                                         ncol = 2,byrow = T)),
                             list(matrix(c(0, 2, 1, 2, 1, 3, 0, 3, 0, 2),
                                         ncol = 2,byrow = T))))

print(polyg)
print(st_cast(polyg,"POLYGON"))

但当multi-类型的sfc和sf对象转换为single的sfc和sf对象时，则会拆分multi 类型至多个single类型:

polyg = st_multipolygon(list(list(matrix(c(1, 5, 2, 2, 4, 1, 4, 4, 1, 5),
                                         ncol = 2,byrow = T)),
                             list(matrix(c(0, 2, 1, 2, 1, 3, 0, 3, 0, 2),
                                         ncol = 2,byrow = T)))) |> 
  st_sfc()

print(polyg)

## Geometry set for 1 feature 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 1 xmax: 4 ymax: 5
## CRS:           NA

print(st_cast(polyg,"POLYGON"))

## Geometry set for 2 features 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 0 ymin: 1 xmax: 4 ymax: 5
## CRS:           NA

大部分情况下sfg对象的几何类型转化与sfc、sf对象工作一致，一个重要的区别就是前面提到的multi-类型转换为single类型:sfg对象仅会保留multi-类型中的第一个而sfc、sf对象会拆分multi类型至多个single类型.

3.3 `sfc`与`sf`对象的几何类型转化

sfc对象与sf对象的几何类型转化原理一致，只是sf对象相较于sfc对象多了属性列，涉及到了属性几何关系(AGR),我们下面以sf对象为例详细解释其几何类型转化,假设我们有下面这7个sf对象:

POI : 1个点
MPOI : 具有4个点的multipoint
LIN : 5个点的linestring
MLIN : 1个5个点的linestring和1个4个点的linestring组成的multilinestring
POL : 5个点组成的polygon
MPOL : 2对5个点的polygon组成的multipolygon
GC : 1个4个点的multipoint和1个5个点的linestring组成的geometry collection

这7个sf对象可能的几何类型变化如下表所示(NA代表操作不能实现,数字代表相应转化后产生的目标对象类型的个数)

3.4 一个具体的例子

下面我们以sf自带的nc.gpkg数据为例演示st_cast()的使用:

nc = system.file('gpkg/nc.gpkg',package = 'sf') |> 
  read_sf()
p = nc[1:3,]
p = st_geometry(p) |> st_as_sf()
plot(st_geometry(p))

st_cast(p,'MULTILINESTRING')

## Simple feature collection with 3 features and 0 fields
## Geometry type: MULTILINESTRING
## Dimension:     XY
## Bounding box:  xmin: -81.74107 ymin: 36.23388 xmax: -80.43531 ymax: 36.58965
## Geodetic CRS:  NAD27
##                                x
## 1 MULTILINESTRING ((-81.47276...
## 2 MULTILINESTRING ((-81.23989...
## 3 MULTILINESTRING ((-80.45634...

st_cast(p,'MULTIPOINT')

## Simple feature collection with 3 features and 0 fields
## Geometry type: MULTIPOINT
## Dimension:     XY
## Bounding box:  xmin: -81.74107 ymin: 36.23388 xmax: -80.43531 ymax: 36.58965
## Geodetic CRS:  NAD27
##                                x
## 1 MULTIPOINT ((-81.47276 36.2...
## 2 MULTIPOINT ((-81.23989 36.3...
## 3 MULTIPOINT ((-80.45634 36.2...

st_cast(p,'MULTIPOINT') |> 
  st_cast('POLYGON')

## Simple feature collection with 3 features and 0 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.74107 ymin: 36.23388 xmax: -80.43531 ymax: 36.58965
## Geodetic CRS:  NAD27
##                                x
## 1 POLYGON ((-81.47276 36.2343...
## 2 POLYGON ((-81.23989 36.3653...
## 3 POLYGON ((-80.45634 36.2425...

sf tutorial

Wenbo Lv

2024-10-13

1. `sf`中的空间对象

1.1 `sf`与`simple features`

1.1.1 `sf`包简介

1.1.2 `simple features` 简单要素类

1.2 常用的7种几何类型

1.3 `sf`中的空间数据组织

1.4 `sf`包空间数据的读入和写出

1.5 关于`sf`中`s2`引擎的说明

1.6 常用`sf`构造数据的一些小tips

WGS84经纬度坐标点转矢量点对象:

2. `sf`中的空间运算

2.1 `predicates` 谓词运算

2.1.1 `unary predicates` 一元谓词运算

2.1.2 `binary predicates` 二元谓词运算

2.2 `measures` 测量运算

2.2.1 `unary measures` 一元测量运算

2.2.2 `binary measures`

2.3 `transformations` 转换运算

2.3.1 `unary transformations` 一元转换运算

2.3.2 `binary transformations` 二元转换运算

2.3.3 `n-ary transformations` 多元转换运算

3. `sf`中几何类型转化

3.1 问题引入:

3.2 从`sfg`对象开始的几何类型转化

3.3 `sfc`与`sf`对象的几何类型转化

3.4 一个具体的例子

1. sf中的空间对象

1.1 sf与simple features

1.1.1 sf包简介

1.1.2 simple features 简单要素类

1.2 常用的7种几何类型

1.3 sf中的空间数据组织

1.4 sf包空间数据的读入和写出

1.5 关于sf中s2引擎的说明

1.6 常用sf构造数据的一些小tips

WGS84经纬度坐标点转矢量点对象:

2. sf中的空间运算

2.1 predicates 谓词运算

2.1.1 unary predicates 一元谓词运算

2.1.2 binary predicates 二元谓词运算

2.2 measures 测量运算

2.2.1 unary measures 一元测量运算

2.2.2 binary measures

2.3 transformations 转换运算

2.3.1 unary transformations 一元转换运算

2.3.2 binary transformations 二元转换运算

2.3.3 n-ary transformations 多元转换运算

3. sf中几何类型转化

3.1 问题引入:

3.2 从sfg对象开始的几何类型转化

3.3 sfc与sf对象的几何类型转化

3.4 一个具体的例子

1. `sf`中的空间对象

1.1 `sf`与`simple features`

1.1.1 `sf`包简介

1.1.2 `simple features` 简单要素类

1.3 `sf`中的空间数据组织

1.4 `sf`包空间数据的读入和写出

1.5 关于`sf`中`s2`引擎的说明

1.6 常用`sf`构造数据的一些小tips

2. `sf`中的空间运算

2.1 `predicates` 谓词运算

2.1.1 `unary predicates` 一元谓词运算

2.1.2 `binary predicates` 二元谓词运算

2.2 `measures` 测量运算

2.2.1 `unary measures` 一元测量运算

2.2.2 `binary measures`

2.3 `transformations` 转换运算

2.3.1 `unary transformations` 一元转换运算

2.3.2 `binary transformations` 二元转换运算

2.3.3 `n-ary transformations` 多元转换运算

3. `sf`中几何类型转化

3.2 从`sfg`对象开始的几何类型转化

3.3 `sfc`与`sf`对象的几何类型转化