Lambert正形投影转经纬度及空间点距离计算技术问询
2026-6-13
计算跨坐标系点与多边形质心的距离（R语言实现）
问题背景
我是空间分析纯新手，之前看过类似问题的解答但完全看不懂。现在手里有一个经纬度格式的单点，还有一组Lambert正形投影（双标准纬线兰伯特投影）格式的多边形，需要计算多边形质心到这个单点的距离。以下是可复现的初始R代码：
library(sf) #> Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE library(tidygeocoder) library(tidyverse) address <- tibble(street="935 Ramsey Lake Road", city="Sudbury", province="ON", postalcode="P3E 2C6", country="Canada") laurentian <- geocode(.tbl=address, street="street", city="city", state="province", postalcode="postalcode", country="country") #> Passing 1 address to the Nominatim single address geocoder #> Query completed in: 1 seconds laurentian <- st_as_sf(laurentian, coords=c('lat', 'long')) poll <- structure(list(OBJECTID = 63404L, ED_ID = 103, PD_NUMBER = 413, PD_LABEL = "PD 413", ED_NAME_EN = "Sudbury", ED_NAME_FR = "Sudbury", SHAPE_Leng = 452.723045332, SHAPE_Area = 12088.4481776, Year = "2018", geometry = structure(list(structure(c(1229668.80703722, 5795765.97362501 ), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = 1229668.80703722, ymin = 5795765.97362501, xmax = 1229668.80703722, ymax = 5795765.97362501 ), class = "bbox"), crs = structure(list(input = "EO_Lambert_Conformal_Conic", wkt = "PROJCRS[\"EO_Lambert_Conformal_Conic\", BASEGEOGCRS[\"NAD83\", DATUM[\"North American Datum 1983\", ELLIPSOID[\"GRS 1980\",6378137,298.257222101, LENGTHUNIT[\"metre\",1]], ID[\"EPSG\",6269]], PRIMEM[\"Greenwich\",0, ANGLEUNIT[\"Degree\",0.0174532925199433]]], CONVERSION[\"unnamed\", METHOD[\"Lambert Conic Conformal (2SP)\", ID[\"EPSG\",9802]], PARAMETER[\"Latitude of false origin\",0, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8821]], PARAMETER[\"Longitude of false origin\",-84, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8822]], PARAMETER[\"Latitude of 1st standard parallel\",44.5, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8823]], PARAMETER[\"Latitude of 2nd standard parallel\",54.5, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8824]], PARAMETER[\"Easting at false origin\",1000000, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8826]], PARAMETER[\"Northing at false origin\",0, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8827]]], CS[Cartesian,2], AXIS[\"(E)\",east, ORDER[1], LENGTHUNIT[\"metre\",1, ID[\"EPSG\",9001]]], AXIS[\"(N)\",north, ORDER[2], LENGTHUNIT[\"metre\",1, ID[\"EPSG\",9001]]]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, -1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(OBJECTID = NA_integer_, ED_ID = NA_integer_, PD_NUMBER = NA_integer_, PD_LABEL = NA_integer_, ED_NAME_EN = NA_integer_, ED_NAME_FR = NA_integer_, SHAPE_Leng = NA_integer_, SHAPE_Area = NA_integer_, Year = NA_integer_), class = "factor", levels = c("constant", "aggregate", "identity"))) poll <- st_centroid(poll)
解决步骤
核心问题是两个空间对象的坐标系不匹配，必须先统一坐标系才能计算有效距离，具体操作如下：
1. 修正经纬度点的坐标顺序
st_as_sf()默认是经度（x轴）在前，纬度（y轴）在后，初始代码里写反了，先修正：
laurentian <- st_as_sf(laurentian, coords=c('long', 'lat'))
2. 统一坐标系
将经纬度点转换为和多边形一致的Lambert投影，这样计算的距离单位是米，结果更直观：
# 提取多边形的坐标系参数，转换经纬度点 laurentian_proj <- st_transform(laurentian, st_crs(poll))
3. 计算两点距离
用st_distance()直接计算转换后的两点距离：
# 计算距离，结果带单位（米） distance <- st_distance(laurentian_proj, poll) # 查看结果 print(distance) # 如果需要纯数值，可转换 distance_num <- as.numeric(distance)
完整整合代码
library(sf) library(tidygeocoder) library(tidyverse) # 地址解析得到经纬度点 address <- tibble(street="935 Ramsey Lake Road", city="Sudbury", province="ON", postalcode="P3E 2C6", country="Canada") laurentian <- geocode(.tbl=address, street="street", city="city", state="province", postalcode="postalcode", country="country") # 修正坐标顺序，转为sf对象 laurentian <- st_as_sf(laurentian, coords=c('long', 'lat')) # 加载多边形数据并计算质心 poll <- structure(list(OBJECTID = 63404L, ED_ID = 103, PD_NUMBER = 413, PD_LABEL = "PD 413", ED_NAME_EN = "Sudbury", ED_NAME_FR = "Sudbury", SHAPE_Leng = 452.723045332, SHAPE_Area = 12088.4481776, Year = "2018", geometry = structure(list(structure(c(1229668.80703722, 5795765.97362501 ), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = 1229668.80703722, ymin = 5795765.97362501, xmax = 1229668.80703722, ymax = 5795765.97362501 ), class = "bbox"), crs = structure(list(input = "EO_Lambert_Conformal_Conic", wkt = "PROJCRS[\"EO_Lambert_Conformal_Conic\", BASEGEOGCRS[\"NAD83\", DATUM[\"North American Datum 1983\", ELLIPSOID[\"GRS 1980\",6378137,298.257222101, LENGTHUNIT[\"metre\",1]], ID[\"EPSG\",6269]], PRIMEM[\"Greenwich\",0, ANGLEUNIT[\"Degree\",0.0174532925199433]]], CONVERSION[\"unnamed\", METHOD[\"Lambert Conic Conformal (2SP)\", ID[\"EPSG\",9802]], PARAMETER[\"Latitude of false origin\",0, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8821]], PARAMETER[\"Longitude of false origin\",-84, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8822]], PARAMETER[\"Latitude of 1st standard parallel\",44.5, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8823]], PARAMETER[\"Latitude of 2nd standard parallel\",54.5, ANGLEUNIT[\"Degree\",0.0174532925199433], ID[\"EPSG\",8824]], PARAMETER[\"Easting at false origin\",1000000, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8826]], PARAMETER[\"Northing at false origin\",0, LENGTHUNIT[\"metre\",1], ID[\"EPSG\",8827]]], CS[Cartesian,2], AXIS[\"(E)\",east, ORDER[1], LENGTHUNIT[\"metre\",1, ID[\"EPSG\",9001]]], AXIS[\"(N)\",north, ORDER[2], LENGTHUNIT[\"metre\",1, ID[\"EPSG\",9001]]]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, -1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(OBJECTID = NA_integer_, ED_ID = NA_integer_, PD_NUMBER = NA_integer_, PD_LABEL = NA_integer_, ED_NAME_EN = NA_integer_, ED_NAME_FR = NA_integer_, SHAPE_Leng = NA_integer_, SHAPE_Area = NA_integer_, Year = NA_integer_), class = "factor", levels = c("constant", "aggregate", "identity"))) poll <- st_centroid(poll) # 统一坐标系 laurentian_proj <- st_transform(laurentian, st_crs(poll)) # 计算距离 distance <- st_distance(laurentian_proj, poll) print(distance)
关键提示
- 坐标系统一是空间计算的核心前提，不同坐标系下直接计算距离没有意义。
st_transform()会自动根据两个坐标系的参数完成投影转换，不需要手动推导公式。st_distance()返回的是带单位的结果，若需要纯数值可使用as.numeric()转换。
内容的提问来源于stack exchange，提问作者spindoctor
相关产品推荐
相关产品推荐
