做地理信息这行经常会被“地学大数据”这个词唬住。听起来很宏观落到具体项目上核心其实就五件事数据从哪来、坐标系对不对、处理的顺序合不合理、成果能不能顺利导出、导出之后别人能不能正常打开。这段时间我在整理山东区域的DEM、POI和NDVI数据时把这五件事挨个折腾了一遍中间还干过一件糗事——把Word模板里的图表数据改了之后生成的文档怎么都打不开。这篇文章把完整流程、参数选择和踩坑记录都写出来给正在捣鼓GIS工具和地学数据的小伙伴做个参考。1. 从山东30米DEM下载到NDVI出图这几类数据的源头先把住1.1 地理空间数据云的通用下载路径DEM数据下载这件事绕不开地理空间数据云。我这次需要山东省30米分辨率的DEM进去之后在“数字高程模型”分类下能找到SRTMDEM 30M和ASTER GDEM 30M两个常见选项两者都用一种叫“图幅”的方式组织数据。页面上会画出一张经纬度网格每一格对应一个下载瓦片山东横跨的经纬度范围不小常规得勾选好几片相邻瓦片再一次性批量下载不然后面拼接时会发现边界缺一块。下载之前尽量先注册好账号登录状态比游客模式下载效率高得多文件也会完整很多。拿到手通常是一个压缩包里面是IMG或TIF格式的高程栅格。我的习惯是先把压缩包解压到全英文且不带空格的目录里文件名就保留原始标识不改成“山东DEM最终版”这种中文长名。原因很简单ArcGIS、QGIS或者GDAL在读取中文路径时偶尔会出奇怪问题不是找不到文件就是编码乱码排查起来相当浪费时间。下载完成后不要急着拿去算坡向、算河网先做两件基础工作。第一把涉及的瓦片统一加载到GIS软件里肉眼看一下瓦片之间的重叠和高程范围有没有突变。第二在栅格属性里确认坐标系SRTMDEM 30M常见的初始坐标系是WGS84地理坐标也就是经纬度格式。先知道这一点后面处理和投影时心里就有数。要是你只是做一张好看的背景图直接用原始经纬度坐标也没多大问题但一旦涉及到面积统计、格网裁剪、叠加分析和测距就必须先把坐标系处理到位。1.2 除了DEMPOI和遥感影像的数据源常识再说POI数据集。很多人以为POI就是从某个地图平台直接导出的表格拿来就能用。实际拿到的POI往往带一堆隐藏问题坐标系可能是加密偏移后的火星坐标或百度坐标类目字段混乱到一塌糊涂重复数据多到让你怀疑人生。POI最典型的来源是地图开放平台你导出时会看到名字、地址、经纬度、类目字段这些字段看似齐全但经纬度到底基于哪个坐标系不同平台并不一致。如果项目里还有别的空间数据要一起叠加分析这步不处理好后面所有结果都不可信。NDVI数据侧又是另一套逻辑。NDVI是归一化植被指数核心输入是含红波段和近红外波段的遥感影像。常见可选的有Landsat系列、Sentinel-2、MODIS等等。我这次用的是30米级别的影像和DEM分辨率恰好能对应上。下载这类影像时建议优先选做过大气校正的L2级产品少很多预处理麻烦。如果没有直接现成NDVI产品就需要自己从原始波段计算计算前看清楚元数据里的波段编号和缩放系数否则算出来的植被指数范围会离谱到没法用。2. DEM预处理链投影、填洼、坡度顺序错了结果差别极大2.1 统一投影坐标系这一步别偷懒我项目里用的DEM原始数据是经纬度坐标直接拿去计算坡度坡向会得到一个让新手完全摸不着头脑的结果。原因是经纬度不是等距坐标系纬度方向上的一个“度”和经度方向上的一个“度”在实地距离上差别很大。如果不先投影坡度算法会把经纬度当作普通平面坐标去算高差和距离算出来的坡度值几乎全是错的。这一步我的习惯是先把DEM统一投影到适合山东区域的投影坐标系。山东大体在东经114度到122度北纬34度到38度之间常用CGCS2000三度带高斯投影或UTM 50N都能很好覆盖。在QGIS里操作就是选中图层右键导出指定目标坐标系用GDAL就是一句gdalwarp命令加上-t_srs参数指定投影后的坐标系。投影完成后顺手检查一下像元大小原本的经纬度像元大概零点零零几度投影后应该变成30米左右的规则栅格。如果像元尺寸变成了一长条或特别奇怪的数值多半是投影参数没选对。2. 2 填洼、流向、汇流累加的处理顺序DEM拿到手之后很多人第一步就做坡度或坡向不能说完全错但如果你后面还要做水文分析、提取河网或者算汇水区这个顺序就会翻车。填洼必须放在流向分析前面原因是原始DEM里有很多凹陷区域可能是真实地形也可能是数据本身的高程噪声。如果不把洼地填平水流方向会在这些地方断掉后面的“流向”栅格和“汇流累积量”栅格会碎成一格一格的错误结果。正确的顺序是先填洼再计算流向再做汇流累积量。填洼工具在各类GIS软件里一般叫Fill Sinks或Fill参数里有个“最大填充阈值”默认值通常足够用但如果你的区域里有明显的道路、桥梁、水库这些人工地形就要谨慎确认。山东中西部有大量平原和农田DEM分辨率只有30米农田边界、土坡和人工沟渠本身就可能造成一些看似洼地的噪声点直接大刀阔斧填掉未必符合实际情况。还有一点特别容易忽略提取河流或汇流累积量之前最好先把栅格裁剪到研究区边界同时让河流边缘在边界处断干净。如果不裁剪汇流累积量会把边界外的大范围地形全部算进去最后生成的“河流”在边界处歪歪扭扭明显不自然。我这次是把填洼后的DEM裁剪到山东省界范围再做流向和汇流最后提取河网效果就正常很多。2.3 坡度坡向计算的参数陷阱坡度坡向计算时有一个特别经典的坑叫z系数。如果DEM还保持经纬度坐标直接在GIS软件里计算坡度软件会提示或要求你输入z因子含义是“高程单位与平面单位之间的换算比例”。假如你的高程是米平面单位是度那z因子至少得设成实际数值左右才能让坡度勉强像个样子。但不同纬度、不同尺度下这个系数并不恒定非常容易出错。所以最稳妥的做法就是我前文说的先投影成米制坐标系再算坡度。投影之后平面单位和高程单位都是米z因子默认变成1软件不会再乱猜。这一步看似简单实际上能排除至少一半的DEM计算异常问题。投影之后再表面分析生成坡度、坡向、山体阴影才能得到符合预期层次的结果。3. NDVI数据实操缩放系数、云掩膜与行列裁剪3.1 拿到影像先看元数据别急着算比值NDVI的计算公式几乎人人都会写近红外波段减红波段再除以近红外波段加红波段取值在-1到1之间。可很多人拿到影像后直接把原始DN值套进公式算出来的NDVI范围从负几千到正几千完全没法用。原因就在波段数值的“物理含义”上。以Sentinel-2 L2A数据为例它的红波段和近红外波段虽然是反射率但储存值往往是原始数值乘以10000。也就是说你在栅格属性里看到的3000、5000这样的数字要除以10000才是真实的反射率。哪怕有的影像看起来已经是浮点型也建议先读元数据里的“Scale factor”字段确认。Landsat系列也类似L2级产品的反射率通常带一个固定的缩放系数。拿到一个陌生的影像文件第一步绝对不是打开栅格就开始算而是看元数据或XML文件搞清楚波段编号、单位、有效值范围和缩放系数。另一个注意事项是波段序号。Landsat 8红波段一般是B4近红外是B5Sentinel-2红波段是B4近红外是B8。要是用红和近红外之外的波段算出个像模像样的结果那真的是“方向错了跑得越快偏得越远”。我一般在脚本里写一个辅助函数先打印所有可用波段名和对应波长再决定后续计算用哪个波段的索引。3.2 云掩膜和异常值怎么处理NDVI计算最容易出的问题不是公式而是数据质量。遥感影像里云、云影、雪还有部分传感器异常值都会让NDVI结果出现极端异常。山东这种区域春天部分时段有云和薄雾直接用有云影像去统计植被覆盖得出的结论往往跟实际差别很大。处理思路有两个一个是用影像自带的云掩膜波段另一个是根据云量条件选影像。以Sentinel-2 L2A为例通常有QA60波段能标识云和卷云掩膜Landsat L2则可能有质量评估波段数值对应不同质量。把云像元判定出来后把这些区域赋予“无数据”或者直接排除再做NDVI才会得到干净指数。有些影像为了标记无数据会使用-9999、0之类的无效值计算时必须把这些值相应处理掉不能让它们参与“近红外加红波段”的运算。处理完后还要做一步裁剪。计算NDVI通常是在整个幅影像范围上做但研究区只占其中一部分。这时候把结果裁剪到山东省界或者具体县界能减少无效计算量也让后续分区统计跑得快。如果研究区跨多个影像幅需要先把两幅相邻影像拼接或镶嵌再统一裁剪。裁剪完成后检查一下空值区域的位置避免把“研究区外的无数据”误读成“植被指数低”。3.3 用NDVI统计区域植被覆盖状况NDVI算出来以后很多人会直接拿它出图。但我个人建议先做一个分区统计比直接上渲染图更有说服力。做法是把NDVI栅格按行政区划边界做zonal统计计算每个区域的平均NDVI、中位数、标准差、像元数输出成一张属性表。这样项目汇报时就能直接说“某县平均NDVI是0.42比邻县高多少”而不是笼统地指着一张图说“这块颜色比较绿”。我这次还在NDVI基础上做了植被覆盖度估算。常见的近似公式是植被覆盖度等于NDVI减去裸土NDVI再除以纯植被NDVI减去裸土NDVI。裸土和纯植被的NDVI阈值可以根据区域实际情况定也可以取统计直方图的5%和95%分位数。这是一种简化处理但对于大多数应用场景已经足够。统计完成后结果可以用自然断点法分层设色出图山地和平原的层次一下就出来了。4. POI数据集的空间化坐标纠偏、去重、类目归并的实践记录4.1 先解决坐标系不是WGS84的问题POI数据集从来不是拿到就能直接画的点。常见地图开放平台导出的POI经纬度坐标系一般是GCJ-02也就是大家常说的火星坐标部分平台还可能输出BD-09。这两个坐标系和标准WGS84之间不只是一个常数偏移而是一种非线性偏移尤其在山东沿海地区直接叠加到WGS84底图上点会明显偏离道路和建筑几十米到几百米。解决办法有两个方向。第一在数据导出时选择“标准GPS坐标”或“WGS84坐标”字段这样后续最省心。第二如果导出时没得选就需要自己写转换工具比如从GCJ-02转到WGS84或者从BD-09转到WGS84。网上有很多现成算法原理是先做经纬度到平面坐标的投影再用一个迭代逼近方式反解偏移误差。我这里更推荐直接调成熟的地理库少自己造轮子。转换完必须做一次范围校验。山东的经度范围大致在114到123度纬度在34到39度之间凡是不在这个范围内的记录基本可以判定是坐标系转换失败或者地址录入错误。另外还要看单条记录是否落在江河湖海等明显异常区域有这类问题的数据打上标记后续核对。4.2 去重的土办法缓冲区和属性权重POI数据集另一个让人头疼的问题是重复记录多。同一家餐饮店可能被多个信息源都收进去名字写法还不太一样坐标相差几米到几十米。如果不去重统计“某区域有多少家餐饮”时结果会明显虚高。去重我一般分两步走。第一步是几何去重给所有点做一个几十米的缓冲区然后把缓冲区有重叠的点找出来当作候选重复集合。第二步是属性去重在候选集合里比对名称相似度、地址相似度和类目字段只要名称相似度高且坐标落在缓冲区范围内就保留质量最好的一条记录其余标记为重复。这种土办法在企业级POI清洗中比单纯按“名字完全相同”去重要可靠得多因为现实里同名不同址、同址不同名的情况很常见。这里我不建议直接用“名称相似度”作为唯一条件否则会把“李记饺子馆”和“李记饺子馆分店”都误判成重复最好把几何距离和名称相似度两个维度组合起来判断误杀率会低很多。4.3 类目归并与字段设计导出的POI原始类目往往是三级甚至四级的嵌套文本比如“美食-中餐厅-川菜”这种结构。这种原始字段适合展示但不适合做空间统计分析。我一般会把它拆成一级类目和二级类目两列比如一级类目统一定位为“餐饮”“购物”“住宿”“医疗”等二级类目再细分。遇到一个点同时归属多个类目的情况优先保留其主要功能类目避免重复统计。字段设计上也要提前想清楚。经过清洗后的POI数据我会统一保留以下字段标识号、名称、经度、纬度、省份、城市、区县、一级类目、二级类目、来源平台、更新时间。在GIS软件里加载CSV时注意指定经度字段和纬度字段并设置正确的几何坐标系导入后顺手转成与项目底图一致的目标坐标系。这一步完成POI数据才算真正能用于空间叠加、核密度分析、最近设施查找等GIS操作。5. 面积平差、Word图表模板以及“改完打不开”的排查记录5.1 面积平差工具的思路与公式GIS面积平差是一个非常日常但又特别容易糊弄的需求。比如我用山东的行政区边界去做地块统计把所有地块图斑面积加总得到的值跟该区的实际控制面积往往对不上。原因很直接图斑边界手绘精度有限投影变形也存在于边缘地块再加上拓扑编辑过程中的微小误差这些误差累积起来就是面积差额。平差的经典思路是按比例分配。先算出一个比例系数等于“控制面积”除以“所有图斑汇总面积”然后用每个图斑的原始面积乘以这个系数得到调整后的面积同时把微量差额记在最大图斑上保证调整后总和与控制面积严格一致。这个操作在ArcGIS里可以用“计算几何”加字段计算器完成或者在QGIS里用字段计算器写表达式本质上就是一个Python或表达式脚本。我这次的做法比较直接先在属性表里算各图斑原始面积然后用汇总工具得到总面积手动计算比例系数再新增一个“平差面积”字段填入“原始面积乘以比例系数”的表达式。计算完成后用“汇总”再验证一下平差面积总和。如果差额还在就把最后一点余数调整到面积最大的图斑上。这种土办法没有花哨界面但结果稳定、容易解释给甲方看也说得清。提到“面积平差工具”时我还想说一句别迷信插件很多问题用属性表字段计算就能解决关键是理解控制面积和汇总面积的关系。真正要写成工具的只是把这几步串成一个自动化流程输入是图斑要素和控制面积字段输出是平差结果字段。5.2 把POI成果批量生成Word的正确姿势项目成果里POI分类统计表、DEM渲染图、NDVI分区统计表往往都要整理进Word文档。我第一次做这类批量报告时直接手动复制图片和表格效率低不说图例、表格格式还容易乱。后来改用脚本生成才打开了新思路。核心思路是先用matplotlib或其他绘图库把所有空间分析结果渲染成图片再用python-docx库生成Word把图片和统计表按顺序插入模板。整个过程可以写成脚本每次数据更新后重新跑一遍就能稳定生成一份格式统一的报告。这么做的好处是图片样式、表格边框、标题字体都可以在脚本里统一控制不会出现上次手工调了半天、这次又全部重来的局面。在生成Word时还有一个容易被忽视的细节如果文档里需要“看似动态”的图表建议直接用matplotlib生成一张静态图片插进去尽量不要用Word内嵌的原生图表对象。原生图表看着高端但脚本操作它的数据源要解析docx内部XML复杂度一下子提高很多而且很容易把文件搞坏。5.3 “改完打不开”的根因与自查清单回到文章开头提到的糗事。我当时想在一个Word模板里直接修改内嵌图表的数值为了不重新排版就去改了docx包里面的图表数据节点。改的时候挺顺利保存后也没报错但双击打开Word时系统直接提示文件损坏要我选择是否修复。那种感觉非常糟糕因为文档内容是好好的却无论如何就是打不开。事后排查了很久总结出几个根因。第一Word的docx本质上是个zip压缩包内部有严格的目录结构图表数据既存在chart XML里也存在内嵌Excel的缓存数据里。你只改了图表XML而没有同步内嵌Excel缓存Word打开时会发现两处数据不一致判断文件损坏。第二如果你解压后改了文件又重新压缩打包格式稍微不规范比如压缩算法选错、根目录多了一层Word也会直接拒认。第三中途操作时产生了临时文件或锁文件没清理也可能导致解压出来的docx结构异常。最后我的解决方案其实很朴素彻底放弃“直接改Word内嵌图表数据”这条路线改成“脚本重新生成Word”。反正原始模板的样式和数据都在每次用python-docx按模板生成新文档该插入的图插图片该填的表格用文本写入完全不用碰XML。这一套方法看起来没那么“高级”但极其稳定已经连续多个月没再出现打不开文件的情况。如果大家已经有“打开即提示修复”的损坏docx我建议先复制一份备份然后用解压工具检验zip结构再用文件对比工具找出损坏的xml片段如果受损严重直接放弃原始文件用模板重新生成。日常操作时养成三个习惯能避免大部分事故一是别直接修改docx内部对象能加图片加图片二是每次修改前保留一份备份三是把生成的文档导入到WPS或Word里都试一次确认双端都能打开再交付。这样“改完打不开”的坑基本就能绕开了。
