如何在R中通过FreeTDS驱动正确读取SQL Server空间数据？

核心问题分析 #

你遇到的问题主要来自两个点：一是FreeTDS将SQL Server的geometry列返回为blob类型，st_read无法直接解析；二是GDAL尝试查询PostGIS风格的spatial_ref_sys表（SQL Server没有这个表），导致SRID识别错误（把25832错判成了70647808）。下面分步骤解决：

方法1：手动转换Blob为SF对象（最可靠） #

绕过st_read的自动解析，直接用SQL Server的空间函数返回标准WKB格式，再用sf包手动转换：

library(odbc)
library(DBI)
library(sf)

# 建立数据库连接
username <- 'SO-user'
sql_server_ip <- '35.214.169.110'
password <- 'SQLaskingfortrouble?!'
con <- dbConnect(odbc(), 
                 Driver = "FreeTDS", 
                 Server = sql_server_ip, 
                 Database = "stackoverflow-example", 
                 UID = username, 
                 PWD = password, 
                 Port = 1433)

# 查询时返回WKB格式的几何数据+SRID
data_raw <- dbGetQuery(con, '
  SELECT TOP 2 
    objectid, 
    geom.STAsBinary() AS geom_wkb,  -- 转换为标准WKB
    geom.STSrid AS srid 
  FROM spatialdata
')

# 转换为SF数据框
data_sf <- st_sf(
  objectid = data_raw$objectid,
  geometry = st_as_sfc(data_raw$geom_wkb),  # 解析WKB
  crs = data_raw$srid[1]  # 应用SRID（假设所有记录SRID一致）
)

# 验证结果
class(data_sf)
st_crs(data_sf)

这个方法完全避开了GDAL对SQL Server空间系统的依赖，兼容性最好。

方法2：配置st_read适配SQL Server #

如果想直接用st_read，需要指定GDAL的layer_options来适配SQL Server的空间格式，同时手动指定CRS避免查询不存在的spatial_ref_sys表：

data_sf <- st_read(
  con,
  query = "SELECT TOP 10 objectid, geom.STAsBinary() AS geom FROM spatialdata",
  geometry_column = "geom",
  crs = 25832,  # 直接指定EPSG编码
  layer_options = "SQLSERVER_SPATIAL=YES"  # 告诉GDAL这是SQL Server空间数据
)

这里的关键是用geom.STAsBinary()返回WKB，让st_read能直接解析，同时手动设置CRS跳过GDAL的SRID查询。

附加问题：高效添加边界框过滤 #

在SQL Server中做空间边界框查询，推荐用STFilter（利用空间索引，效率远高于STIntersects），步骤如下：

# 定义你的边界框（示例：EPSG 25832下的坐标范围）
bbox <- st_bbox(
  c(xmin = 500000, ymin = 6000000, xmax = 510000, ymax = 6010000),
  crs = st_crs(25832)
)

# 转换为SQL Server能识别的WKT格式
bbox_wkt <- st_as_text(st_as_sfc(bbox))

# 构造带边界框过滤的查询
bbox_query <- sprintf("
  SELECT objectid, geom.STAsBinary() AS geom 
  FROM spatialdata 
  WHERE geom.STFilter(geometry::STGeomFromText('%s', 25832)) IS NOT NULL
", bbox_wkt)

# 读取过滤后的空间数据
data_bbox_sf <- st_read(con, query = bbox_query, geometry_column = "geom", crs = 25832)

如果你的表没有创建空间索引，先执行这个SQL创建索引（大幅提升查询速度）：

CREATE SPATIAL INDEX SI_spatialdata_geom ON spatialdata(geom)

内容的提问来源于stack exchange，提问作者fifthace