Basic

library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")

map <- st_read(pathshp, quiet = TRUE)
class(map)

print(map)

st_crs(map)

head(st_drop_geometry(map))

st_geometry(map)

summary(map)

st_geometry_type(map)

plot(map[1]) # plot first attribute

st_bbox(map)

dim(map) # 数据框的行数（代表几何体的数量）和列数（包括属性列和几何列）。

str(map)

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- terra::rast(pathraster)
r
plot(r)

dim(r)

ext(r)

res(r)

crs(r)

minmax(r)

global(r, fun = "mean", na.rm = TRUE)

st_crs("EPSG:4326")$Name

st_crs("EPSG:4326")$proj4string

st_crs("EPSG:4326")$epsg

st_crs(map)$epsg

st_crs(map)$proj4string

st_crs(map)$Name

map2 <- st_transform(map, crs = "EPSG:4326")
st_crs(map2)

st_crs(map2)$Name

crs(r)

r2 <- terra::project(r, "EPSG:2169")
crs(r2)

# r2 is existing raster
# r is raster we project
r3 <- terra::project(r, r2)

library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(pathshp, quiet = TRUE)
class(nc)

dim(nc)

print(nc)

plot(nc)

colnames(nc)

nc$AREA

class(nc$AREA)

nc$AREA[1]

class(nc$AREA[1])

head(map)

geom_column <- st_geometry(map) # extract the sfc objects from the sf object
print(geom_column)

# extract the first geometry
single_geometry <- geom_column[[1]] 
# sfg (simple feature geometry): a single geometry
print(single_geometry)

class(geom_column)
class(single_geometry)

# Single point (point as a vector)
p1_sfg <- st_point(c(2, 2))
p2_sfg <- st_point(c(2.5, 3))

# Set of points (points as a matrix)
p <- rbind(c(6, 2), c(6.1, 2.6), c(6.8, 2.5),
           c(6.2, 1.5), c(6.8, 1.8))
mp_sfg <- st_multipoint(p)

# Polygon. Sequence of points that form a closed,
# non-self intersecting ring.
# The first ring denotes the exterior ring,
# zero or more subsequent rings denote holes in the exterior ring
p1 <- rbind(c(10, 0), c(11, 0), c(13, 2),
            c(12, 4), c(11, 4), c(10, 0))
p2 <- rbind(c(11, 1), c(11, 2), c(12, 2), c(11, 1))
pol_sfg <- st_polygon(list(p1, p2))

# Create sf object
p_sfc <- st_sfc(p1_sfg, p2_sfg, mp_sfg, pol_sfg)
df <- data.frame(TEST = c("A", "B", "C", "D"))
p_sf <- st_sf(df, geometry = p_sfc)

class(df)

class(p_sf)

print(p_sf)

plot(p_sf)

library(ggplot2)
ggplot(p_sf) + geom_sf(aes(col = TEST), size = 3) + theme_bw()

library(sf)
library(ggplot2)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

map$FIPS

# Delete polygon
map <- map[-which(map$FIPS %in% c("37125", "37051", "37077")), ]
# 是查找 map 数据框中 FIPS 列的值是 "37125" 或 "37051" 的行索引
# -which(...) 的作用是排除掉这些行
dim(map)
ggplot(map) + geom_sf(aes(fill = SID79))
# 用 SID79 列的值来填充颜色
# geom_sf() 是用于绘制 sf 对象的几何体（如多边形、线、点）的函数

# Combine geometries
new_map <- st_union(map, by_feature = FALSE)
ggplot(new_map %>% st_sf()) + geom_sf()
dim(new_map)
# st_union() 函数将多个几何体合并为一个
# by_feature = FALSE 表示将所有的几何体合并为单个几何对象，而不是保持原有的多个单独几何对象
# by_feature = TRUE 表示将所有的几何体合并为单个几何对象，并保留原有的多个单独几何对象
new_map1 <- st_union(map, by_feature = TRUE)
ggplot(new_map1 %>% st_sf()) + geom_sf()
dim(new_map1)
# Simplify
ggplot(st_simplify(new_map, dTolerance = 10000)) + geom_sf()
ggplot(st_simplify(new_map1, dTolerance = 10000)) + geom_sf()
dim(new_map1)

library(sf)
library(mapview)

