ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

土壤空间数据成果包:mxd、shape、TIF三件套协同与校验实战

土壤空间数据成果包:mxd、shape、TIF三件套协同与校验实战 简介北京六种土壤类型空间分布数据包面向地理信息系统、城市规划、农业与环境管理等专业人士可用于土壤资源调查、农业适宜性评价、环境影响评估等场景。数据基于联合国粮农组织土壤分类体系划分制图单元涵盖土类、亚类、土属、土种等层级属性表中的土壤编号与编码表中的亚类一一对应便于精确识别与对照。资源包含标准矢量格式文件提供空间坐标与属性信息同时配有可编辑地图文档和标准成图文件前者方便用户灵活调整图层样式后者可直接用于打印或出版。资源包共二十个文件涵盖矢量数据文件、表格格式的土壤分类编码表以及地理标注文件、元数据等辅助文件压缩包整体约一百五十一兆。目前已有177人学习适合需要开展土壤空间分析或制作专题地图的中高级地理信息系统使用者。随包附带的显示样式修改示例支持自由选择配色方案进一步降低了制图门槛。1. 北京6种土壤类型空间分布mxd、shape、TIF三件套到底怎么配合拿到一个“北京6种土壤类型空间分布”成果包最值钱的不是那张标准成图TIF而是那个mxd可编辑文件和配套的标准shape文件。mxd负责把图层组织、符号、标注和出图版式固定下来shape负责回答“哪块地是哪种土壤”的矢量问题TIF负责对外分发和快速加载。对于做土壤调查、国土空间规划、农业区划、环境影响评价的从业者这套组合能省掉大量重建图层和重新制图的功夫。但实际落地时多数人第一步打开mxd就卡住数据源红感叹号、坐标对不上、面积统计跑偏。这篇笔记按验收工程文件、核对矢量属性、验证栅格成图、处理真实踩坑的顺序把这套成果包用到位。2. 先修mxd工程文件里全是路径引用不是数据本身2.1 mxd的红感叹号不是玄学为什么数据明明在却加载不出来mxdArcMap Document只记录图层名、路径、符号、标注、布局和工具设置本身不复制任何数据。ArcMap 打开工程文件时按路径去“找”数据找到就显示找不到就在图层名前挂一个红感叹号。这个行为让很多刚接触mxd的人误以为文件损坏了其实只是路径关系断了。最常见的断路径场景是成果包经过压缩、发送、解压后原来的目录层级变了。比如压缩包里的结构是soil_mxd/soil.mxd和soil_mxd/data/soil_type.shp你解压时只把soil.mxd拽到桌面shape文件留在原目录里那打开mxd必然满屏红感叹号。另一个常见原因是网盘同步时把.mxd和.shp分到了不同机器数据源自然失效。处理思路要分两层。第一层是防mxd在ArcMap里保存时勾选相对路径Customize → ArcMap Options → General → Make relative paths只要mxd和data的相对位置不变拷到任何机器都能打开。第二层是救已经断链的用Catalog一个图层一个图层重新指定数据源也行但图层一多效率很低。我处理这类成果包的习惯是直接写脚本重连10分钟之内把几十个图层的断链修完。提示拿到mxd后别急着双击打开先按原目录结构把整个成果包放好mxd和shape、tif的相对位置不变90%的红感叹号根本不会出现。2.2 用arcpy把失效数据源一次性重连脚本逻辑与参数说明如果你的mxd已经红感叹号了而且不想手工点可以按下面这个脚本批量重连。它假设你已经把所有shape文件和TIF统一放到了一个新的数据目录里。import arcpy, os mxd_path rD:\soil_beijing\soil.mxd data_dir rD:\soil_beijing\data mxd arcpy.mapping.MapDocument(mxd_path) results [] for lyr in arcpy.mapping.ListLayers(mxd): if not lyr.supports(dataSource): continue old_src lyr.dataSource # 数据源还能正常解析就跳过 if arcpy.Exists(old_src): continue base os.path.basename(old_src) candidate os.path.join(data_dir, base) if not os.path.exists(candidate): print(f未找到匹配文件: {base}) continue # shapefile按工作目录重连传不带扩展名的数据集名 if base.lower().endswith(.shp): lyr.replaceDataSource(data_dir, SHAPEFILE, base[:-4], False) # TIF/IMG 等栅格按RASTER类型重连 elif base.lower().endswith((.tif, .img)): lyr.replaceDataSource(data_dir, RASTER, base, False) results.append((old_src, candidate)) mxd.save() mxd.close() for old, new in results: print(f{old} - {new}) print(f共重连 {len(results)} 个图层)这段代码有几个参数要留意。mxd_path必须指向.mxd的绝对路径data_dir是所有 shapefile 和 TIF 的根目录脚本会用文件名去这个目录里寻找同名文件。replaceDataSource的参数含义分别是新工作空间路径、工作空间类型、数据集名称、是否校验。SHAPEFILE和RASTER是固定枚举值不能写成小写或中文。另外这个脚本依赖ArcGIS Desktop或ArcGIS Pro自带的Python环境普通Python解释器没有arcpy。如果你拿到的是地理数据库.gdb里的要素类不能直接用上面的写法。这时更稳妥的是lyr.findAndReplaceWorkspacePaths(old_ws, new_ws)把旧的gdb路径整体替换成新的gdb路径。我一般会在脚本里先判断lyr.workspacePath是文件夹还是gdb再分两条路走。重连完成后打开mxd目测一遍各图层是否正常显示这一步别省。2.3 符号和注记也得一起验6类土壤图例别等出图再核对mxd里真正容易出问题的不只是数据源还有符号系统和注记。土壤分类图对图例要求很高6种土壤通常用6种颜色区分颜色相近的两类如果出现在相邻图斑读图时非常容易混淆。拿到mxd后我建议在Contents面板里逐个展开图层确认符号是按SOIL_CODE字段做的唯一值渲染而不是单一符号渲染。如果mxd里的图例和实际shape字段对不上最典型的表现是打开后所有图斑一个颜色。这时要看图层属性里的符号系统重新指定唯一值字段再手动设置6个类别的颜色。注记部分则要检查字体是否缺失中文字体缺失会把注记显示成方框或乱码。这个坑在转移mxd到没有安装中文字体的机器上时特别常见。修复注记没有太好的批量办法通常是把 mxd 里的注记图层Annotation重新指向原有字体或者干脆在ArcMap里把注记转成图形Graphics再固定住。后者会丢失和属性的联动带属性标注的场景不建议这么做。3. 标准shape文件怎么核对字段、面积和拓扑缺一不可3.1 先认清shape的“四件套”和.prj坐标文件shape与其说是一个文件不如说是一组文件。标准shape文件至少包含.shp几何、.shx空间索引、.dbf属性表三个文件坐标准确的话还会带.prj投影坐标系定义。缺了.shx部分软件能通过重建索引打开但ArcMap会提示缺少文件缺了.dbf属性表直接丢失缺了.prj图形能显示但坐标系统未知。拿到“标准shape文件”这一步第一件事就是用ArcCatalog或者arcpy.Describe()查坐标系。北京地区的土壤数据最常见的是CGCS2000高斯克吕格投影比如CGCS2000_3_Degree_GK_CM_117E少数老资料用Beijing 1954。如果shape和mxd里的TIF坐标不一致叠加时就会出现一个在东一个在西的情况。检查代码如下import arcpy fc rD:\soil_beijing\soil_type.shp desc arcpy.Describe(fc) print(形状类型:, desc.shapeType) print(坐标系:, desc.spatialReference.name) print(投影:, desc.spatialReference.projectionCode) print(范围:, desc.extent.XMin, desc.extent.YMin, desc.extent.XMax, desc.extent.YMax)看输出时重点确认三件事坐标系是否带投影而不是裸的GCS_WGS_1984经纬度、范围是否落在北京附近的大致坐标区间、shapeType是否为Polygon。如果坐标系是经纬度但mxd里用投影坐标出图那所有面积统计都会出问题必须先用Project工具做统一。3.2 用Pyshp统计6类土壤面积字段选择与计算口径有了标准shape文件后第一件要做的事通常就是统计6类土壤各占多少面积。很多人习惯直接在ArcMap里打开属性表右键字段计算统计这一步没问题但要注意两点直接用经纬度度的图层算面积是错的算出来的是“度”的平方要用投影坐标系的shape才能得到平方米或公顷。大多数成果包的shp属性表里会自带面积字段比如AREA_HA。这个字段一般是建库时按投影坐标算好的直接累加即可。如果字段不存在或全部是0就先用Project工具投到CGCS2000高斯克吕格再添加到属性表算面积。用Pyshp读取并统计的方法如下import shapefile from collections import defaultdict path rD:\soil_beijing\soil_type.shp sf shapefile.Reader(path) # fields[0]是删除标记字段从1开始取真实字段 fields [f[0] for f in sf.fields[1:]] soil_col fields.index(SOIL_CODE) area_col fields.index(AREA_HA) soil_count defaultdict(int) soil_area defaultdict(float) for rec in sf.records(): code rec[soil_col] area float(rec[area_col] or 0) soil_count[code] 1 soil_area[code] area for code in sorted(soil_count): print(f类型编码 {code}: 图斑数 {soil_count[code]} f面积 {soil_area[code]:.2f} 公顷) total sum(soil_area.values()) print(f合计面积 {total:.2f} 公顷)这段代码依赖pyshp安装命令是pip install pyshp。字段名必须和dbf里的完全一致dbf字段名一般不带特殊字符和中文如果你打开看到拼音字段如SOIL_CODE、AREA_HA直接用。统计结果出来后要和TIF成图上的图例说明逐个对一遍。尤其是总面积如果和行政区总面积差的太多就要警惕是不是shape里包含了大量空白区或拓扑错误区域。注意Pyshp只读属性不解析投影。它拿到的坐标值原始形态所以统计面积时不要用坐标自己算。就用AREA_HA或提前投影后的字段。3.3 拓扑检查重叠、缝隙、碎面三个常见病土壤类型shape文件里最见不得三个问题一是图斑重叠二是图斑之间有缝隙三是碎小图斑过多。重叠会导致同一块土地被计算两次缝隙会导致统计面积偏小碎面则影响制图美观也在后续空间分析时制造大量额外处理。在ArcGIS里做拓扑检查可以这样操作新建File Geodatabase在里面创建拓扑规则选择“不能重叠Must Not Overlap”然后把shape导进去跑一遍。如果不想建复杂的Topology数据集也可以用下面这套简化流程先用Dissolve把同类图斑合并再和行政边界做Union对比。import arcpy env_ws rD:\soil_beijing\check.gdb fc rD:\soil_beijing\soil_type.shp arcpy.env.workspace env_ws arcpy.env.overwriteOutput True # 第一步按土壤代码融合 dissolved env_ws r\soil_dissolve arcpy.Dissolve_management(fc, dissolved, SOIL_CODE) # 第二步与县域边界做Union boundary env_ws r\beijing_boundary union_out env_ws r\soil_union arcpy.Union_analysis([dissolved, boundary], union_out, NO_FID) # 第三步检查SOIL_CODE为空的区域就是没有土壤类型覆盖的缝隙 fields [SOIL_CODE, AREA_HA] with arcpy.da.SearchCursor(union_out, fields) as cur: for code, area in cur: if code is None or code : print(缝隙区域面积:, area, 公顷)Union之后SOIL_CODE为空的记录代表行政边界内没有被任何土壤类型覆盖的区域基本可以认定为数据缝隙。碎面检查一般通过按面积字段筛选来看统计图斑数量后把面积小于一定阈值的记录挑出来人工判断是制图精度导致的“针眼”还是真实存在的细小地块。阈值常见做法是取1公顷或1万平方米具体看你的研究尺度。4. 标准成图TIF怎么验怎么用分辨率、背景值和掩码提取4.1 用gdalinfo完成TIF元数据自检分辨率、波段、NoData标准成图TIF一般是一张带地理参考的栅格图打开之前先看元数据能省去很多“加载后怪怪的”问题。最简单的检查工具是GDAL自带的gdalinfo命令行执行gdalinfo D:\soil_beijing\soil_type.tif关键信息看这几行Size is后面的宽高像素数Pixel Size 后面每个像元对应地面尺寸Band 1 Block下面的TypeByte表示单波段8位整型NoData Value标的背景值。土壤分类TIF通常像元值就是类别编号1代表褐土、2代表潮土以此类推。如果Pixel Size单位是米说明TIF已经带了投影信息如果单位是度说明是经纬度坐标系后续做面积统计分析前得先投影。再看有无NoData Value很多分类TIF把0或255作为背景如果背景没有标成NoData统计面积时会把这些巨大面积的“背景”也算进去污染结果。用Python读TIF更直观我用rasterio较多import rasterio from collections import Counter path rD:\soil_beijing\soil_type.tif with rasterio.open(path) as ds: arr ds.read(1) profile ds.profile nodata ds.nodata transform ds.transform cell_width abs(transform.a) cell_height abs(transform.e) print(波段数:, profile[count], 数据类型:, profile[dtype]) print(NoData:, nodata, 像元尺寸:, cell_width, cell_height) valid arr[arr ! nodata] for value, count in Counter(valid.tolist()).most_common(): area_ha count * cell_width * cell_height / 10000 print(f像元值 {value}: 数量 {count}面积约 {area_ha:.2f} 公顷)这里用transform.a和transform.e直接拿像元宽高前提是TIF没有旋转参数土壤成图TIF一般不会做旋转直接乘没问题。算出来的面积精度取决于像元大小像元是30米那每个像元是900平方米整体统计可以到公顷级别够用了。4.2 用面图层裁剪DEM和掩码提取有啥区别核心在输出边界与NoData这个是一个非常高频的问题手头有一个面图层比如研究区边界要对DEM栅格做裁剪到底该用Arctoolbox里的Clip工具还是用Spatial Analyst的Extract by Mask这两者在ArcMap里经常被混用实际行为有明显差别。如果选数据管理工具里的栅格Clip不勾选“使用输入要素裁剪几何”输出TIF的范围是输入面要素的最小外接矩形得到的结果是矩形图片面边界以外的区域依然保留像元值只不过范围被切小了。很多人误以为给了面就能按面形状切出来最后看到矩形就懵了。解决办法是勾选Use Input Features for Clipping Geometry输出就按面边界裁剪区域外像元变为NoData。而Extract by Mask从原理上就是“掩膜提取”输入一个面或栅格作为掩膜输出范围严格贴合掩膜边界区域外全部是NoData。它的优势不是“剪得好看”而是后续做像元统计或重分类时直接按有效区域计算不用再额外做一次按属性筛选。实践中两者怎么选如果是纯出图、只要边界范围内的图像我会用Extract by Mask如果要批量处理多个矩形区块、保留完整栅格信息方便后续重采样用Clip。还有一个细节用面裁剪DEM之后边缘会生成NoData像元如果做坡度分析边缘一圈会产生异常数据建议裁剪后先用邻域填充把边缘NoData修掉。4.3 TIF加载后一片黑或全灰检查colormap与符号化标准成图TIF如果本身是分类图通常会带Color Table调色板。ArcMap里加载后如果没有自动套用就会按灰度拉伸显示6类土壤全变成深浅不一的灰色看不出分类。此时右键图层 → Properties → Symbology选择Unique Values按像元值字段渲染再把6个类别分别指定颜色即可。在QGIS里更简单加载TIF后Layer Properties → Symbology在Render type里选Paletted/Unique values然后在下面点击Color map加载就能恢复到制图时的配色。如果加载后背景显示成大块黑色或白色通常是NoData没有设置为透明。ArcMap里在图层属性 → Symbology → Display NoData as transparent勾上即可。最后一个比较隐蔽的问题是成图TIF如果经过重采样或裁剪可能丢失Colormap。这时候显示出来就是单个灰带怎么调都调不回原配色。修复方式只能从原始shape重新用Polygon to Raster转一次或者用GDAL的pct2rgb.py等工具重新挂调色板。所以收到TIF时我建议第一件事就用gdalinfo确认有没有COLORMAP字段没有的话抓紧和发数据方核对原始制图文件。5. 落地踩坑实录mxd、shape、TIF协作时的5个典型案例5.1 满屏红感叹号数据源路径集体失效现象打开mxd后所有图层都带红色感叹号双击图层提示无法加载数据源但同一目录里shape和TIF都好好躺着。原因mxd里记录的是绝对路径成果包解压后目录变了或者mxd在网盘上被单独下载过。解决不要手工一个个图层指定数据源按2.2的脚本批量重连如果脚本跑不通先检查是否使用了ArcGIS自带的Python环境。重连后保存前务必勾选相对路径避免下次再断。5.2 shape和TIF边界明显错位一个偏东一个偏西现象把shp和tif叠加到ArcMap里两类数据边界不重合有的地方错开几百米甚至更远。原因最常见的是shape带.prj但TIF没有投影定义或两者坐标系统不一致如shape是CGCS2000高斯克吕格TIF是WGS84经纬度。解决用gdalinfo分别看两个文件的投影和Origin确认后把其中一个投影到另一个的坐标系用GDAL Warp或ArcGIS Project Raster都行。对齐后再跑一次Boundary对比直到Extent基本重合。5.3 统计面积比行政面积小了一大截现象把6类土壤的面积累加起来比北京行政边界总面积少了几个百分点。原因一是shape本身存在缝隙二是统计用的AREA_HA字段过期或计算口径错误三是TIF和shape使用的边界版本不同。解决先跑3.3的缝隙检查看SOIL_CODE为空的区域有多少再投影后用Geometry Calculator重算面积字段。如果缝隙超过预期确认是不是数据源在制图时把河流、湖泊等水域单独抠出去了这不算错误但要在说明里写清楚。5.4 TIF加载后全黑只有一个值现象tif在ArcMap里加载后全黑或看起来只有两种颜色完全看不到6类土壤的分布。原因分类图的值域只有1到6被ArcMap按RGB连续色带拉伸显示低值全映射到黑色或NoData没设置透明把大片背景显示为黑色。解决图层属性 → Symbology → Unique Values按“Value”字段渲染手动填充6类颜色再把NoData设为透明。如果是QGIS在Symbology里选择Paletted/Unique values模式最省事。5.5 注记和文字全变成方框现象mxd里的土壤类型名称、图例文字都显示成方框数字正常但中文乱掉。原因字体缺失或字体名称在系统里不存在。ArcMap在打开mxd时会按字体名寻找找不到就用默认字体替换中文字形变成空白方框。解决先装好常用中文字体宋体、黑体、微软雅黑再重开mxd如果不能装字体打开mxd后把Annotation选中转成Graphics这样文字就变成可移动的图形元素不再依赖字体。注意转完后文字不再随属性变化更新。6. 一致性复核一套脚本检查mxd、shape、TIF三者对不上6.1 复核脚本把三样东西放在同一套逻辑下检查我习惯在拿到此类成果包后先跑一段一体化检查脚本不打开ArcMap界面也能量化出问题在哪。核心检查三条线mxd所有图层数据源是否全部解析shape和TIF的Extent是否重叠以及TIF像素面积与shape累积面积是否在合理误差范围内。脚本如下import arcpy, os mxd_path rD:\soil_beijing\soil.mxd shp_path rD:\soil_beijing\soil_type.shp tif_path rD:\soil_beijing\soil_type.tif result [] # 检查1: mxd数据源是否完整 mxd arcpy.mapping.MapDocument(mxd_path) missing [lyr.dataSource for lyr in arcpy.mapping.ListLayers(mxd) if lyr.supports(dataSource) and not arcpy.Exists(lyr.dataSource)] result.append((mxd数据源, PASS if not missing else f缺少{len(missing)}个)) # 检查2: shape与TIF范围重叠 shp_desc arcpy.Describe(shp_path) tif_desc arcpy.Describe(tif_path) shp_ext, tif_ext shp_desc.extent, tif_desc.extent dx min(shp_ext.XMax, tif_ext.XMax) - max(shp_ext.XMin, tif_ext.XMin) dy min(shp_ext.YMax, tif_ext.YMax) - max(shp_ext.YMin, tif_ext.YMin) result.append((shp-tif范围, PASS if dx 0 and dy 0 else FAIL)) for name, status in result: print(f{name}: {status})6.2 结果怎么读复核不是挑刺是给自己留后路跑完这段如果mxd数据源显示PASS再打开mxd基本不会被红感叹号绊住shp-tif范围PASS说明两套数据在空间上对得上后续做叠加不会穿帮。如果显示FAIL先别急着修数据回头查是不是gdb里的辅助图层被删了或者TIF经过裁剪后范围本身就小于shp范围。我自己的习惯是每收一套这类成果包不动原始文件先复制到工作目录再做一次复核把mxd重连、shape统计、TIF元数据三条线全部过一遍再提交下游。这套流程能拦住八成的翻车也让拿着数据做分析的人心里有谱。数据成果基本可以信但还是要自己校验一轮。希望帮到你。本文还有配套的精品资源点击获取
返回列表