高效计算地理距离：3万条数据下30英里半径地点数量统计

刚处理过类似的3万级地点数据需求，来给你捋捋可行的方案——首先你提到的小数据集方法确实能用，但直接套在3万条数据上会炸，先给你补全小数据集的完整代码，再重点说大数据量的优化方案。

小数据集的基础实现（不适合3万条） #

如果数据量小（比如几千条），用PROC SQL自连接+GEODIST函数计算距离是最直接的，完整代码如下：

PROC SQL;
CREATE TABLE Distance_Table_1 AS
SELECT 
    MASTER.PlaceID AS PlaceID,
    Master.INTPTLAT AS LAT1,
    Master.INTPTLONG AS LONG1,
    Match.INTPTLAT AS LAT2,
    Match.INTPTLONG AS LONG2,
    -- 用'MILE'参数指定计算英里距离
    GEODIST(Master.INTPTLAT, Master.INTPTLONG, Match.INTPTLAT, Match.INTPTLONG, 'MILE') AS Distance
FROM 
    your_dataset_name AS MASTER
INNER JOIN 
    your_dataset_name AS Match
ON 
    MASTER.PlaceID NE Match.PlaceID; -- 排除地点自身的匹配
QUIT;

-- 统计每个地点30英里内的数量
PROC SQL;
CREATE TABLE Final_Count AS
SELECT 
    PlaceID,
    COUNT(*) AS Within_30Mile_Count
FROM 
    Distance_Table_1
WHERE 
    Distance <= 30
GROUP BY 
    PlaceID;
QUIT;

但注意：3万条数据自连接会生成30000×29999≈9亿条记录，SAS处理起来会慢到离谱，甚至会因为内存不足报错，所以必须优化。

3万条数据的优化方案 #

方案1：先过滤经纬度范围，减少无效计算 #

30英里对应的经纬度范围大概是：

• 纬度：1度≈69英里，所以30英里≈±0.45度
• 经度：不同纬度下差异大，比如北纬40度，1度经度≈53英里，所以取±0.6度的范围（留冗余避免漏判）

先通过经纬度范围过滤掉肯定不在30英里内的点，再计算距离，能大幅减少连接数：

PROC SQL;
CREATE TABLE Distance_Table_Optimized AS
SELECT 
    MASTER.PlaceID AS PlaceID,
    Master.INTPTLAT AS LAT1,
    Master.INTPTLONG AS LONG1,
    Match.INTPTLAT AS LAT2,
    Match.INTPTLONG AS LONG2,
    GEODIST(Master.INTPTLAT, Master.INTPTLONG, Match.INTPTLAT, Match.INTPTLONG, 'MILE') AS Distance
FROM 
    your_dataset_name AS MASTER
INNER JOIN 
    your_dataset_name AS Match
ON 
    MASTER.PlaceID NE Match.PlaceID
    -- 先过滤经纬度范围，缩小计算范围
    AND Match.INTPTLAT BETWEEN Master.INTPTLAT - 0.45 AND Master.INTPTLAT + 0.45
    AND Match.INTPTLONG BETWEEN Master.INTPTLONG - 0.6 AND Master.INTPTLONG + 0.6;
QUIT;

-- 后续统计和基础版一致
PROC SQL;
CREATE TABLE Final_Count_Optimized AS
SELECT 
    PlaceID,
    COUNT(*) AS Within_30Mile_Count
FROM 
    Distance_Table_Optimized
WHERE 
    Distance <= 30
GROUP BY 
    PlaceID;
QUIT;

这个方案不需要额外模块，3万条数据跑起来应该在几十分钟内完成（具体看服务器性能）。

方案2：用SAS空间模块创建索引（更快） #

如果你的SAS有SAS/GIS或SAS/GRAPH模块，用空间索引做邻近查询效率会高很多，代码示例：

-- 先给经纬度创建空间索引
PROC SPATIAL DATA=your_dataset_name;
    CREATE INDEX spatial_idx ON your_dataset_name(INTPTLAT, INTPTLONG) TYPE=GEODIST;
RUN;

-- 直接查找每个点30英里内的所有邻近点
PROC SPATIAL DATA=your_dataset_name OUT=Neighbors;
    NEAR NEIGHBORS=ALL DISTANCE=30 UNITS=MILE;
    ID PlaceID; -- 指定用PlaceID标识每个地点
RUN;

-- 统计数量，记得排除自身匹配
PROC SQL;
CREATE TABLE Final_Count_Spatial AS
SELECT 
    PlaceID,
    COUNT(*) AS Within_30Mile_Count
FROM 
    Neighbors
WHERE 
    PlaceID NE NeighborID
GROUP BY 
    PlaceID;
QUIT;

这个方法利用空间索引的特性，能把计算量从9亿级降到百万级，跑起来会快很多。

最后提醒下：确保你的经纬度是十进制数值型（比如北纬40.7128，西经-74.0060），GEODIST函数不支持度分秒格式的输入哦。

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