让 AI 帮我做空间分析,结果直接翻车了
上周我用 Cursor 帮我写了个面积统计脚本跑完一看结果——差了 10 倍。查了半天发现是坐标系的问题。AI 生成的代码没问题但它不懂“空间思维”有些坑它不会主动提醒你。今天把这次翻车经历拆开讲顺便总结几个 AI 做 GIS 时的高频坑。一、事情经过起因我统计各地块的面积按地类汇总输出 Excel。文件打开看了一眼坐标系GCS_WGS_1984地理坐标系要素数量2,856 个范围112.3°E ~ 112.8°E, 37.4°N ~ 37.9°N太原市周边的数据坐标系是 WGS84 地理坐标。我让 AI 写的代码直接把需求丢给 Cursor帮我写一个 Python 脚本1. 读取 Shapefile2. 计算每个要素的面积3. 按地类字段DLBM汇总面积4. 输出到 ExcelCursor 3 秒就给了代码import geopandas as gpdimport pandas as pd# 读取数据gdf gpd.read_file(“地块.shp”)# 计算面积gdf[“面积”] gdf.geometry.area# 按地类汇总summary gdf.groupby(“DLBM”)[“面积”].sum().reset_index()summary.columns [“地类代码”, “面积”]# 输出summary.to_excel(“面积统计.xlsx”, indexFalse)print(“完成”)代码看起来没问题我直接跑了。结果输出 Excel 一看地类代码面积010.000312020.000458030.000291......面积全是 0.000x 这种小数。我愣了一下——甲方说这块区域大概 50 平方公里怎么加起来才 0.00x二、排查过程第一反应单位问题我看了一眼数据坐标系是 WGS84单位是度。geometry.area 算出来的是平方度不是平方米。但平方度怎么换算我查了一下1 度 ≈ 111 km赤道1 平方度 ≈ 12,321 平方公里那我这个 0.003 平方度换算成平方公里应该是0.003 × 12,321 ≈ 37 平方公里好像对上了不对这个换算只在赤道成立。真正的问题纬度修正地球是椭球经度 1 度对应的实际距离随纬度变化经度 1 度 ≈ 111 km × cos(纬度)太原纬度约 37.5°cos(37.5°) ≈ 0.79。所以在这个纬度1 平方度 ≈ 111 × 111 × 0.79 ≈ 9,700 平方公里但这个换算太粗糙了真正准确的做法是投影到平面坐标系。正确做法import geopandas as gpdgdf gpd.read_file(“地块.shp”)# 关键一步投影到平面坐标系# 太原在 UTM 49N 区EPSG:32649gdf_proj gdf.to_crs(epsg32649)# 现在面积单位是平方米gdf_proj[“面积_平方米”] gdf_proj.geometry.area# 换算成公顷gdf_proj[“面积_公顷”] gdf_proj[“面积_平方米”] / 10000# 汇总summary gdf_proj.groupby(“DLBM”)[“面积_公顷”].sum().reset_index()投影之后面积单位才是米结果才准确。三、AI 为什么没提醒我我回头问 Cursor“你为什么不提醒我要先投影”它回答“抱歉你的需求里没有提到坐标系我假设数据已经是合适的坐标系。”问题就在这里。AI 懂 Python懂 geopandas 语法但它不懂 GIS 的空间思维。它不知道- 地理坐标系的单位是度不能直接算面积- 不同纬度1 度对应的实际距离不同- 面积/距离计算必须用投影坐标系这些是 GIS 的领域知识不是编程知识。AI 没学过。四、AI 做 GIS 的 9 个高频坑这次翻车之后我总结了 AI 容易踩的 9 个坑。实际上这些坑我们在用传统方法的时候也会犯但传统方法直面数据结构和属性能够直观地找到错误本身用 AI 后这些问题就会被掩盖不仔细认真很难发现坑 1坐标系不检查表现AI 直接对地理坐标系做面积/距离计算后果结果单位是“平方度”或“度”数值差 10~100 倍正确做法# 先检查坐标系print(gdf.crs) # 看一眼是地理坐标还是投影坐标# 如果是地理坐标EPSG:4326先投影if gdf.crs.to_epsg() 4326:gdf gdf.to_crs(epsg合适的目标坐标系)坑 2投影坐标系选错表现AI 随便选一个投影比如全用 UTM不管数据在哪后果投影变形大面积/距离不准正确做法- 中国境内用 CGCS2000 投影如 EPSG:4547 是 3°带第 47 区- 或者根据数据范围自动选 UTM 区from pyproj import CRS# 根据中心经度选 UTM 区lon gdf.geometry.centroid.x.mean()utm_zone int((lon 180) / 6) 1utm_epsg 32600 utm_zone # 北半球坑 3Shapefile 字段名截断表现AI 用中文字段名或字段名超过 10 字符后果导出 Shapefile 时字段名被截断数据丢失正确做法# 导出前检查字段名长度for col in gdf.columns:if len(col) 10:print(f“警告字段 {col} 超过 10 字符将被截断”)# 用英文字段名gdf gdf.rename(columns{“地类代码”: “DLBM”, “面积”: “AREA”})坑 4几何有效性不检查表现AI 直接做叠加分析不检查几何是否有效后果报错 TopologyException或者结果有空洞正确做法# 检查几何有效性invalid gdf[~gdf.geometry.is_valid]if len(invalid) 0:print(f“发现 {len(invalid)} 个无效几何正在修复...”)gdf.geometry gdf.geometry.make_valid()坑 5大文件直接读进内存表现AI 用 gpd.read_file() 读几 GB 的文件后果内存爆掉程序崩溃正确做法- 用分块读取gpd.read_file(..., rows10000)- 或用 arcpy 的游标arcpy.da.SearchCursor- 或转成 GeoPackage 分层存储坑 6坐标顺序搞反经度/纬度 vs 纬度/经度表现AI 写 Point(lat, lon) 或 Point(y, x)顺序反了后果点跑到地球另一边或者报错“坐标超出范围”正确做法# GeoPandas / Shapely 标准顺序经度在前纬度在后from shapely.geometry import Point# 正确经度在前x, ypoint Point(112.5, 37.8) # 太原附近# 错误纬度在前# point Point(37.8, 112.5) # 跑到别的地方了# 检查坐标是否在合理范围lon, lat point.x, point.yif not (-180 lon 180 and -90 lat 90):print(“警告坐标可能反了”)我踩过一次200 个点经纬度列反了结果全跑到太平洋里去了。坑 7高德/百度坐标直接叠加标准底图表现AI 用高德 POI 坐标直接叠加到 OSM/天地图底图后果点位偏移几百米看起来“对不上”原因高德/百度用的是 GCJ-02/BD-09 坐标系不是 WGS84正确做法# 需要手动实现坐标转换函数import mathdef gcj02_to_wgs84(lon, lat):GCJ-02 转换为 WGS84a 6378245.0 # 长半轴ee 0.00669342162296594323 # 偏心率平方if _out_of_china(lon, lat):return lon, latdlat _transform_lat(lon - 105.0, lat - 35.0)dlon _transform_lon(lon - 105.0, lat - 35.0)radlat lat / 180.0 * math.pimagic math.sin(radlat)magic 1 - ee * magic * magicsqrtmagic math.sqrt(magic)dlat (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * math.pi)dlon (dlon * 180.0) / (a / sqrtmagic * math.cos(radlat) * math.pi)return lon - dlon, lat - dlat# 使用wgs_lon, wgs_lat gcj02_to_wgs84(gcj_lon, gcj_lat)# 然后才能叠加到标准底图gdf gpd.GeoDataFrame(df,geometrygpd.points_from_xy(df[“wgs_lon”], df[“wgs_lat”]),crs“EPSG:4326”)提醒国内地图服务分三派——- 天地图/OSMWGS84- 高德/腾讯GCJ-02- 百度BD-09混用之前必须转换。坑 8缓冲区距离单位搞错表现AI 写 gdf.buffer(1000)以为生成 1 公里缓冲区后果如果数据是地理坐标系参数单位是“度”而不是“米”生成的缓冲区和你想的完全不一样正确做法# 先确认坐标系单位print(gdf.crs.axis_info[0].unit_name) # “degree” 还是 “metre”# 如果是地理坐标先投影gdf_proj gdf.to_crs(epsg32649) # 投影后单位是米# 再做缓冲区buffer gdf_proj.buffer(1000) # 1000 米 1 公里我踩过一次做“河流 500 米缓冲区”忘了看坐标系单位数据是 WGS84 地理坐标——结果生成的不是 500 米的缓冲区偏差极大。坑 9ArcGIS Pro 和 ArcMap 的 API 混用表现AI 用 arcpy.mappingArcMap API但你用的是 ArcGIS Pro后果报错 AttributeError: module arcpy has no attribute mapping原因ArcGIS Pro 用 arcpy.mp新 APIArcMap 用 arcpy.mapping旧 API正确做法# ArcGIS Pro3.ximport arcpy.mp as mp # 注意是 mp 不是 mappingaprx mp.ArcGISProject(CURRENT)m aprx.listMaps()[0]lyr m.listLayers()[0]# ArcMap10.x— 已过时# import arcpy.mapping as mapping# mxd mapping.MapDocument(CURRENT)提醒告诉 AI 你用的是 ArcGIS Pro 3.x别让它生成 ArcMap 的代码。五、怎么避免这些坑方法 1把领域知识写进 Skill我现在用 WorkBuddy把 GIS 的坑都写进了 gis-assistant Skill 里## 数据质量检查分析前必做- 坐标系是否定义地理坐标不能直接算面积- 几何是否有效无效几何会导致分析失败- 字段名是否超长Shapefile 限 10 字符- 编码是否正确GBK/UTF-8 乱码问题AI 加载这个 Skill 之后每次分析前会自动提醒我检查这些项。方法 2建立检查清单每次拿到数据先跑一遍def quick_check(gdf):print(f“坐标系: {gdf.crs}”)print(f“是否地理坐标: {gdf.crs.is_geographic}”)print(f“要素数量: {len(gdf)}”)print(f“无效几何: {(~gdf.geometry.is_valid).sum()}”)print(f“字段名超长: {sum(len(c)10 for c in gdf.columns)}”)跑完再决定用什么分析方法。方法 3让 AI 解释代码拿到 AI 生成的代码多问一句“这段代码在什么情况下会出错”或者“坐标系是 WGS84 时这个计算结果是什么单位”逼 AI 去思考领域问题而不是只给语法正确的代码。六、这次翻车的教训1. AI 懂编程不懂 GIS——语法正确 ≠ 结果正确2. 坐标系是 GIS 的命根子——面积/距离计算前必须确认3. 领域知识要固化——写进 Skill让 AI 每次都带着4. 结果要校验——算出来先看数量级对不对别直接交5. 坑会重复出现——踩过的坑要记下来形成检查清单现在我用 AI 做 GIS第一件事就是让 AI 先跑一遍数据质量检查再开始分析。多花 1 分钟检查少花 1 小时排查。写在最后AI 确实能大幅提升效率但在刚开始用于分析空间数据时它还不是 GIS 专家需要多轮对话确定他的工作程度从简单任务到复杂任务训练其逐渐成为专家就像咱们刚开始接触空间地理信息数据一样逐步形成空间分析的逻辑闭环。它不知道“1 度在赤道和北极不一样长”不知道“面积计算必须用投影坐标”不知道“Shapefile 字段名会截断”——这些是你才知道的事。AI 是工具核心还是人。把你的领域知识教给它它才能真的帮上忙。互动你用 AI 做 GIS 时踩过什么坑评论区说说看看谁的最离谱。【GIS与AI】深度了解GISAI 实战经验、源码模板、避坑指南