d <- data.frame(
    place = c("London", "Paris", "Madrid", "Rome"),
    long = c(-0.118092, 2.349014, -3.703339, 12.496366),
    lat = c(51.509865, 48.864716, 40.416729, 41.902782),
    value = c(200, 300, 400, 500)
)
class(d)

dsf <- st_as_sf(d, coords = c("long", "lat"))
st_crs(dsf) <- 4326
class(dsf)

print(dsf)

plot(dsf[0])
plot(dsf[1])

library(sf)
library(ggplot2)

# Map with divisions (sf object)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

# Points over map (simple feature geometry list-column sfc)
points <- st_sample(map, size = 100) # 生成 100 个随机点

# Map of points within polygons
ggplot() + geom_sf(data = map) + geom_sf(data = points)

# Intersection (first argument map, then points)
inter <- st_intersects(map, points)

# Add point count to each polygon
map$count <- lengths(inter)

# Map of number of points within polygons
ggplot(map) + geom_sf(aes(fill = count))

library(sf)
library(ggplot2)

# Map with divisions (sf object)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

# Points over map (sf object)
points <- st_sample(map, size = 3) %>% st_as_sf()

# Intersection (first argument points, then map)
inter <- st_intersects(points, map)
print(inter)

# Adding column areaname with the name of
# the areas containing the points
points$areaname <- map[unlist(inter), "NAME",
                       drop = TRUE] # drop geometry
print(points)
print(points$areaname)

# Map
ggplot(map) + geom_sf() + geom_sf(data = points) + 
 geom_sf_label(data = map[unlist(inter), ],
               aes(label = NAME), nudge_y = 0.2)

library(sf)
# 创建两个示例数据集：一个多边形和一个点
polys <- st_as_sf(data.frame(id = 1:2, wkt = c("POLYGON((0 0, 0 2, 2 2, 2 0, 0 0))", 
                                               "POLYGON((3 0, 3 2, 5 2, 5 0, 3 0))")),
                  wkt = "wkt", crs = 4326)
points <- st_as_sf(data.frame(id = 3:4, wkt = c("POINT(1 1)", "POINT(4 1)")),
                   wkt = "wkt", crs = 4326)

# 进行空间连接，将点分配到所在的多边形
joined <- st_join(points, polys)
print(joined)

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

library(sf)
# 创建一个表示某个地点的点
point <- st_point(c(1, 1)) %>% st_sfc(crs = 4326) %>% st_sf()
# %>% 是 管道操作符，用于简化代码的写法
# 来自 magrittr 包（dplyr 包也包括这个操作符）
# 它允许你将一个表达式的结果传递给下一个函数，避免嵌套多个函数调用，提高代码的可读性。

# 创建一个围绕这个点的缓冲区，半径为1单位
buffered <- st_buffer(point, dist = 1)
plot(buffered, main = "Buffer around point")

library(sf)

# 创建多边形 (城市边界)
polys <- st_as_sf(data.frame(
    city = c("CityA", "CityB"),
    geometry = st_sfc(
        st_polygon(list(rbind(c(0, 0), c(0, 3), c(3, 3), c(3, 0), c(0, 0)))),
        st_polygon(list(rbind(c(3, 0), c(3, 3), c(6, 3), c(6, 0), c(3, 0))))
    )
))

# 创建点 (学校位置)
points <- st_as_sf(data.frame(
    school = c("School1", "School2"),
    geometry = st_sfc(st_point(c(1, 1)), st_point(c(4, 1)))
))

# 左连接：将每个学校连接到所在的城市
left_join <- st_join(points, polys, left = TRUE)

print(left_join)

# 右连接：保留所有的城市（多边形），即使其中没有学校
right_join <- st_join(polys, points, left = FALSE)

print(right_join)

# 内连接：仅保留点与多边形有交集的部分
inner_join <- st_join(points, polys, left = FALSE)

print(inner_join)

# 创建线 (河流)
river <- st_as_sf(data.frame(name = "River", 
                             geometry = st_sfc(st_linestring(rbind(c(0,0), c(5,5))))))

# 创建点 (学校)
points <- st_as_sf(data.frame(school = c("School1", "School2", "School3"), 
                              geometry = st_sfc(st_point(c(1,1)), st_point(c(4,1)), st_point(c(6,6)))))

# 为河流创建缓冲区 (例如，1单位范围内)
buffered_river <- st_buffer(river, dist = 1)

# 查找哪些学校位于河流缓冲区内
buffer_join <- st_join(points, buffered_river, join = st_within)

