尧图网页设计 
这个标题我一看就挺有共鸣——Acgis中实现栅格经纬度和行政区关联这里的Acgis我理解就是ArcGIS桌面或ArcGIS Pro这套GIS平台。这类需求在项目实施里特别常见手里有一张栅格数据遥感影像、DEM高程、降雨量插值、土地分类结果要弄清楚每个像元的经纬度落在哪个省、市、县或者反过来要把行政区边界内的栅格值统计出来。别看问题描述简单实际处理时坐标系、字段类型、数据量、边界归属任何一个环节不留意结果就会对不上。这篇就按我实际做过的思路把几种靠谱的做法和踩过的坑一次讲清楚。适合正在做空间数据分析、GIS数据加工、遥感应用的同学参考。1. 先搞清楚需求你要的关联到底是哪一种1.1 三种常见需求拆解拿到栅格经纬度和行政区关联这个需求我第一件事不是打开ArcGIS而是问清楚对方到底要什么结果。因为这个说法至少对应三种完全不同的技术路径。第一种最常见给栅格的每个像元打上行政区标签。比如你有一张土地利用分类栅格每个像元是一个地类代码现在领导要求每个像元旁边标记出它属于哪个县后续好按县统计地类面积。这种需求的核心是空间位置判读需要把栅格像元转成可参与空间分析的对象再和行政区面做叠加。第二种需求是把点状数据或者栅格像元中心点的经纬度坐标直接换成行政区名称或区划代码。经常出现在外业采集场景手里有一堆GPS点或者采样点坐标是Excel里的两列经纬度要快速判断每个点落在哪个行政村。这种情况下栅格反而是辅助核心是坐标点落面判断。第三种是按行政区范围做栅格值统计。比如你有一张全国降雨量栅格要算每个县的平均降雨量、最大降雨量、累计降雨量。这种需求根本不关心每个像元具体属于哪个县只需要以县界为范围对栅格做聚合统计对应ArcGIS里的分区统计工具。这三种需求如果不先分清楚很容易做着做着就发现方案选错了。我见过有同事用第一种的转点空间连接方法去做第三种数据量一上去直接卡死最后换成分区统计几分钟就出结果。1.2 为什么不能用普通的表关联来解决很多刚接触GIS的朋友会把关联理解成数据库里的Join两张表有个共同字段按照字段值把属性拼在一起。但经纬度和行政区之间没有这个共同字段你不能说某个经纬度等于某个区划代码它们是两种完全不同的信息。从底层逻辑讲经纬度是坐标值描述的是地球上某个点的位置行政区是面要素描述的是一个由边界围起来的区域。要把两者关联起来必须通过空间关系判断——也就是判断这个点是否落在某个面内部。这个判断在GIS里叫Point-in-Polygon点在多边形内判断ArcGIS里对应的工具是空间连接Spatial Join和按位置选择Select By Location。还有一个容易被忽视的点一个Excel表格里写着经纬度的两列数据在ArcGIS里还不是要素它只是普通表格。你至少要用添加XY数据或者创建要素类功能把它变成真正的点图层才能参与空间关联。这就像你手里有一张写着门牌号的纸条但你没有走到门前系统不知道这个门牌号对应的是哪栋房子。2. 数据准备坐标系不统一后面全是白干2.1 经纬度到底算什么坐标系栅格又算什么坐标系先说一个基础但特别容易栽跟头的知识点。经纬度本身不是坐标系它是坐标值的一种表示方式单位是度。在ArcGIS里经纬度坐标默认对应的是地理坐标系比如WGS84、CGCS2000、西安80、北京54这些属于同一个大类只是参考椭球和基准面不同。栅格数据则可能是地理坐标系也可能是投影坐标系。如果栅格是地理坐标系那像元中心坐标直接就是经纬度如果栅格是投影坐标系像元中心坐标就是投影平面上的坐标单位是米比如Web Mercator下的X、Y值或者是高斯投影下的公里网值。行政区边界数据通常是矢量面坐标系五花八门常见的有CGCS2000、西安80、北京54也有GPS外业直接拿回来的WGS84。做关联之前如果坐标系不统一哪怕所有操作步骤都对结果也是偏移的。这就好比你把一张旧地图和一个GPS定位叠在一起地图上明明显示在河南实际人已经走到河北了。2.2 定义投影和投影转换是两个概念别弄混这个坑我在项目里见过不止一次。有个同事拿到一个没有坐标系信息的DEM栅格他以为用定义投影工具选了CGCS2000就能让它坐标正确。但实际上定义投影只是给数据打一个标签告诉ArcGIS这个数据的坐标是CGCS2000它不会改变数据里真实存储的坐标值。如果这个数据本身就是CGCS2000坐标只是缺少坐标系信息那么定义投影是对的如果这个数据本来是WGS84经纬度只是没写坐标系你直接定义为CGCS2000它还是经纬度的数值但被ArcGIS误认为是CGCS2000的经纬度后续无论是显示还是转换都会错。正确的做法是先搞清楚数据来源。如果是国家基础地理信息中心下载的数据大概率是CGCS2000如果是GPS设备采集大概率是WGS84如果是老项目数据可能是西安80甚至北京54。确认之后再决定用投影工具做转换还是用定义投影补标签。转换的时候要注意WGS84转CGCS2000在多数地区差异不大但西安80转CGCS2000或者北京54转CGCS2000通常需要七参数或四参数不能直接忽略。2.3 数据格式与字段类型也要提前核对栅格数据本身有两个特性会直接影响关联方案。第一多波段栅格不能直接用栅格转点工具工具只支持单波段栅格。如果你手里是一张多波段的遥感影像需要先用提取单波段或者波段合成处理一下取其中一个波段参与关联。第二栅格的NoData像元在转点时不会被输出这会导致转出的点数量比栅格行列数乘积少如果后续要按点数量核对结果要提前想到这一点。行政区数据这边最关键的字段是行政区代码。中国行政区划代码是6位数字比如北京东城区是110101。这个字段在Shapefile里通常是字符串类型但如果你用面转栅格把行政区面转成栅格时选了它作为Value字段ArcGIS可能会把字符串转成数值型这样110101还是110101但有些以0开头的代码就会丢掉前导0后面关联的时候就对不上了。稳妥的做法是在转换前先把行政区代码复制成一个整型字段或者干脆用已有的唯一ID字段作为Value之后再用ID关联回行政区属性表。3. 实操方案一栅格转点 空间连接中小数据量首选3.1 整体思路一看就懂这个方案的思路特别直白栅格不是由一个个像元组成的吗我先把每个像元转成一个点要素让这个点带上像元的中心坐标和像元值然后再用空间连接把这些点和行政区面做叠加行政区的属性自然就落到了每个点上。最后得到的是一个带行政区名称、行政区代码、栅格值的点图层。这个方案最大的优点是直观、好排查。任何一步结果都可以打开属性表检查比如某个点为什么没有关联上行政区是落在了边界缝隙里还是坐标系偏了一眼就能看出来。缺点是数据量大了以后性能差像元一多转点动辄几百万上千万ArcGIS处理起来非常吃力。我一般建议像元数量在几十万到一两百万以内用这个方法超过这个量级直接看后面的脚本方案。3.2 栅格转点工具的参数细节在ArcToolbox里栅格转点位于 Conversion Tools转换工具下面的 From Raster从栅格转换里ArcGIS Pro里直接在地理处理窗格搜索中文名也行。工具参数并不多输入栅格、字段、输出点要素、是否只转内部像元。字段那里一般选Value也就是把像元值写到点的属性表里。如果你栅格的值本身就是分类代码这个字段后续会很有用如果只是高程或温度等连续值也可以选反正以后可以删。只转内部像元这个选项默认是勾选的意思是栅格边缘不完整像元不参与转换。如果你的栅格范围本身不是整数像元边界转出来的点会少一圈这时候要注意核对。还有一个重要限制多波段栅格必须提前拆成单波段。实际操作中如果输入是遥感影像可以先复制一个波段或者用提取单波段工具把需要的波段单独导出成一个单波段栅格再做栅格转点。3.3 空间连接与结果检查拿到栅格转出来的点要素后下一步用空间连接Spatial Join位置在 Analysis Tools分析工具下的 Overlay叠加分析里。目标要素选点连接要素选行政区面输出要素类自己命名连接操作选JOIN_ONE_TO_ONE匹配选项选INTERSECT。这里有一个容易踩的坑如果点恰好落在行政区面的边界上ArcGIS可能匹配不到或者匹配出多个面。默认情况下JOIN_ONE_TO_ONE只保留第一个匹配结果但到底是哪一个它的规则并不绝对可控。如果你的数据里这种情况很多可以把匹配选项改成CLOSEST这样至少能保底匹配到最近的一个行政区。匹配完之后一定要做检查。打开结果点图层的属性表看行政区代码字段有没有空值。有空的用按属性选择把空值记录选出来看看它们的坐标位置多半是落在行政区的缝隙里或者数据范围外。如果空值比例很高优先检查坐标系是否统一不要急着手动补数据。整个流程结束后如果最终交付需要的是一个带行政区属性的栅格而不是点图层可以用点转栅格再做一次栅格化字段选行政区代码像元大小设成和原始栅格一致。不过我个人建议除非必须尽量保留点图层交付因为栅格化之后属性表就没那么灵活了而且点转栅格容易在边缘出现空白缝隙。4. 实操方案二分区统计与面积制表按行政区汇总4.1 分区统计的适用场景只要统计值不需要逐像元标签如果需求是每个县的平均降雨量是多少每个乡镇的DEM最高海拔是多少每个地市的土地利用面积怎么分布那就完全没必要把每个像元都转成点去做空间连接。这时候该用分区统计Zonal Statistics。分区统计的底层逻辑是以行政区面为分区对落在每个分区内的栅格像元值做统计计算输出一张统计表。统计类型支持均值、最大值、最小值、总和、标准差、中值、众数、极差等。工具在ArcToolbox里是 Spatial Analyst Tools空间分析工具下的 Zonal分区工具集。有一个细节很多人不知道分区数据可以直接用行政区面要素不一定非要先转栅格。工具内部会自动把矢量面转成栅格参与运算。但这里有个潜在问题工具内部转栅格时像元大小怎么定你是控制不了的如果你直接用它默认设置统计结果可能和你预期栅格匹配得不好。所以我一般建议手动先把行政区面转成栅格设置好像元大小再去做分区统计全程自己掌握。4.2 面转栅格的像元大小和范围匹配手动转行政区面的时候工具是面转栅格Polygon to Raster参数包括输入要素、值字段、输出栅格、像元大小。值字段建议选行政区代码或一个唯一的ID。做这一步之前先在环境设置里把捕捉栅格Snap Raster设成你的原始数据栅格把像元大小也设成和原始栅格一致。这样做的好处是转出来的行政区栅格和原始栅格在像元网格上完全对齐后续分区统计每个像元都精确对应不会出现边缘像元错位。范围设置也很关键。如果行政区范围比栅格范围大统计时多出来的区域全是NoData如果栅格范围比行政区范围大统计结果会把行政区之外的像元也算进去。最稳妥的做法是在分区统计的环境设置里把处理范围设成两者的交集或者干脆设成原始栅格的范围再把行政区和栅格的NoData都处理好。4.3 Tabulate Area分类栅格统计面积的好帮手和分区统计配套的还有一个工具叫面积制表Tabulate Area专门用来做分类栅格的面积统计。比如你有一张土地利用分类栅格地类代码有建设用地、耕地、林地、草地等你想知道每个县各类别面积有多少这时候直接用Tabulate Area最合适。它不需要你指定统计类型工具会统计每个分区内每个类别值的像元数量然后乘以像元面积输出一张交叉表行是行政区列是地类值是面积。这个工具我在地类变化监测项目里用得非常多。做的时候同样要先把行政区面和分类栅格的像元网格对齐否则面积算出来会有偏差。有一点要注意分类栅格的值字段必须是整型如果原栅格是浮点型的分类代码先转成整型再做不然可能报错或者结果异常。5. 实操方案三Python脚本化处理数据量大、批处理5.1 为什么数据量大必须上脚本栅格如果是一个全省范围、30米分辨率的DEM像元数量轻松上亿。栅格转点空间连接在这种数据量下基本跑不动ArcMap很可能直接转圈一个下午最后还报内存不足。即便勉强跑完生成几个G的点要素打开都费劲。这时候就要换思路。我的做法是不去ArcGIS界面里点工具而是用Python脚本直接读取栅格的像元坐标结合行政区面的几何信息做空间判断。这样有几大好處第一可以控制内存分段读取第二可以跳过中间过程产生的几G临时数据只保留最终关联结果第三方便在其他项目里复用改一改路径和字段名就能跑。5.2 用arcpy NumPy读取栅格坐标并关联行政区先看一个常用的思路用arcpy把栅格读成NumPy数组然后根据栅格原点和像元大小推算出每个像元中心点的经纬度坐标再把这些坐标点交给空间判断逻辑处理。import arcpy import numpy as np # 设置工作环境 arcpy.env.workspace rC:\gis_data arcpy.env.overwriteOutput True # 读取栅格为numpy数组 raster_path rC:\gis_data\dem.tif raster arcpy.Raster(raster_path) arr arcpy.RasterToNumPyArray(raster, nodata_to_value-9999) # 获取栅格信息 desc arcpy.Describe(raster) origin_x desc.extent.XMin origin_y desc.extent.YMin cell_w desc.meanCellWidth cell_h desc.meanCellHeight rows, cols arr.shape # 计算每一个像元的中心点坐标 # 注意numpy行号向下对应Y减小列号向右对应X增大 x_coords origin_x (np.arange(cols) 0.5) * cell_w y_coords origin_y (rows - np.arange(rows) - 0.5) * cell_h X, Y np.meshgrid(x_coords, y_coords)上面的X和Y就对应每个像元中心点的平面坐标。如果栅格是地理坐标系比如WGS84或CGCS2000那这些坐标就是经纬度。接下来有两种做法。一种是把坐标点批量转成点要素再交给arcpy的空间连接工具处理。另一种是直接用arcpy的几何对象做判断。批量创建点要素在大数据量下也不是特别高效所以更常见的是配合空间索引或者只用必要的点做判断。5.3 用geopandas shapely做空间关联如果电脑上装了Python的geopandas库有另一套更简洁的路线把栅格的每个像元中心点构造成GeoDataFrame把行政区shp读成另一个GeoDataFrame然后直接用空间连接函数sjoin一次完成关联。这个方法不需要ArcGIS许可跑起来也比较快尤其适合批量处理多个栅格文件。import geopandas as gpd import pandas as pd import numpy as np from shapely.geometry import Point # 假设已经用前面类似的方式得到 x_flat, y_flat, value_flat points pd.DataFrame({ x: x_flat, y: y_flat, value: value_flat }) # 把经纬度坐标转成点要素 gdf_points gpd.GeoDataFrame( points, geometrygpd.points_from_xy(points.x, points.y), crsEPSG:4326 ) # 读取行政区面并统一坐标系 gdf_adm gpd.read_file(rC:\gis_data\county.shp) gdf_adm gdf_adm.to_crs(EPSG:4326) # 空间连接每个点匹配所在的行政区面 joined gpd.sjoin( gdf_points, gdf_adm, howleft, predicatewithin ) # 保存结果 joined.to_file(rC:\gis_data\points_with_admin.shp, encodingutf-8)这段代码里最需要留意的是内存。如果像元数量超过几百万一次性把全部点放进GeoDataFrame再sjoin内存会吃得非常狠机器不好的会直接崩。我实际处理时是按行数分块的比如每50万行读取一次栅格切片转成点关联完保存一批再处理下一批。这样做虽然代码稍长但稳定。5.4 两种脚本方案的选型建议我自己比较习惯的分界线是这样的如果栅格数量不大但行政区面很细碎用arcpy方案更顺手因为可以直接复用ArcGIS自带的空间关系和坐标系管理能力。如果栅格本身很大或者你有多个栅格要批量关联用geopandas方案更好因为它写起来简洁而且不依赖ArcGIS许可可以在自己的服务器或笔记本上跑。无论是哪种脚本方案开工前都建议先做一次小范围的测试裁剪一小块栅格跑通全流程确认坐标、字段、输出都对再放全量数据去跑。否则全量跑完才发现坐标系没对齐几小时就白费了。6. 常见问题与排查技巧实录6.1 关联结果发生偏移该怎么办如果关联出来的点落到相邻的行政区里或者反过来无法匹配第一件事检查坐标系是否完全统一。在ArcGIS里把点图层和行政区图层都打开属性列表里的源标签看坐标系到底是不是同一个。如果一个是WGS84经纬度一个是CGCS2000投影坐标它们在地图上可能显示位置差不多但精度完全不一样空间连接的结果就不可靠。建议全部统一到一个坐标系再做关联。多数情况下我统一到CGCS2000地理坐标系单位是度便于直接用经纬度描述结果。如果坐标系看起来是统一的但偏移仍然存在检查是否做了定义投影而不是投影转换。最简单的验证方法是在原始栅格的属性中看范围和坐标值假如文件是WGS84经纬度范围应该大约在73度到135度、18度到54度之间如果文件被错误定义为CGCS2000投影坐标范围可能会变成几百万米的巨大数值。这种情况只能先把坐标系定义改正再重新投影。6.2 栅格转点后数量对不上少了或多了一些点点数量比像元总数少最常见的原因是NoData。栅格里的NoData像元不会转成点如果你不确定哪些是NoData可以在栅格属性里看统计信息或者用栅格计算器把NoData像元单独提取出来看一下。如果NoData很多先做填洼、插值或者其他预处理再去转点。点数量比预期多可能是因为栅格被重采样过像元大小发生变化或者输入栅格是多波段误把多个波段都转了。使用工具时注意输入必须是一个单波段栅格并且输出的像元大小和原始栅格一致。6.3 行政区代码前导0丢失导致无法关联行政区代码如果是字符串面转栅格时如果字段被自动转成整型前导0就会丢。最典型的是某些地区代码以0开头比如有些乡镇级代码是010开头转成整型后变成了10后面想用这个代码关联回行政区属性表就完全对不上了。解决办法有两个一是在面转栅格之前用添加字段字段计算器生成一个整型的唯一ID面转栅格时用这个ID作为Value字段后面统计结果再用ID关联回行政区属性表取得完整代码。二是如果一定要保留字符串代码建议不要转栅格而是走栅格转点的方案因为点属性表里字符串字段不会丢前导0。6.4 数据量太大卡死如何优化如果栅格像元数量在百万级以下栅格转点和空间连接还能撑住一旦到千万级以上就要考虑三点优化。第一用环境设置里的处理范围先把分析区域裁剪到目标行政区内不要全图跑第二先用行政区面的外接矩形裁剪栅格缩小数据量第三用脚本方案按行分块处理避免一次占用太多内存。另外空间连接时尽量精简输出字段。默认情况下Spatial Join会把连接要素的所有属性都带过去如果行政区面有几十个字段输出点要素的表会非常臃肿拖慢后续操作。在工具参数里可以设置字段映射只保留行政区代码和名称其他全删掉运行速度会明显提升。6.5 边界上的点归属判断不一致落在行政区边界线上的点用INTERSECT匹配时结果不确定尤其两个面共用一个边界时。ArcGIS的JOIN_ONE_TO_ONE默认会选一个但具体选哪个不完全可控。为了结果稳定我习惯先用按位置选择或脚本把这种边界点筛出来再用最近距离归属到一个面。这里可以设置一个很小的空间容差比如1米让边界点优先匹配到面积更小或者行政中心所在的那一侧具体需要结合项目规则判断。6.6 常见问题速查表现象可能原因排查方向关联后大量空值坐标系不统一检查点图层和面图层坐标系关联结果偏移到相邻区域做错了投影对比栅格范围确认坐标值是否符合预期转点数量少于像元数存在NoData统计NoData数量确认是否可忽略行政区代码对不上前导0丢失检查面转栅格字段类型改用ID关联工具运行卡死数据量过大裁剪范围、精简字段、改脚本分块处理边界点归属不确定匹配规则问题改用CLOSEST或设置空间容差关于这类需求我个人的习惯是只要不是必须保留栅格成果我都优先走栅格转点空间连接这条路因为点数据随时可以核对、转换、导Excel出了问题也容易定位。如果最终非要栅格成果我一般先把行政区面转成栅格像元大小和范围都在环境设置里锁死再基于这个栅格做后续操作这样每个环节都在自己掌控中。坐标系方面新做的项目我一般统一到CGCS2000再开始处理。这个流程踩过不少坑之后总结出来的希望对做类似需求的同行有帮助。 ArcGIS中实现栅格经纬度与行政区关联的三种实用方法
这个标题我一看就挺有共鸣——Acgis中实现栅格经纬度和行政区关联这里的Acgis我理解就是ArcGIS桌面或ArcGIS Pro这套GIS平台。这类需求在项目实施里特别常见手里有一张栅格数据遥感影像、DEM高程、降雨量插值、土地分类结果要弄清楚每个像元的经纬度落在哪个省、市、县或者反过来要把行政区边界内的栅格值统计出来。别看问题描述简单实际处理时坐标系、字段类型、数据量、边界归属任何一个环节不留意结果就会对不上。这篇就按我实际做过的思路把几种靠谱的做法和踩过的坑一次讲清楚。适合正在做空间数据分析、GIS数据加工、遥感应用的同学参考。1. 先搞清楚需求你要的关联到底是哪一种1.1 三种常见需求拆解拿到栅格经纬度和行政区关联这个需求我第一件事不是打开ArcGIS而是问清楚对方到底要什么结果。因为这个说法至少对应三种完全不同的技术路径。第一种最常见给栅格的每个像元打上行政区标签。比如你有一张土地利用分类栅格每个像元是一个地类代码现在领导要求每个像元旁边标记出它属于哪个县后续好按县统计地类面积。这种需求的核心是空间位置判读需要把栅格像元转成可参与空间分析的对象再和行政区面做叠加。第二种需求是把点状数据或者栅格像元中心点的经纬度坐标直接换成行政区名称或区划代码。经常出现在外业采集场景手里有一堆GPS点或者采样点坐标是Excel里的两列经纬度要快速判断每个点落在哪个行政村。这种情况下栅格反而是辅助核心是坐标点落面判断。第三种是按行政区范围做栅格值统计。比如你有一张全国降雨量栅格要算每个县的平均降雨量、最大降雨量、累计降雨量。这种需求根本不关心每个像元具体属于哪个县只需要以县界为范围对栅格做聚合统计对应ArcGIS里的分区统计工具。这三种需求如果不先分清楚很容易做着做着就发现方案选错了。我见过有同事用第一种的转点空间连接方法去做第三种数据量一上去直接卡死最后换成分区统计几分钟就出结果。1.2 为什么不能用普通的表关联来解决很多刚接触GIS的朋友会把关联理解成数据库里的Join两张表有个共同字段按照字段值把属性拼在一起。但经纬度和行政区之间没有这个共同字段你不能说某个经纬度等于某个区划代码它们是两种完全不同的信息。从底层逻辑讲经纬度是坐标值描述的是地球上某个点的位置行政区是面要素描述的是一个由边界围起来的区域。要把两者关联起来必须通过空间关系判断——也就是判断这个点是否落在某个面内部。这个判断在GIS里叫Point-in-Polygon点在多边形内判断ArcGIS里对应的工具是空间连接Spatial Join和按位置选择Select By Location。还有一个容易被忽视的点一个Excel表格里写着经纬度的两列数据在ArcGIS里还不是要素它只是普通表格。你至少要用添加XY数据或者创建要素类功能把它变成真正的点图层才能参与空间关联。这就像你手里有一张写着门牌号的纸条但你没有走到门前系统不知道这个门牌号对应的是哪栋房子。2. 数据准备坐标系不统一后面全是白干2.1 经纬度到底算什么坐标系栅格又算什么坐标系先说一个基础但特别容易栽跟头的知识点。经纬度本身不是坐标系它是坐标值的一种表示方式单位是度。在ArcGIS里经纬度坐标默认对应的是地理坐标系比如WGS84、CGCS2000、西安80、北京54这些属于同一个大类只是参考椭球和基准面不同。栅格数据则可能是地理坐标系也可能是投影坐标系。如果栅格是地理坐标系那像元中心坐标直接就是经纬度如果栅格是投影坐标系像元中心坐标就是投影平面上的坐标单位是米比如Web Mercator下的X、Y值或者是高斯投影下的公里网值。行政区边界数据通常是矢量面坐标系五花八门常见的有CGCS2000、西安80、北京54也有GPS外业直接拿回来的WGS84。做关联之前如果坐标系不统一哪怕所有操作步骤都对结果也是偏移的。这就好比你把一张旧地图和一个GPS定位叠在一起地图上明明显示在河南实际人已经走到河北了。2.2 定义投影和投影转换是两个概念别弄混这个坑我在项目里见过不止一次。有个同事拿到一个没有坐标系信息的DEM栅格他以为用定义投影工具选了CGCS2000就能让它坐标正确。但实际上定义投影只是给数据打一个标签告诉ArcGIS这个数据的坐标是CGCS2000它不会改变数据里真实存储的坐标值。如果这个数据本身就是CGCS2000坐标只是缺少坐标系信息那么定义投影是对的如果这个数据本来是WGS84经纬度只是没写坐标系你直接定义为CGCS2000它还是经纬度的数值但被ArcGIS误认为是CGCS2000的经纬度后续无论是显示还是转换都会错。正确的做法是先搞清楚数据来源。如果是国家基础地理信息中心下载的数据大概率是CGCS2000如果是GPS设备采集大概率是WGS84如果是老项目数据可能是西安80甚至北京54。确认之后再决定用投影工具做转换还是用定义投影补标签。转换的时候要注意WGS84转CGCS2000在多数地区差异不大但西安80转CGCS2000或者北京54转CGCS2000通常需要七参数或四参数不能直接忽略。2.3 数据格式与字段类型也要提前核对栅格数据本身有两个特性会直接影响关联方案。第一多波段栅格不能直接用栅格转点工具工具只支持单波段栅格。如果你手里是一张多波段的遥感影像需要先用提取单波段或者波段合成处理一下取其中一个波段参与关联。第二栅格的NoData像元在转点时不会被输出这会导致转出的点数量比栅格行列数乘积少如果后续要按点数量核对结果要提前想到这一点。行政区数据这边最关键的字段是行政区代码。中国行政区划代码是6位数字比如北京东城区是110101。这个字段在Shapefile里通常是字符串类型但如果你用面转栅格把行政区面转成栅格时选了它作为Value字段ArcGIS可能会把字符串转成数值型这样110101还是110101但有些以0开头的代码就会丢掉前导0后面关联的时候就对不上了。稳妥的做法是在转换前先把行政区代码复制成一个整型字段或者干脆用已有的唯一ID字段作为Value之后再用ID关联回行政区属性表。3. 实操方案一栅格转点 空间连接中小数据量首选3.1 整体思路一看就懂这个方案的思路特别直白栅格不是由一个个像元组成的吗我先把每个像元转成一个点要素让这个点带上像元的中心坐标和像元值然后再用空间连接把这些点和行政区面做叠加行政区的属性自然就落到了每个点上。最后得到的是一个带行政区名称、行政区代码、栅格值的点图层。这个方案最大的优点是直观、好排查。任何一步结果都可以打开属性表检查比如某个点为什么没有关联上行政区是落在了边界缝隙里还是坐标系偏了一眼就能看出来。缺点是数据量大了以后性能差像元一多转点动辄几百万上千万ArcGIS处理起来非常吃力。我一般建议像元数量在几十万到一两百万以内用这个方法超过这个量级直接看后面的脚本方案。3.2 栅格转点工具的参数细节在ArcToolbox里栅格转点位于 Conversion Tools转换工具下面的 From Raster从栅格转换里ArcGIS Pro里直接在地理处理窗格搜索中文名也行。工具参数并不多输入栅格、字段、输出点要素、是否只转内部像元。字段那里一般选Value也就是把像元值写到点的属性表里。如果你栅格的值本身就是分类代码这个字段后续会很有用如果只是高程或温度等连续值也可以选反正以后可以删。只转内部像元这个选项默认是勾选的意思是栅格边缘不完整像元不参与转换。如果你的栅格范围本身不是整数像元边界转出来的点会少一圈这时候要注意核对。还有一个重要限制多波段栅格必须提前拆成单波段。实际操作中如果输入是遥感影像可以先复制一个波段或者用提取单波段工具把需要的波段单独导出成一个单波段栅格再做栅格转点。3.3 空间连接与结果检查拿到栅格转出来的点要素后下一步用空间连接Spatial Join位置在 Analysis Tools分析工具下的 Overlay叠加分析里。目标要素选点连接要素选行政区面输出要素类自己命名连接操作选JOIN_ONE_TO_ONE匹配选项选INTERSECT。这里有一个容易踩的坑如果点恰好落在行政区面的边界上ArcGIS可能匹配不到或者匹配出多个面。默认情况下JOIN_ONE_TO_ONE只保留第一个匹配结果但到底是哪一个它的规则并不绝对可控。如果你的数据里这种情况很多可以把匹配选项改成CLOSEST这样至少能保底匹配到最近的一个行政区。匹配完之后一定要做检查。打开结果点图层的属性表看行政区代码字段有没有空值。有空的用按属性选择把空值记录选出来看看它们的坐标位置多半是落在行政区的缝隙里或者数据范围外。如果空值比例很高优先检查坐标系是否统一不要急着手动补数据。整个流程结束后如果最终交付需要的是一个带行政区属性的栅格而不是点图层可以用点转栅格再做一次栅格化字段选行政区代码像元大小设成和原始栅格一致。不过我个人建议除非必须尽量保留点图层交付因为栅格化之后属性表就没那么灵活了而且点转栅格容易在边缘出现空白缝隙。4. 实操方案二分区统计与面积制表按行政区汇总4.1 分区统计的适用场景只要统计值不需要逐像元标签如果需求是每个县的平均降雨量是多少每个乡镇的DEM最高海拔是多少每个地市的土地利用面积怎么分布那就完全没必要把每个像元都转成点去做空间连接。这时候该用分区统计Zonal Statistics。分区统计的底层逻辑是以行政区面为分区对落在每个分区内的栅格像元值做统计计算输出一张统计表。统计类型支持均值、最大值、最小值、总和、标准差、中值、众数、极差等。工具在ArcToolbox里是 Spatial Analyst Tools空间分析工具下的 Zonal分区工具集。有一个细节很多人不知道分区数据可以直接用行政区面要素不一定非要先转栅格。工具内部会自动把矢量面转成栅格参与运算。但这里有个潜在问题工具内部转栅格时像元大小怎么定你是控制不了的如果你直接用它默认设置统计结果可能和你预期栅格匹配得不好。所以我一般建议手动先把行政区面转成栅格设置好像元大小再去做分区统计全程自己掌握。4.2 面转栅格的像元大小和范围匹配手动转行政区面的时候工具是面转栅格Polygon to Raster参数包括输入要素、值字段、输出栅格、像元大小。值字段建议选行政区代码或一个唯一的ID。做这一步之前先在环境设置里把捕捉栅格Snap Raster设成你的原始数据栅格把像元大小也设成和原始栅格一致。这样做的好处是转出来的行政区栅格和原始栅格在像元网格上完全对齐后续分区统计每个像元都精确对应不会出现边缘像元错位。范围设置也很关键。如果行政区范围比栅格范围大统计时多出来的区域全是NoData如果栅格范围比行政区范围大统计结果会把行政区之外的像元也算进去。最稳妥的做法是在分区统计的环境设置里把处理范围设成两者的交集或者干脆设成原始栅格的范围再把行政区和栅格的NoData都处理好。4.3 Tabulate Area分类栅格统计面积的好帮手和分区统计配套的还有一个工具叫面积制表Tabulate Area专门用来做分类栅格的面积统计。比如你有一张土地利用分类栅格地类代码有建设用地、耕地、林地、草地等你想知道每个县各类别面积有多少这时候直接用Tabulate Area最合适。它不需要你指定统计类型工具会统计每个分区内每个类别值的像元数量然后乘以像元面积输出一张交叉表行是行政区列是地类值是面积。这个工具我在地类变化监测项目里用得非常多。做的时候同样要先把行政区面和分类栅格的像元网格对齐否则面积算出来会有偏差。有一点要注意分类栅格的值字段必须是整型如果原栅格是浮点型的分类代码先转成整型再做不然可能报错或者结果异常。5. 实操方案三Python脚本化处理数据量大、批处理5.1 为什么数据量大必须上脚本栅格如果是一个全省范围、30米分辨率的DEM像元数量轻松上亿。栅格转点空间连接在这种数据量下基本跑不动ArcMap很可能直接转圈一个下午最后还报内存不足。即便勉强跑完生成几个G的点要素打开都费劲。这时候就要换思路。我的做法是不去ArcGIS界面里点工具而是用Python脚本直接读取栅格的像元坐标结合行政区面的几何信息做空间判断。这样有几大好處第一可以控制内存分段读取第二可以跳过中间过程产生的几G临时数据只保留最终关联结果第三方便在其他项目里复用改一改路径和字段名就能跑。5.2 用arcpy NumPy读取栅格坐标并关联行政区先看一个常用的思路用arcpy把栅格读成NumPy数组然后根据栅格原点和像元大小推算出每个像元中心点的经纬度坐标再把这些坐标点交给空间判断逻辑处理。import arcpy import numpy as np # 设置工作环境 arcpy.env.workspace rC:\gis_data arcpy.env.overwriteOutput True # 读取栅格为numpy数组 raster_path rC:\gis_data\dem.tif raster arcpy.Raster(raster_path) arr arcpy.RasterToNumPyArray(raster, nodata_to_value-9999) # 获取栅格信息 desc arcpy.Describe(raster) origin_x desc.extent.XMin origin_y desc.extent.YMin cell_w desc.meanCellWidth cell_h desc.meanCellHeight rows, cols arr.shape # 计算每一个像元的中心点坐标 # 注意numpy行号向下对应Y减小列号向右对应X增大 x_coords origin_x (np.arange(cols) 0.5) * cell_w y_coords origin_y (rows - np.arange(rows) - 0.5) * cell_h X, Y np.meshgrid(x_coords, y_coords)上面的X和Y就对应每个像元中心点的平面坐标。如果栅格是地理坐标系比如WGS84或CGCS2000那这些坐标就是经纬度。接下来有两种做法。一种是把坐标点批量转成点要素再交给arcpy的空间连接工具处理。另一种是直接用arcpy的几何对象做判断。批量创建点要素在大数据量下也不是特别高效所以更常见的是配合空间索引或者只用必要的点做判断。5.3 用geopandas shapely做空间关联如果电脑上装了Python的geopandas库有另一套更简洁的路线把栅格的每个像元中心点构造成GeoDataFrame把行政区shp读成另一个GeoDataFrame然后直接用空间连接函数sjoin一次完成关联。这个方法不需要ArcGIS许可跑起来也比较快尤其适合批量处理多个栅格文件。import geopandas as gpd import pandas as pd import numpy as np from shapely.geometry import Point # 假设已经用前面类似的方式得到 x_flat, y_flat, value_flat points pd.DataFrame({ x: x_flat, y: y_flat, value: value_flat }) # 把经纬度坐标转成点要素 gdf_points gpd.GeoDataFrame( points, geometrygpd.points_from_xy(points.x, points.y), crsEPSG:4326 ) # 读取行政区面并统一坐标系 gdf_adm gpd.read_file(rC:\gis_data\county.shp) gdf_adm gdf_adm.to_crs(EPSG:4326) # 空间连接每个点匹配所在的行政区面 joined gpd.sjoin( gdf_points, gdf_adm, howleft, predicatewithin ) # 保存结果 joined.to_file(rC:\gis_data\points_with_admin.shp, encodingutf-8)这段代码里最需要留意的是内存。如果像元数量超过几百万一次性把全部点放进GeoDataFrame再sjoin内存会吃得非常狠机器不好的会直接崩。我实际处理时是按行数分块的比如每50万行读取一次栅格切片转成点关联完保存一批再处理下一批。这样做虽然代码稍长但稳定。5.4 两种脚本方案的选型建议我自己比较习惯的分界线是这样的如果栅格数量不大但行政区面很细碎用arcpy方案更顺手因为可以直接复用ArcGIS自带的空间关系和坐标系管理能力。如果栅格本身很大或者你有多个栅格要批量关联用geopandas方案更好因为它写起来简洁而且不依赖ArcGIS许可可以在自己的服务器或笔记本上跑。无论是哪种脚本方案开工前都建议先做一次小范围的测试裁剪一小块栅格跑通全流程确认坐标、字段、输出都对再放全量数据去跑。否则全量跑完才发现坐标系没对齐几小时就白费了。6. 常见问题与排查技巧实录6.1 关联结果发生偏移该怎么办如果关联出来的点落到相邻的行政区里或者反过来无法匹配第一件事检查坐标系是否完全统一。在ArcGIS里把点图层和行政区图层都打开属性列表里的源标签看坐标系到底是不是同一个。如果一个是WGS84经纬度一个是CGCS2000投影坐标它们在地图上可能显示位置差不多但精度完全不一样空间连接的结果就不可靠。建议全部统一到一个坐标系再做关联。多数情况下我统一到CGCS2000地理坐标系单位是度便于直接用经纬度描述结果。如果坐标系看起来是统一的但偏移仍然存在检查是否做了定义投影而不是投影转换。最简单的验证方法是在原始栅格的属性中看范围和坐标值假如文件是WGS84经纬度范围应该大约在73度到135度、18度到54度之间如果文件被错误定义为CGCS2000投影坐标范围可能会变成几百万米的巨大数值。这种情况只能先把坐标系定义改正再重新投影。6.2 栅格转点后数量对不上少了或多了一些点点数量比像元总数少最常见的原因是NoData。栅格里的NoData像元不会转成点如果你不确定哪些是NoData可以在栅格属性里看统计信息或者用栅格计算器把NoData像元单独提取出来看一下。如果NoData很多先做填洼、插值或者其他预处理再去转点。点数量比预期多可能是因为栅格被重采样过像元大小发生变化或者输入栅格是多波段误把多个波段都转了。使用工具时注意输入必须是一个单波段栅格并且输出的像元大小和原始栅格一致。6.3 行政区代码前导0丢失导致无法关联行政区代码如果是字符串面转栅格时如果字段被自动转成整型前导0就会丢。最典型的是某些地区代码以0开头比如有些乡镇级代码是010开头转成整型后变成了10后面想用这个代码关联回行政区属性表就完全对不上了。解决办法有两个一是在面转栅格之前用添加字段字段计算器生成一个整型的唯一ID面转栅格时用这个ID作为Value字段后面统计结果再用ID关联回行政区属性表取得完整代码。二是如果一定要保留字符串代码建议不要转栅格而是走栅格转点的方案因为点属性表里字符串字段不会丢前导0。6.4 数据量太大卡死如何优化如果栅格像元数量在百万级以下栅格转点和空间连接还能撑住一旦到千万级以上就要考虑三点优化。第一用环境设置里的处理范围先把分析区域裁剪到目标行政区内不要全图跑第二先用行政区面的外接矩形裁剪栅格缩小数据量第三用脚本方案按行分块处理避免一次占用太多内存。另外空间连接时尽量精简输出字段。默认情况下Spatial Join会把连接要素的所有属性都带过去如果行政区面有几十个字段输出点要素的表会非常臃肿拖慢后续操作。在工具参数里可以设置字段映射只保留行政区代码和名称其他全删掉运行速度会明显提升。6.5 边界上的点归属判断不一致落在行政区边界线上的点用INTERSECT匹配时结果不确定尤其两个面共用一个边界时。ArcGIS的JOIN_ONE_TO_ONE默认会选一个但具体选哪个不完全可控。为了结果稳定我习惯先用按位置选择或脚本把这种边界点筛出来再用最近距离归属到一个面。这里可以设置一个很小的空间容差比如1米让边界点优先匹配到面积更小或者行政中心所在的那一侧具体需要结合项目规则判断。6.6 常见问题速查表现象可能原因排查方向关联后大量空值坐标系不统一检查点图层和面图层坐标系关联结果偏移到相邻区域做错了投影对比栅格范围确认坐标值是否符合预期转点数量少于像元数存在NoData统计NoData数量确认是否可忽略行政区代码对不上前导0丢失检查面转栅格字段类型改用ID关联工具运行卡死数据量过大裁剪范围、精简字段、改脚本分块处理边界点归属不确定匹配规则问题改用CLOSEST或设置空间容差关于这类需求我个人的习惯是只要不是必须保留栅格成果我都优先走栅格转点空间连接这条路因为点数据随时可以核对、转换、导Excel出了问题也容易定位。如果最终非要栅格成果我一般先把行政区面转成栅格像元大小和范围都在环境设置里锁死再基于这个栅格做后续操作这样每个环节都在自己掌控中。坐标系方面新做的项目我一般统一到CGCS2000再开始处理。这个流程踩过不少坑之后总结出来的希望对做类似需求的同行有帮助。