坐标参考系与投影
使用 st_crs() 和 st_transform() 正确地在 CRS 之间转换
坐标参考系与投影 是 CoddyKit 上的免费 R Academy 课时。 这是第 2 节课,共 4 节。 你可以在下方免费阅读本课时的完整内容 — 然后在浏览器中使用内置代码编辑器和全天候 AI 导师进行实践。 这是 R Academy 学习路径的一部分,你的进度在网页和 CoddyKit 应用中同步。 R Academy 课程共包含 4 节课。
CRS 为什么重要
坐标参考系(CRS)告诉 R,几何列中的数字代表什么——纬度/经度、投影网格中的米,或其他含义。如果不进行重新投影就混用不同 CRS,会导致空间连接静默失败,并使距离计算结果错误。
library(sf)
# Two sf objects in the same CRS: operations work
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# sf will error or warn when CRS do not match during operations
# This is the CRS metadata attached to the object
crs <- st_crs(nc)
cat('EPSG:', crs$epsg, '\n')
cat('Input is geographic (lon/lat)?', st_is_longlat(nc), '\n')
cat('Units:', crs$units_gdal, '\n')st_crs:读取 CRS 信息
st_crs(sf_obj) 会返回一个 crs 对象,其中包含 EPSG 代码、WKT 字符串、Proj4 字符串和单位信息。您也可以直接传入整数 EPSG 代码,以获取 CRS 定义。
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# Inspect the CRS object
crs <- st_crs(nc)
cat('Class:', class(crs), '\n')
cat('EPSG:', crs$epsg, '\n')
cat('Proj4:', substr(crs$proj4string, 1, 60), '...\n')
# Lookup a CRS by EPSG code
wgs84 <- st_crs(4326)
web_mercator <- st_crs(3857)
cat('\nWGS84 name:', wgs84$Name, '\n')
cat('Web Mercator name:', web_mercator$Name, '\n')EPSG:4326——WGS84
EPSG:4326 是1984 世界大地测量系统,也是 GPS 坐标的标准。坐标使用十进制度数表示:经度(X)范围为 −180 到 180,纬度(Y)范围为 −90 到 90。它是地理 CRS,而不是投影 CRS。
library(sf)
# Create points in WGS84 (GPS coordinates)
cities <- data.frame(
city = c('London', 'Tokyo', 'Sydney', 'New York'),
lon = c(-0.118, 139.692, 151.209, -74.006),
lat = c(51.509, 35.689, -33.869, 40.713)
)
cities_sf <- st_as_sf(cities,
coords = c('lon', 'lat'),
crs = 4326
)
cat('CRS:', st_crs(cities_sf)$Name, '\n')
cat('Units:', st_crs(cities_sf)$units_gdal, '\n')
# Naive distance in degrees (not meaningful for real distances)
d <- st_distance(cities_sf[1, ], cities_sf[2, ])
cat('London-Tokyo distance (degrees):', round(as.numeric(d), 2), '\n')EPSG:3857——Web Mercator
EPSG:3857 是 Google Maps、OpenStreetMap 和 Leaflet 使用的 Web Mercator 投影。坐标单位为米,以本初子午线和赤道为中心。它会使极地附近的面积产生变形,但适合用于网页可视化。
library(sf)
cities <- data.frame(
city = c('London', 'Tokyo'),
lon = c(-0.118, 139.692),
lat = c(51.509, 35.689)
)
cities_sf <- st_as_sf(cities, coords = c('lon', 'lat'), crs = 4326)
# Transform to Web Mercator
cities_merc <- st_transform(cities_sf, 3857)
cat('WGS84 coords (degrees):\n')
print(st_coordinates(cities_sf))
cat('\nWeb Mercator coords (metres):\n')
print(st_coordinates(cities_merc))
# Distance is now in metres
d <- st_distance(cities_merc[1, ], cities_merc[2, ])
cat('\nLondon-Tokyo distance (km):', round(as.numeric(d) / 1000, 0), '\n')st_set_crs:指定 CRS
st_set_crs(sf_obj, crs) 会在不转换坐标的情况下指定 CRS。仅当 CRS 元数据缺失或错误时才应使用它。如果您想将坐标重新投影到新的 CRS,请改用 st_transform()。
library(sf)
# Create geometry with NO CRS
pts <- st_sfc(
st_point(c(-74.006, 40.713)),
st_point(c(-87.629, 41.878))
)
cat('CRS before:', is.na(st_crs(pts)), '(NA = no CRS)\n')
# Assign WGS84 (no coordinate transformation)
pts_wgs84 <- st_set_crs(pts, 4326)
cat('CRS after assignment:', st_crs(pts_wgs84)$epsg, '\n')
# Compare with st_transform (DOES transform coordinates)
pts_sf <- st_sf(id = 1:2, geometry = pts_wgs84)
pts_merc <- st_transform(pts_sf, 3857)
cat('Original lon:', st_coordinates(pts_sf)[1, 'X'], '\n')
cat('Mercator X: ', st_coordinates(pts_merc)[1, 'X'], '\n')st_transform:重新投影
st_transform(sf_obj, crs) 会将所有几何对象从当前 CRS 重新投影到目标 CRS,并执行完整的大地测量转换。crs 参数可以接受 EPSG 整数、Proj4 字符串或 WKT 字符串。
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
cat('Original CRS:', st_crs(nc)$epsg, '\n')
# Re-project to WGS84
nc_wgs84 <- st_transform(nc, 4326)
cat('Transformed CRS:', st_crs(nc_wgs84)$epsg, '\n')
# Re-project to a local projected CRS (NC State Plane)
nc_stateplane <- st_transform(nc, 32119) # NAD83 / North Carolina
cat('State Plane CRS:', st_crs(nc_stateplane)$Name, '\n')
# Area calculation is more accurate in a projected CRS
area_deg <- st_area(nc[1, ])
area_m2 <- st_area(nc_stateplane[1, ])
cat('Area (degrees^2):', format(area_deg, big.mark = ','), '\n')
cat('Area (m^2): ', format(area_m2, big.mark = ','), '\n')选择合适的投影
不同投影会针对不同属性进行优化:等积投影(例如 Mollweide)、保形投影(Mercator)和等距投影。进行局部分析时,请使用以米为单位的投影 CRS;查看全球概览时,WGS84 通常适用;制作网页地图时,请使用 Web Mercator。
library(sf)
# Common EPSG codes cheat sheet
epsg_table <- data.frame(
EPSG = c(4326, 3857, 4269, 32637, 27700),
Name = c(
'WGS84 (GPS, global)',
'Web Mercator (web maps)',
'NAD83 (North America geographic)',
'UTM zone 37N (Middle East)',
'British National Grid (UK)'
),
Units = c('degrees', 'metres', 'degrees', 'metres', 'metres')
)
print(epsg_table, row.names = FALSE)
# Look up any EPSG online: https://epsg.io/<code>
cat('\nTip: use st_crs(epsg_code)$Name to verify a code\n')Proj4 字符串
在 EPSG 代码成为标准之前,CRS 由 Proj4 字符串定义,这些是传递给 PROJ 库的紧凑参数字符串。现在仍然可以将它们与 st_transform() 一起使用,但为了清晰明了并兼容未来版本,优先使用 EPSG 代码。
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# WGS84 as Proj4 string
wgs84_proj4 <- '+proj=longlat +datum=WGS84 +no_defs'
nc_wgs84 <- st_transform(nc, crs = wgs84_proj4)
cat('Using Proj4 string:\n', wgs84_proj4, '\n')
cat('Result CRS EPSG:', st_crs(nc_wgs84)$epsg, '\n')
# Compare: equivalent to st_transform(nc, 4326)
cat('Same as EPSG 4326?',
identical(st_crs(nc_wgs84), st_crs(st_transform(nc, 4326))), '\n')检查 CRS 是否相同
在执行任何空间连接或叠加操作之前,请使用 st_crs(a) == st_crs(b) 验证两个对象是否共享相同的 CRS。如果不同,请先将其中一个转换为与另一个相同的 CRS,再继续操作。
library(sf)
# Two objects in different CRS
nc_nad27 <- read_sf(system.file('shape/nc.shp', package = 'sf'))
nc_wgs84 <- st_transform(nc_nad27, 4326)
cat('Same CRS?', st_crs(nc_nad27) == st_crs(nc_wgs84), '\n')
# Always check before join/overlay
align_crs <- function(a, b) {
if (st_crs(a) != st_crs(b)) {
message('CRS mismatch! Transforming b to match a.')
b <- st_transform(b, st_crs(a))
}
list(a = a, b = b)
}
aligned <- align_crs(nc_nad27, nc_wgs84)
cat('After alignment - same CRS?',
st_crs(aligned$a) == st_crs(aligned$b), '\n')使用 CRS 计算距离
在地理 CRS(单位为度)上计算的距离没有实际意义。请先重新投影到等距或等积投影 CRS,或者使用 sf 内置的大地测量距离计算功能;它通过 PROJ 库正确处理 WGS84。
library(sf)
# City points in WGS84
cities <- data.frame(
city = c('Paris', 'Berlin', 'Madrid', 'Rome'),
lon = c(2.350, 13.404, -3.702, 12.496),
lat = c(48.865, 52.520, 40.417, 41.902)
)
cities_sf <- st_as_sf(cities, coords = c('lon', 'lat'), crs = 4326)
# sf computes geodesic distances on WGS84 automatically
dist_matrix <- st_distance(cities_sf)
rownames(dist_matrix) <- cities$city
colnames(dist_matrix) <- cities$city
# Convert to km
dist_km <- round(as.numeric(dist_matrix) / 1000)
mat_km <- matrix(dist_km, 4, 4, dimnames = list(cities$city, cities$city))
print(mat_km)完整的 CRS 工作流程
下面是一套完整的工作流程:加载原始 GPS 数据,指定为 WGS84,重新投影到本地米制 CRS 以准确计算面积和距离,执行分析,然后重新投影回 WGS84,以便输出到网页地图。
library(sf)
# Load and assign CRS
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
nc_wgs84 <- st_transform(nc, 4326)
# Re-project to NC State Plane (metres) for accurate areas
nc_proj <- st_transform(nc, 32119)
# Compute area per county in km2
nc_proj$area_km2 <- as.numeric(st_area(nc_proj)) / 1e6
# Top 5 largest counties
top5 <- nc_proj[order(-nc_proj$area_km2), c('NAME', 'area_km2')]
cat('5 Largest NC Counties (km2):\n')
print(head(as.data.frame(top5)[, c('NAME', 'area_km2')], 5))
# Back to WGS84 for output
nc_output <- st_transform(nc_proj, 4326)
cat('\nOutput CRS:', st_crs(nc_output)$epsg, '\n')快速检查
您需要为当前处于 WGS84(EPSG:4326)中的多边形要素计算精确面积(单位为平方米)。您应该怎么做?
回顾:CRS 与投影
要点:
st_crs(sf_obj)可读取 CRS 元数据;指定 CRS 时优先使用 EPSG 代码- EPSG:4326 = WGS84(用于 GPS 的度数);EPSG:3857 = Web Mercator(用于网页地图的米)
st_set_crs()会指定 CRS 但不进行转换;st_transform()会重新投影坐标- 执行空间操作前,始终使用
st_crs(a) == st_crs(b)验证 CRS 是否相同 - 要准确计算距离和面积,请先重新投影到米制投影 CRS
- Proj4 字符串仍然可用,但 EPSG 代码更简洁、可读性更好
library(sf)
nc <- read_sf(system.file('shape/nc.shp', package = 'sf'))
# CRS pipeline
cat('Original:', st_crs(nc)$epsg, '\n')
nc_wgs84 <- st_transform(nc, 4326)
nc_proj <- st_transform(nc_wgs84, 32119)
cat('WGS84:', st_crs(nc_wgs84)$epsg, '\n')
cat('Projected:', st_crs(nc_proj)$epsg, '\n')用 AI 导师学习 R — 免费
在浏览器中编写并运行真实代码,获得全天候 AI 导师的即时帮助,并在网页或应用中继续学习。
- 课程
- 43
- 课程
- 159
常见问题解答
「坐标参考系与投影」课时是免费的吗?
是的 — 「坐标参考系与投影」的完整文本可在网页上免费阅读。要进行交互式练习(内置代码编辑器和全天候 AI 导师)并解锁 R Academy 课程的其余内容,请升级到 CoddyKit PRO。 R Academy 课程共包含 4 节课。
「坐标参考系与投影」这节课中我会学到什么?
使用 st_crs() 和 st_transform() 正确地在 CRS 之间转换 你通过在浏览器中直接运行的动手代码来练习 R Academy,全天候 AI 导师会在你学习这节课的过程中回答你的问题。
学习 R Academy 需要有经验吗?
无需任何先前经验。CoddyKit 上的 R Academy 课程适合初学者到高级学习者,你可以从这里开始或从头开始,按照自己的节奏学习。 这是第 2 节课,共 4 节。
「坐标参考系与投影」课时需要多长时间?
大多数 CoddyKit 课程大约需要 5–10 分钟。每节课都很精短且互动,所以你能稳步进步,并在网页和应用中从离开的地方继续。
我能在这节 R Academy 课中编写并运行代码吗?
能。每节 R Academy 课都包含内置代码编辑器,你可以在浏览器中直接编写并运行真实代码,并获得即时 AI 反馈 — 无需本地设置。