print(buffer_join)

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

# values(r) # 可能会存在很多的NA值

r <- rast(ncol = 10, nrow = 10,
          xmin = -150, xmax = -80, ymin = 20, ymax = 60)
r

nrow(r) # number of rows
ncol(r) # number of columns
dim(r) # dimension
ncell(r) # number of cells

values(r) <- 1:ncell(r)

r2 <- r * r
s <- c(r, r2)

plot(s[[2]]) # layer 2

plot(min(s))
plot(r + r + 10)
plot(round(r))
plot(r == 1)

pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)
v
# dimensions：数据包含 12 个几何对象（多边形）和
# 6 个属性字段（如 ID_1, NAME_1, AREA, POP）。

# Longitude and latitude values
long <- c(-0.118092, 2.349014, -3.703339, 12.496366)
lat <- c(51.509865, 48.864716, 40.416729, 41.902782)
longlat <- cbind(long, lat)

# CRS
crspoints <- "+proj=longlat +datum=WGS84"

# Attributes for points
d <- data.frame(
place = c("London", "Paris", "Madrid", "Rome"),
value = c(200, 300, 400, 500))

# SpatVector object
pts <- vect(longlat, atts = d, crs = crspoints)

pts
plot(pts)

geometries <- geom(pts)
print(geometries)

attributes <- values(pts)
print(attributes)

library(terra)
r <- geodata::worldclim_country(country = "Spain", var = "tavg",
                                res = 10, path = tempdir())
plot(r)

dim(r)

r <- mean(r)
dim(r)
plot(r)
ext(r)
r <- terra::mask(r, vect(map))
plot(r)

# Aggregating
r <- terra::aggregate(r, fact = 20, fun = "mean", na.rm = TRUE)
# fact = 20：表示将栅格分辨率减少 20 倍，
# 即每 20x20 的像元块被聚合为一个新的像元。
# fun = "mean"：在聚合过程中，使用平均值作为新的像元值。

# na.rm = TRUE 是一个参数，它的作用是忽略缺失值 (NA)
# 在计算栅格数据的聚合时不将 NA 值计入计算
plot(r)

library(terra)
# Raster (SpatRaster)
r <- rast(system.file("ex/elev.tif", package = "terra"))
# Polygons (SpatVector)
v <- vect(system.file("ex/lux.shp", package = "terra"))

plot(r)
plot(v)

if (FALSE) {
  # 这是多行注释
  # 使用 if (FALSE) 包围的代码块不会被执行
  z <- 30
  print(z)
}

points <- crds(centroids(v))
# 这行代码的作用是计算矢量对象 v 中几何体的质心（centroid），并提取它们的坐标

plot(r)
plot(v, add = TRUE) # add = TRUE 表示在已有的图形上添加新的图层
points(points)  # 这行代码在已经绘制好的栅格图和矢量图上添加质心点

points

# data frame with the coordinates
points <- as.data.frame(points)
points
valuesatpoints <- extract(r, points)
cbind(points, valuesatpoints)

# Extracted raster cells within each polygon
head(extract(r, v, na.rm = TRUE))

# Extracted raster cells and percentage of area
# covered within each polygon
head(extract(r, v, na.rm = TRUE, weights = TRUE))

# Average raster values by polygon
v$avg <- extract(r, v, mean, na.rm = TRUE)$elevation
# Area-weighted average raster values by polygon (weights = TRUE)
v$weightedavg <- extract(r, v, mean, na.rm = TRUE,
                 weights = TRUE)$elevation

print(v)

library(ggplot2)
library(tidyterra)

# Plot average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = avg)) +
  scale_fill_terrain_c()

# Plot area-weighted average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = weightedavg)) +
  scale_fill_terrain_c()

library(terra)

# 创建一个示例栅格数据 SpatRaster
r <- rast(ncol = 10, nrow = 10, xmin = 0, xmax = 10, ymin = 0, ymax = 10)
values(r) <- runif(ncell(r), min = 15, max = 25)  # 随机填充温度值

# 创建一个示例多边形 SpatVector
v <- vect(cbind(c(2, 4, 6, 8), c(2, 4, 6, 8)), type = "polygons")

# 提取多边形内的像元值，并计算平均温度（默认仅提取中心点被覆盖的像元）
extracted_values <- extract(r, v, fun = mean)
head(extracted_values)