2.3 坐标参考系统
空间数据的坐标参考系统（CRS）指定了空间坐标的原点和测量单位。CRS 对于空间数据的操作、分析和可视化非常重要，允许通过将多个数据转换为相同的 CRS 来处理它们。地球上的位置可以使用未投影（也称为地理）或投影的 CRS 进行参照。未投影或地理 CRS 使用纬度和经度值表示地球三维椭球表面上的位置（例如厦门市大致位于东经118度，北纬24.5度）。
投影 CRS 使用笛卡尔坐标在二维平面上展示地球上的位置。所有的投影坐标参考系都会以某种方式扭曲地球表面，无法同时保留面积、方向、形状和距离的所有属性。
最常见的 CRS 可以通过提供 EPSG（欧洲石油勘探集团）代码或 Proj4 字符串来指定。常见的空间投影可以在 https://spatialreference.org/ref/ 找到。例如，EPSG 代码 4326 指的是 WGS84 经纬度坐标系（适用于GPS，有时被称为等经纬度投影，实际是未投影坐标系，但是绘制地图时，会自动应用简单投影或将经纬度假设为平面坐标）。而Proj4 字符串则通过一系列+号相连的键与键值对来制定投影方式的参数（例如名称、区域、基准、距离单位等）。可以使用 sf 包的 st_crs() 函数检查给定坐标系的详细信息。
# 坐标系的名称
st_crs("EPSG:4326")$Name
[1] "WGS 84"
# 坐标系的投影字符串标识
st_crs("EPSG:4326")$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
# 坐标系的EPSG代码
st_crs("EPSG:4326")$epsg
[1] 4326
可以使用 sf 和 terra 对坐标系进行转换。通过 st_crs(x) <- value（如果 x 是 sf 对象）或 crs(r) <- value（如果 r 是栅格对象）可以设置空间数据的 CRS。通过 sf::st_transform() 和 terra::project() 可以转换坐标系。
library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
map <- st_read(pathshp, quiet = TRUE)
# 获取 CRS
st_crs(map)
Coordinate Reference System:
User input: NAD27
wkt:
GEOGCRS["NAD27",
DATUM["North American Datum 1927",
ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4267]]
# 转换 CRS
map2 <- st_transform(map, crs = "EPSG:4326")
# 获取新的 CRS
st_crs(map2)
Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
MEMBER["World Geodetic System 1984 (G2296)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
terra 读取并获取 CRS，转换新的 CRS类似。
library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
# 获取 CRS
crs(r)
[1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
# 转换 CRS
r2 <- terra::project(r, "EPSG:2169")
# 获取 新的CRS
crs(r2)
[1] "PROJCRS[\"LUREF / Luxembourg TM\",\n BASEGEOGCRS[\"LUREF\",\n DATUM[\"Luxembourg Reference Frame\",\n ELLIPSOID[\"International 1924\",6378388,297,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4181]],\n CONVERSION[\"Luxembourg TM\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",49.8333333333333,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",6.16666666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",1,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",80000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",100000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"northing (X)\",north,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"easting (Y)\",east,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Luxembourg.\"],\n BBOX[49.44,5.73,50.19,6.53]],\n ID[\"EPSG\",2169]]"
R语言发展迅速，建议尽量使用新的包，旧的包一般不再维护，安装中容易出错，调试起来也比较麻烦。在 sf 包开发之前，sp 包被用来表示和处理向量空间数据。sp 以及 rgdal、rgeos 和 maptools包已不再维护并将退出使用。