ArcPy实战：坐标系统一、渔网创建与空间选择全流程

前阵子接了一批数据处理原始数据来源五花八门——有从公开平台下载的有甲方直接甩过来的CAD转换成果还有从旧系统里导出的历史库。量不算大但坐标系乱得让人头疼有的是地理坐标GCS有的是投影坐标PCS还有的根本没有任何空间参考信息。最麻烦的是这批数据一会儿要跟影像叠加一会儿要跟行政区边界做匹配如果不在前期把坐标系理顺后面每一步都可能出莫名其妙的错位。后来我用ArcPy把整套流程串起来了先统一坐标系再建渔网作为分析格网然后按位置关系做选择与导出最后回填字段属性。四个动作一气呵成半小时搞定了原本要手动折腾一下午的活儿。这篇就围绕这四个环节把我实际跑下来的思路、参数细节和踩过的坑都写出来给正在被同类问题折磨的朋友一个现成的参考。1. 坐标系读取先搞清楚数据“户口”再动手1.1 为什么第一步必须是读坐标系很多人在处理矢量数据时会下意识跳过坐标系检查尤其是数据本身看起来“能对上”的时候。但坐标系不是玄学它是数据坐标值的“解释器”同样的X、Y数值在WGS84下表示的是经纬度在投影坐标系下则是米。如果你不清楚一个数据的坐标系后续做缓冲、算面积、叠加分析都有可能产生量级上的偏差。真正让我意识到这步重要的是一次实战手头有份某区域的规划边界甲方说是“2000国家大地坐标系”但加载后发现和另一份标准底图偏离了大约几百米。逐项排查后发现这份边界数据虽然属性里写了CGCS2000但实际坐标值却像是基于旧基准采集的必须重新定义坐标系才能纠正偏移。这种情况如果不先读取、再核对仅靠肉眼判断几乎不可能发现。1.2 Describe读取空间参考代码与解读ArcPy里获取矢量文件坐标系非常直接核心就是Describe对象。几乎所有GIS数据源都支持这个操作包括shapefile、File Geodatabase要素类、CAD要素等只是有些源的返回信息会稍有差别。import arcpy fc rD:\work\data\boundary.shp desc arcpy.Describe(fc) sr desc.spatialReference print(名称: , sr.name) print(类型: , sr.type) print(工厂代码: , sr.factoryCode) print(投影坐标系: , sr.PCSName) print(地理坐标系: , sr.GCSName) print(线性单位: , sr.linearUnitName) print(WKID: , sr.factoryCode)这段代码输出多少信息取决于这个文件的空间参考是否完整。如果数据本身没有定义坐标系sr.name一般会返回Unknownsr.type也不是Geographic或Projected而是Uninitialized。这种情况并不代表坐标值没意义只是软件不知道该用什么规则去解释它们。1.3 从.prj文本看懂坐标系构成对于shapefile来说坐标系信息存放在同名.prj文件中读出来是一串WKT文本。很多人看到这串文本就头疼字母全认识、放一起就不认识了。我简单解释下常见关键字WKT关键字含义实际作用GEOGCS地理坐标系基于经纬度表示的坐标参照DATUM大地基准面决定椭球体模型例如WGS84、CGCS2000SPHEROID椭球体定义地球半径与扁率PRIMEM本初子午线一般指格林尼治UNIT单位角度度或距离米PROJECTION / PROJCS投影坐标系在GCS基础上做数学投影得到平面坐标你不需要背这串文本但至少要能分辨文件是GEOGCS开头还是PROJCS开头。前者代表地理坐标后者代表投影坐标。很多初学者把投影坐标里带的一串大数值X、Y当成经纬度去处理就是因为没分清这两个概念。1.4 数据源不同读取结果也会不同读取坐标系时还有几个容易忽略的差别。一是同一份数据在ArcGIS Pro和ArcMap中显示的坐标系名称有时不一样但这只是显示名称的差异WKID工厂代码才是真正意义上的“身份证”。二是从CAD软件转出来的数据经常没有.prj文件Describe返回的空间参考是空的。三是个人地理数据库与文件夹型shapefile在坐标系的存储机制上略有差异但通过Describe读取结果基本一致不需要特殊处理。所以我的建议是拿到任何一批矢量数据第一步先写个循环把每个文件的坐标系打印出来归档。这不是浪费时间而是在给后续所有操作“上保险”。2. 坐标系更改定义投影与投影转换完全是两件事2.1 先分清“定义”和“转换”坐标系更改在ArcPy里有两个完全不同的工具新人经常混淆Define Projection定义投影只是给没有坐标系的数据“补户口”。它的本质是修改数据自身的空间参考描述不改变坐标数值。例如一份数据坐标值本身就是CGCS2000的经纬度只是文件缺失了坐标系统信息用它去定义即可。Project投影转换真正改变坐标数值。比如把WGS84地理坐标转成Web Mercator投影坐标数值会从“度”变为“米”。用生活的方式理解定义投影是给你的地址贴上正确的省市区标签地址内容不动投影转换是把“北京市朝阳区”翻译成“北纬39.9度、东经116.4度”表达方式变了。2.2 定义投影的经典应用场景实际项目中定义投影最常见的场景是CAD数据导入GIS后没有任何坐标系信息但设计方口头告诉你采用的国家坐标标准。这时候你不能直接做投影转换因为软件不知道源坐标系是什么必须先告诉它“我的坐标是基于哪个空间参考下的”。import arcpy cad_fc rD:\work\data\land.dwg\Polygon sr arcpy.SpatialReference(4528) # 以CGCS2000 3度分带为例的WKID arcpy.DefineProjection_management(cad_fc, sr) print(坐标系定义完成, sr.name)这里有个细节DefineProjection会直接修改源数据的空间参考信息操作前务必备份。即使你只是想“试一试”一旦写错基准面数据本身的坐标不会被改写但下次打开时它会按错误的规则去解释坐标值往往会带来灾难性的后果。2.3 投影转换的参数与注意点相比之下Project会在原数据基础上生成一份新数据更安全。它主要包含四个关键参数输入要素、输出要素、输出坐标系、以及可选的地理变换。import arcpy in_fc rD:\work\data\boundary_wgs84.shp out_fc rD:\work\data\boundary_cgcs2000.shp # 如果源是WGS84目标是CGCS2000需要指定地理变换 out_sr arcpy.SpatialReference(4528) # 某些数据源之间的基准面差异不需要变换即可转换但地下数据转换通常需要 # 如果提示地理变换缺失需要查找并使用相应变换参数 arcpy.Project_management(in_fc, out_fc, out_sr) print(投影转换完成)当然这儿必须补充一点如果源地理坐标系与目标地理坐标系相同系统会自动跳过基准面变换只有跨基准面时才需要额外传入地理变换参数如WGS84与CGCS2000之间通常需要使用七参数或三参数模型。根据我实际项目的经验判断基准面是否一致最可靠的方式是查看两个坐标系各自的GCSName是否相同而不要只看它们的名称长得像不像。2.4 更改坐标系后必须做的核对坐标系一旦更改不能只看数据能加载就完事要核实转换结果的合理性。我个人的核对办法是用Describe读取输出数据的坐标系名称确认factoryCode符合预期再手动抽查几个关键点的坐标值判断单位量级是否符合预期比如是不是从度变成了米有条件的话叠加一份已知正确底图做目视比对。desc arcpy.Describe(out_fc) sr_out desc.spatialReference print(输出坐标系: , sr_out.name) print(WKID: , sr_out.factoryCode)坐标系统一之后再往后走渔网创建、空间选择、属性回填这些步骤就踏实多了。3. 渔网创建让分析格网跟着数据范围走3.1 渔网的本质给复杂分布建立统一参考框在做空间统计、抽样调查、统一管理这类任务时经常需要把目标区域切分成规则格网每个格网作为一个独立单元参与后续计算或查询。ArcPy中的CreateFishnet就是干这个的生成一个由矩形或正方形组成的渔网要素类。渔网创建逻辑不复杂但参数多、排列方向容易搞混。在开始写代码之前先明确一件事渔网的方向与原点、行列宽高密切相关如果只指定范围而没有给定合适的原点生成的网格可能与你的期望位置偏出半个格宽。3.2 CreateFishnet的核心参数拆解工具的参数大致可以分为三组范围、格网尺寸、输出类型。范围通常用一个Extent对象指定也可以用另一个要素类的范围来确定。渔网左上角默认从范围的左上角开始铺但如果尺寸不能被范围完整整除则边缘会有个不完整的格子。格网尺寸cell_width和cell_height分别指定单个格子的宽度和高度单位由渔网的坐标系决定。这两个参数也可以不写改为指定行列数二选一即可。几何类型geometry_type可以是POLYLINE或POLYGON一般分析用面状画格网线用线状。输出坐标系这个参数容易被忽略。如果不显式指定坐标系渔网默认使用arcpy.env.outputCoordinateSystem未设置则按目标区域环境或输入范围所在坐标系处理。import arcpy arcpy.env.workspace rD:\work\data arcpy.env.overwriteOutput True arcpy.env.outputCoordinateSystem arcpy.SpatialReference(4528) # 使用某个要素范围作为渔网的边界 template_fc rD:\work\data\study_area.shp desc arcpy.Describe(template_fc) extent desc.extent out_fc rD:\work\data\fishnet.shp arcpy.CreateFishnet_management( out_fcout_fc, origin_coordf{extent.XMin} {extent.YMin}, y_axis_coordf{extent.XMin} {extent.YMin 10}, cell_width1000, cell_height1000, number_rowsNone, number_columnsNone, corner_coordf{extent.XMax} {extent.YMax}, labelsNO_LABELS, templateextent, geometry_typePOLYGON ) print(渔网创建完成)y_axis_coord这个参数决定渔网的“垂直方向”它不直接影响生成结果但会影响旋转方向正常情况只要给它一个原点正北方向上方的坐标即可。template参数用来传递范围信息也可以是另一条数据。3.3 指定行数/列数的模式有时候你不知道单个格子的尺寸定多少合适但希望整个区域恰好被N×M等分。这时就用到另一种模式不填cell_width与cell_height而是填入number_rows与number_columns。arcpy.CreateFishnet_management( out_fcrD:\work\data\fishnet.shp, origin_coordf{extent.XMin} {extent.YMin}, y_axis_coordf{extent.XMin} {extent.YMin 10}, cell_widthNone, cell_heightNone, number_rows10, number_columns20, corner_coordf{extent.XMax} {extent.YMax}, labelsNO_LABELS, templateextent, geometry_typePOLYGON )需要注意的是这个模式下实际生成的格网尺寸是按照范围与行列数自动均分的每个格子都是完整的矩形不会出现半个格的情况相较指定尺寸更容易保证边界贴合。3.4 渔网创建的三个常见坑第一渔网范围如果不与目标数据保持一致后期空间关联时会有一堆边缘网格是空的。这里建议直接用要素类的extent对象获取最小外接矩形这样最稳。第二输出坐标系要在创建前设好否则可能出现“范围是投影坐标、但输出坐标系是地理坐标”这种矛盾设置工具虽然不报错但生成的格网变形得厉害。第三渔网创建的origin_coord如果随意指定哪怕只与范围左上角差半个格宽也会让全部格网错位。尤其当你后续要用格网的中心点或编号做匹配时这种错位几乎无法排查。4. 按位置选择导出空间关系的选择与批量处理4.1 为什么用“选择”而不是硬编码查坐标接下来这步很多新手喜欢直接写SQL条件按字段筛选但现实中的需求往往不是基于字段而是“看这张网点落在哪个格网里”或者“这条道路穿过哪些行政区”。这种带空间关系的需求用SQL写不出来需要用空间查询。ArcPy在这方面对应的核心工具是SelectLayerByLocation。它基于两个数据集之间的空间关系进行选择常见的关系词包括INTERSECT相交、WITHIN完全在内部、CONTAINS包含、HAVE_THEIR_CENTER_IN质心在内等。import arcpy arcpy.env.workspace rD:\work\data arcpy.env.overwriteOutput True target rD:\work\data\points.shp # 待筛选的要素 overlay rD:\work\data\fishnet.shp # 用于筛选的要素 # 先用目标数据创建要素图层 arcpy.MakeFeatureLayer_management(target, points_layer) # 按位置选择点要素与哪个网格相交 arcpy.SelectLayerByLocation_management( points_layer, overlap_typeINTERSECT, select_featuresoverlay, selection_typeNEW_SELECTION ) # 导出选中要素 out_fc rD:\work\data\points_selected.shp arcpy.CopyFeatures_management(points_layer, out_fc) print(按位置选择导出完成, arcpy.GetCount_management(out_fc))这一段代码就是标题中“按位置选择导出”最直接的实现。需要注意的是SelectLayerByLocation操作对象必须是一个要素图层Feature Layer而不是直接对要素类操作所以前面要通过MakeFeatureLayer先把要素类转为图层。4.2 不同关系词的选择逻辑实际项目中选择逻辑直接决定结果对不对。我整理了一下常用关系词适合的场景关系词判断标准适用场景INTERSECT几何有交集点线面混合判断、多对多匹配WITHIN目标完全位于另一要素内部点在多边形内、线在区域内部CONTAINS目标要素完全包含另一要素多边形包含点/线/面HAVE_THEIR_CENTER_IN质心落在另一要素内判断一个面的“重心位置”归属BOUNDARY_TOUCHES仅边界接触邻接关系分析例如你要筛选每个网格内面积占比最大的地块单靠“质心是否落入网格”来计算很可能会漏掉边界处的大地块。这种情况下更适合把目标要素“与网格相交”的格子全部选出来再进行后续的面积计算。4.3 批量处理时要时刻记住“清理选择”SelectLayerByLocation默认的selection_type是NEW_SELECTION这个参数意味着每次运行会刷新上一次的选择。但如果你在循环里反复使用同一个图层并且在某些条件下没有执行新的选择上一轮的选择结果会残留下来。我有一段踩过的坑循环40个格网为每个格网选择落在其中的要素并统计数量结果第15个格网开始数值异常偏大。排查后发现是某一次选择的要素较少程序提前退出了但图层里的选择集没清空下一次继续使用该图层时统计的还是旧选择结果。务必在每个循环的开头加上一行arcpy.SelectLayerByLocation_management(points_layer, CLEAR)在循环结构里这一步不仅能保证流程可复现还能避免脏数据污染统计结果。4.4 大批量导出时的替代方案如果只是偶尔选一次再导出SelectLayerByLocation CopyFeatures足够。但如果要按几十个网格分别导出独立文件这种方案会反复创建图层性能很差。这种情况我更建议直接用搜索游标按空间关系筛选后写入新文件或者一次性导出所有要素后在字段中标记网格编号再按字段批量拆分。空间查询没法直接做到“每一个目标要素检索它的归属网格编号”所以我的实际策略通常是先按位置选择再给选中的要素写入网格编号最后用Split by Attributes按编号拆分。# 按位置选择并写入网格ID with arcpy.da.SearchCursor(fishnet, [FID], spatial_filterarcpy.SpatialFilter(points, INTERSECT)) as cur: ...使用SearchCursor时可以传入spatial_filter参数这也是一个不依赖图层的空间查询方式。和图层选择相比游标方案内存占用低、适合批量循环但代码逻辑稍复杂两种方式各有适用场景。5. 字段赋值从计算字段到游标映射5.1 三种常见的赋值方式位置选择完成之后往往还要对成果做一个属性回填的小尾巴。比如给网格内的每个点写入网格编号或者根据位置关系打上分区标签。字段赋值在ArcPy中主要有三种方式CalculateField计算字段适合用表达式直接计算单个字段例如根据已有字段值生成新值或者做简单的条件判断。UpdateCursor更新游标适合批量逐行读取、判断并写入灵活性最好。字典映射本质上还是游标但提前把对应关系查成字典在循环里查字典赋值效率较高。5.2 使用UpdateCursor完成字段回填一个很典型的实际需求渔网创建后要把各网格内所含点要素的数量统计并回填到渔网属性表里。如果用SelectLayerByLocation逐个网格选择统计再去更新渔网的对应行逻辑上很绕。我更倾向于用UpdateCursor配合空间关系去实现。import arcpy arcpy.env.workspace rD:\work\data fishnet rD:\work\data\fishnet.shp points rD:\work\data\points.shp # 给渔网添加一个统计字段 arcpy.AddField_management(fishnet, PT_COUNT, LONG) # 空间查询的写法用SearchCursor查出每个点落在哪个网格 # 这里把查询结果先记在字典里再更新到渔网属性表 from collections import defaultdict count_map defaultdict(int) with arcpy.da.SearchCursor(points, [SHAPE, FID]) as cur: for row in cur: point_geom row[0] # 用一个简单的空间查询判断该点落在哪个渔网内部 # 最简单的办法是使用SelectLayerByLocation但如果数据量大建议直接做join pass如果要直接按空间关系做计数统计ArcPy没有直接提供“空间表关联”的工具接口但可以通过SpatialJoin工具完成把渔网作为目标要素点要素作为连接要素选择JOIN_ONE_TO_ONE统计每个网格内的点数。这一步得到的连接结果直接就是带统计字段的网格比手动循环高效得多。arcpy.SpatialJoin_analysis( target_featuresfishnet, join_featurespoints, out_feature_classrD:\work\data\fishnet_count.shp, join_operationJOIN_ONE_TO_ONE, join_typeKEEP_ALL, match_optionINTERSECT )说到字段赋值CalculateField在普通场景下也很顺手。比如给网格按行列号生成编码名称可以直接用Python表达式拼接arcpy.CalculateField_management( fishnet, GRID_CODE, G str(int(!LEFT!/1000)) _ str(int(!TOP!/1000)), PYTHON3 )这里的!LEFT!、!TOP!是要素类内置的几何属性分别代表要素外包络的左上角坐标算网格行列号非常方便省去了再建两个字段的步骤。5.3 游标赋值时的“顺手”与“陷阱”使用UpdateCursor时光标默认会持有数据集上的一份锁。如果赋值结束后没有及时释放游标后续对同图层执行重命名或删除操作就会报错。正确的写法是用上下文管理器with arcpy.da.UpdateCursor(fishnet, [GRID_ID, AREA_HA]) as cursor: for row in cursor: grid_id row[0] # 原值 # 根据原值计算出新值再写回去 row[1] round(float(row[1]) / 10000, 4) cursor.updateRow(row)这里有个不算违反直觉、但很容易出错的点cursor.updateRow更新的是当前row的引用如果你在循环里修改了row的某些元素必须把更新后的row整体传回去。直观上很多新手会误以为修改row后只要继续循环数据就会自动保存结果跑完了字段值没有变化。5.4 字段名与表达式的兼容性中文表头、特殊字符字段名在ArcPy里容易踩坑。有些字段名是中文或者以数字开头直接写表达式时会报错或返回空值。稳妥的办法是先用ListFields确认字段名再在表达式里用!字段名!括起来。fields [f.name for f in arcpy.ListFields(fishnet)] if TYPE_NAME not in fields: arcpy.AddField_management(fishnet, TYPE_NAME, TEXT)拼接字段、中文赋值也是一样CalculateField的表达式里直接写中文字符串很容易因编码问题引发错误可以先赋值给一个变量再用格式化字符串写入。6. 一个完整的处理链路读取、统一、建网、选择、回填6.1 串起来看整条流程单独看每个步骤都不难但真正的问题是怎么把它们组合成一条高效的流水线。这里给出一个完整的示例脚本骨架覆盖了标题里的所有环节读取坐标系、统一坐标系、按范围创建渔网、按位置选择导出、字段赋值回填。import arcpy from collections import defaultdict arcpy.env.workspace rD:\work\data arcpy.env.overwriteOutput True # 1. 读取源数据坐标系 source rD:\work\data\source.shp desc arcpy.Describe(source) print(源坐标系, desc.spatialReference.name) # 2. 统一定义坐标系 sr arcpy.SpatialReference(4528) if desc.spatialReference.factoryCode ! sr.factoryCode: out_source rD:\work\data\source_cgcs2000.shp arcpy.Project_management(source, out_source, sr) else: out_source source # 3. 使用要素范围创建渔网 desc2 arcpy.Describe(out_source) extent desc2.extent fishnet rD:\work\data\fishnet.shp arcpy.CreateFishnet_management( out_fcfishnet, origin_coordf{extent.XMin} {extent.YMin}, y_axis_coordf{extent.XMin} {extent.YMin 10}, cell_width500, cell_height500, number_rowsNone, number_columnsNone, corner_coordf{extent.XMax} {extent.YMax}, labelsNO_LABELS, templateextent, geometry_typePOLYGON ) # 4. 按位置选择并导出 arcpy.MakeFeatureLayer_management(out_source, source_layer) # 选择并导出所有与渔网相交的点位 arcpy.SelectLayerByLocation_management( source_layer, INTERSECT, fishnet, selection_typeNEW_SELECTION ) selected_fc rD:\work\data\selected.shp arcpy.CopyFeatures_management(source_layer, selected_fc) # 5. 字段赋值给渔网添加网格编号 arcpy.AddField_management(fishnet, GRID_CODE, TEXT) arcpy.CalculateField_management( fishnet, GRID_CODE, G str(int(!LEFT!/500)) _ str(int(!TOP!/500)), PYTHON3 ) print(全流程完成)这段脚本只是一个最基础的骨架。实际项目里我还会加上数据格式检查、坐标系判断分支、异常捕获和日志记录确保中途失败时能快速定位出错环节。6.2 脚本化的执行效率与内存管理在处理几百条要素的小数据时上述脚本几乎瞬间跑完。但如果面对的是几十万条数据有两个地方性能差异很大一是SelectLayerByLocation配合CopyFeatures在大数据量下会明显变慢尤其是被选中的数据量很庞大时CopyFeatures的写盘压力很大。这时可以尝试用arcpy.Select_analysis在图层上进行选择输出或者干脆用属性查询代替空间查询。二是CreateFishnet生成百万级网格时建议改为先按范围切分再合并而不是一次生成全部。否则后续的字段赋值和空间连接都会变得卡顿。6.3 关于运行环境与许可的小提示ArcPy相关脚本只能运行在安装了ArcGIS Desktop或Pro的机器上并且要有对应等级的许可。如果平时习惯使用开源GIS栈这个脚本对你可能不具备直接运行条件但“坐标系统一→规则网格→空间选择→属性回填”这套思路仍然值得借鉴。具体到ArcPy这一侧几个小细节有助于减少失败率路径中不要带有中文或空格尤其不要放在桌面或临时目录下防止字符编码引发隐性问题。每次运行前设置arcpy.env.overwriteOutput True否则输出已存在时脚本会中断报错。脚本执行前先用ArcGIS Pro或ArcMap手动加载一次数据确认数据本身能正常打开。如果连界面里都加载不了再查脚本意义也不大。6.4 一个并非偶然的用法把网格编号写回源要素做完整流程时我经常会多走一步把渔网的网格编号再空间关联回源要素这样就不需要反复做按位置选择导出以后只要用属性筛选就能随时取得某个格网内的数据。arcpy.SpatialJoin_analysis( target_featuresarcpy.MakeFeatureLayer_management(out_source, src_lyr), join_featuresfishnet, out_feature_classrD:\work\data\source_with_grid.shp, join_operationJOIN_ONE_TO_ONE, join_typeKEEP_ALL, match_optionHAVE_THEIR_CENTER_IN )这样源要素与格网编号建立起了长期可用的关联后续报表统计、专题图制作都可以直接基于这个结果不必再重复空间查询。这套流程跑下来我对“批量数据处理”这件事的体会是坐标系、网格、空间关系这几件事看似是独立的小工具组合起来却几乎是矢量数据预处理的标准底座。别嫌麻烦前期花十分钟把坐标系理清楚建网时多确认一下范围和坐标系的匹配后面所有的空间选择、属性赋值都会顺很多。最怕的反而是每个环节单独看起来都“差不多”结果最后叠在一起偏差大到无从下手。