ARTICLE DETAIL

资讯详情

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

上海市12.5m DEM数据处理全流程:投影、裁剪与踩坑实践

上海市12.5m DEM数据处理全流程:投影、裁剪与踩坑实践 简介面向GIS从业者及城市规划、环境研究人员这份数据包提供上海市12.5米分辨率数字高程模型并附行政边界矢量文件。借助格网高程与shp边界可以完成坡度坡向提取、地形曲率计算、地形渲染、汇水与水流路径模拟、日照时数估算及灾害风险评估适合GIS课程实践教学、城市规划前期分析、环境影响评价、防汛与地质灾害研究等场景。压缩包内含15个文件除核心tif和shp外还有tfw、prj等坐标参考文件dbf属性表、ovr金字塔、adf辅助栅格、xml元数据等便于在ArcGIS、QGIS等专业软件中直接调用整包大小约48MB。数据采用WGS84坐标系兼容性好便于与其他地理数据配准叠加。目前已有806人学习下载对需要高精度市域地形并结合行政边界开展研究、制图或教学演示的读者尤其具有实用价值。1. 上海市DEM 12.5m高程包一张能直接拿去分析的地形底图做城市内涝模拟、三维选址或者基站视距分析的人大概率都经历过被30m、90m分辨率DEM“磨平”的无奈河道拐弯糊成一条粗线几米高的小土丘直接消失。上海这种地势平缓的地区三五米的起伏就能决定雨水往哪边流地形分析对分辨率格外敏感。这份12.5m分辨率的上海市DEM数字高程数据改写了这个尴尬——它正好卡在30m和5m之间兼顾覆盖范围和细节还自带上海市范围shp文件省去了自己找行政区边界、再对齐数据范围的功夫。拿到手就能按区县边界裁剪、重投影、算坡度填洼。这套流程不只是打开看一眼而是要实实在在地落到项目里所以下面几章会按“先看数据、再裁剪、后投影、边踩坑边出活”的顺序把整条路径讲透。2. 先看清数据再动手DEM与shp的坐标系、范围和文件组织2.1 先看数据组成12.5m DEM和shp到底放在哪里解压之后先别急着拖进软件。我一般的做法是先打开压缩包目录看清楚里面到底有哪几类文件。常见组织方式是一个tif或img格式的DEM栅格文件加一套上海市区县或街镇边界的shp矢量文件偶尔还会带一个坐标系说明的txt或xml元数据。栅格文件通常只有一个主文件但tif格式会伴生tfw世界文件、ovr金字塔文件如果少了tfw加载时地理参考可能会丢失。shp则是一组文件主文件shp、索引shx、属性表dbf还有可选的prj投影文件。dbf决定属性表能不能正常打开prj决定坐标系描述两个谁都不能少。数据本身的组织方式会影响你先做什么如果DEM是分幅的比如按图幅切成多块tif那第一步应该把分幅数据合并成一个完整覆盖市的栅格如果只有一个整幅tif就省事多了。合并的常规做法是用ArcGIS的栅格数据集镶嵌工具或者QGIS里的合并工具把相邻图幅拼起来。有一个容易被忽略的点分幅tif之间可能有重叠区镶嵌时缝合线处理不好会出现明显色差或高程断裂用带羽化的混合镶嵌比直接覆盖更稳妥。路径上建议别用中文目录ArcGIS老版本对中文路径有时表现出玄学问题报错都看不出原因。2.2 加载后第一件事检查坐标系与高程基准双击打开ArcMap或QGIS把shp和DEM都拖进去之后我习惯先看属性——右键图层打开图层属性切换到“源”或“信息”选项卡。这里的重点有两处一是坐标系二是像元大小、深度和NO DATA值。这份数据如果做得好shp和DEM的坐标系应该是一致的常见的是CGCS2000 3度带高斯投影中央经线121.5°E正好覆盖上海也有可能是WGS 1984 UTM 51N这两种投影在上海市区范围内偏差在几十米到上百米肉眼可能看不出但和GPS采集点或在线地图叠起来就会飘。如果发现shp是地理坐标系而DEM是投影坐标系必须先统一再裁剪。我一般推荐把矢量重投影到与DEM一致的坐标系因为后边还要拿这个shp去做掩膜提取矢量投影的成本远低于栅格重投影。重投影矢量用“投影”工具选好目标坐标系后输出一个新的shp别忘了指定正确的输入坐标系否则工具会用默认值给你瞎猜一个结果直接翻车。高程基准同样值得看一眼。国内陆形数据的高程基准通常是1985国家高程基准也有使用吴淞高程的旧数据。做淹没模拟时水位数据要和高程基准对应基准差了零点几米淹没范围就能从“没有”变成“淹了一层”。怎么确认基准看元数据最靠谱没有元数据就看属性里的深度单位以及和已知水准点的高程差值来反推。2.3 用全局和局部视图确认数据范围与重叠关系检查完坐标系还得确认数据范围是否覆盖整个上海市。方法是把shp置于最上层打开属性表的范围统计或者看数据框属性里的范围栅格的覆盖范围在“源”选项卡里也能看到四至坐标。如果DEM范围明显小于shp范围说明数据有缺口后面裁剪时会出现大面积空白。上海崇明岛、长兴岛这类江心岛区域经常在分幅数据里被漏掉裁剪前用全局视图晃一眼就能发现。还可以用交叉重叠分析来判断把shp转成面用“相交”工具和DEM的矩形范围做一次交叉人工看交叉面积占比。常规操作是直接目测但数据多批次时要学会用属性表面积统计。我处理这类数据时有个固定习惯先在地图上打开刻度尺测量崇明岛到中心城区的距离和网上地图对比如果吻合就说明投影关系没问题。这个动作快而且能帮助建立对数据精度的感性认识。确认了范围、坐标系、基准后面裁剪和投影才不会出露馅问题。3. 用上海市界shp把DEM裁出来掩膜提取与栅格裁剪的差异和实操3.1 掩膜提取还是栅格裁剪两个工具的工作逻辑裁剪DEM这一步最常见的选择题是“栅格裁剪”和“按掩膜提取”到底用哪个。先说结论如果只要矩形范围切一刀用栅格裁剪要的是跟shp边界完全吻合的不规则形状用按掩膜提取。很多新手在ArcMap里用“裁剪”工具裁出来的结果外轮廓还是方的边界外的像元只是被赋成了NoData并没有真正切掉。这些NoData黑边在后续镶嵌、坡度计算、转ASCII时都会制造麻烦。掩膜提取的逻辑完全不同它把你传入的shp当作一个掩膜输出栅格只保留落在掩膜内部的像元外部直接丢。对上海市这种临江沿海、边界曲折的地区掩膜提取出来的栅格边界是贴合岸线的后面和别的数据做空间叠加时不会出现越界像元。不过掩膜提取也有它的代价边界处经过重采样一个像元若有一半在界外可能被保留也可能被剔除这取决于工具默认的“中心点法”。所以做淹没模拟和坡度分析时边界一圈的高程值要留个心眼。还有一个小细节热词里经常有人搜“依靠面图层掩码提取和依靠面图层裁剪有啥区别”本质就是上面说的——裁剪是“切矩形”掩膜提取是“抠形状”。实际操作中如果你用面图层做裁剪ArcGIS会先用面的外接矩形去切再按面形状处理。这也解释了为什么同样一个shp两个工具产出的文件大小会差很多掩膜提取的往往小一圈。3.2 用ArcMap按上海市界裁剪DEM的分步操作打开ArcMap加载完整DEM和上海市区县边界shp之后我一般这么操作第一步打开ArcToolbox找到“Spatial Analyst工具—提取分析—按掩膜提取”双击打开对话框。 第二步输入栅格选择DEM输入栅格数据或要素掩膜数据选择shp输出栅格填一个路径和名称比如tif格式直接写“sh_dem_12_5.tif”。 第三步关键一步点开环境设置把“处理范围”设为“与图层shp相同”把“捕捉栅格”设为输入DEM。不设置处理范围的话输出范围可能按掩膜的外接矩形算又回到方形黑边不设置捕捉栅格的话像元网格可能和原始DEM错半个像元后面和别的栅格做运算时对不齐。 第四步点击确定执行。跑完后用“识别”工具在边界处点几下对比原始DEM和裁剪结果的高程值验证没有发生偏移。掩膜提取工具默认的输出像元大小和输入栅格一致但如果你在环境里改了像元大小比如强行改成10m工具会在范围内重采样这等于对12.5m数据做了一次插值。除非你有明确的统一分辨率需求否则别乱改重采样没办法提高真实精度。3.3 用Python批量裁剪区县ArcPy脚本与参数说明一次只裁一个上海市范围shp还不够很多时候要按区县出图比如给每个区的规划科单独交付一份DEM。这时就要上ArcPy脚本了。下面这段脚本是我常用的批量裁剪模板按shp属性表中的“区县”字段循环输出。import arcpy from arcpy import env from arcpy.sa import ExtractByMask # 基本路径和参数 arcpy.env.workspace rD:\shanghai_dem arcpy.env.overwriteOutput True dem rD:\shanghai_dem\sh_dem_12_5.tif # 裁剪前的整幅DEM shp rD:\shanghai_dem\shanghai_district.shp # 上海市区县边界 out_gdb rD:\shanghai_dem\output.gdb # 输出到文件地理数据库 # 环境设置处理范围对齐候选面捕捉栅格保证网格一致 arcpy.env.snapRaster dem arcpy.env.extent shp # 按属性字段“区县名”逐条掩膜提取 with arcpy.da.SearchCursor(shp, [FID, 区县名]) as cursor: for fid, name in cursor: where_clause FID {}.format(fid) out_name dem_ name .tif arcpy.Select_analysis(shp, in_memory/select_fea, where_clause) out_raster ExtractByMask(dem, in_memory/select_fea) out_raster.save(out_gdb \\ out_name) print(已输出, name)这段脚本核心动作有三个Select_analysis把单个区县要素单独拎出来当掩膜ExtractByMask按这个掩膜抠栅格最后save输出。需要注意where_clause里的字段名如果是中文要确保数据库连接层能正确读取否则会报“字段不存在”。另外snapRaster和extent两个环境参数必须在裁剪前设定好它们让每个区的输出栅格网格与原始DEM完全对齐不会一个区一个网格错位。如果你要把输出直接给别人建议tif格式而不是gdb栅格后者在跨软件交换时常遇到无法直接识别的问题。脚本跑完后打开输出目录检查几个典型特征每个文件大小应该跟区县面积成正比边界小区如黄浦区文件应该很小崇明区的覆盖范围应该包含江面岛屿。任何突兀的文件大小变化都值得进一步查看通常意味着属性字段选错或掩膜重叠。4. 12.5m分辨率该配什么投影从CGCS2000到高斯投影的转换参数4.1 为什么12.5m DEM必须重投影分辨率对投影的敏感度拿到手的数据如果坐标系是地理坐标系WGS84或CGCS2000经纬度格式那在ArcMap里看起来像平面图但实际上每个像元对应的地面距离在纬度方向上会变化。上海的纬度约31°经度1°对应约95km纬度1°约111km如果用经纬度直接算距离东西方向和南北方向的比例不一致一切基于距离的统计都会失真。最直接的影响是坡度计算GIS的坡度工具是基于每个像元的水平距离来计算垂直变化水平距离错了坡度值自然错。12.5m这个分辨率本身不高不低它对投影的敏感度比90m大得多——90m栅格一个像元覆盖约8100平方米误差百分之几还看不大出来12.5m一个像元只有约156平方米投影变形带来的几何误差已经能直接扭曲地形形态。所以拿到DEM和shp之后统一投影是比裁剪更前置的一件事。常见的目标投影有两种一个是WGS 1984 UTM Zone 51N一个是CGCS2000 3度带高斯投影中央经线121.5°E。上海的大部分范围正好落在UTM 51N带内和高斯3度带121.5°E这套几乎一致但两者基准不同UTM是WGS84椭球CGCS2000是自成一体的地心坐标系。转来转去时要选对地理变换方法否则坐标能差出几十米。对于只能接受一套坐标系的项目我推荐统一到CGCS2000 3度带因为国内规划院、测绘院的成果基本都是这套日后交接省心。4.2 投影转换的参数设置与重采样方法选择在ArcMap里用“栅格投影”工具Project Raster做重投影对话框里有几个参数必须搞明白。输入栅格选原始DEM输出坐标系选目标投影重采样技术有三项可选最近邻法、双线性内插、三次卷积内插。对于高程数据基本原则是不要用最近邻法去插值因为它会把某个像元的值整块复制让地形出现台阶感双线性是稳妥默认速度适中平滑度也可以三次卷积的保真度更高边缘会稍微出现振铃效应但高程数据本身是连续表面用它也不会出错。我实际处理12.5m数据时常用双线性原因很务实计算量小、边界不太过冲、后续坡度分析取值稳定。参数表我整理如下照着设置就行参数项推荐设置说明输出坐标系CGCS2000_3_Degree_GK_CM_120E上海所在3度带若跨带选择121.5E地理变换CGCS2000_To_WGS_1984 若从WGS84转缺失时可能导致几十米偏移重采样BILINEAR高程连续表面推荐双线性输出像元大小12.5保持原分辨率不要顺手改成整数环境处理范围与shp相同保证输出范围与行政区边界一致这里特别想提醒一个坑很多人转完投影后发现输出的tif像元大小变成了12.499999或12.500001这类不整的数字这是投影前后椭球体变形导致的正常现象。不要手动把它改成“看着整齐”的12强行修改会让栅格网格发生扭曲后续和原始数据叠加时错位。正确的做法是在工程环境里设置输出像元大小让软件按统一规则计算保持全范围一致。QGIS用户可以用“导出—另存为”对话框中的CRS设置勾选“重采样”并选双线性道理一样。Global Mapper用户也可以直接做重投影它的投影库更新得比较全很多老坐标系在ArcGIS里要额外安装坐标转换文件才能识别Global Mapper往往开箱即用。但不管哪个软件做完之后都要抽查一个已知坐标点比如人民广场把栅格上读取的值与坐标对照验证投影过程没有产生意外偏移。5. 踩坑记录坐标偏移、NoData黑边、范围对不上和属性表丢失5.1 坐标偏移几百米数据源或投影不一致引起的经典故障现象把裁剪后的DEM叠到在线影像图上山体轮廓、岸线和影像明显错位偏移量不是几米而是几百米。我遇到过拿上海市界shp去裁一个以WGS84经纬度存储的DEM直接在ArcMap里拖进去看时还好一做掩膜提取边界漂出足足两个街区的距离。原因两个数据集的基准面不一致。shp可能投影在CGCS2000而DEM是WGS84或者反过来。坐标系的名称在属性里都写着“GCS_WGS_1984”但元数据里没有注明基准转换参数软件不知道两者间需要一个地理变换于是硬把它当成同一个基准。解决先分别查清两个数据集的坐标系如果发现基准不同在重投影时明确设置地理变换。CGCS2000转WGS84的参数在ArcGIS里通常是“CGCS2000_To_WGS_1984”选择后偏移量能压到1m级别。如果找不到这个变换宁可另外找一份坐标系统一的shp也不要勉强手动平移数据那属于饮鸩止渴。这类问题发生的频率说实话比想象中高得多尤其在不同来源的数据拼盘时。5.2 裁剪后出现NoData黑边掩膜提取的边界像元处理问题现象用掩膜提取裁剪出来的DEM在边界处有一圈黑色的NoData带有的甚至出现大面积黑块叠到影像上时黑边盖住了岸线。原因一是掩膜shp在边界处与栅格像元不完全对齐导致部分像元被判定为“不在掩膜内”赋成了NoData二是如果有多个面要素要素之间存在微小裂缝裂缝处的栅格算出来是NoData三是输出格式用了不支持NoData的格式比如某些img设置导致无效区显示成黑块。解决用“按掩膜提取”前先对shp做一次“修复几何”操作把自相交和裂缝修掉。输出时选tif格式并在栅格属性里把NoData值明确设置为一个不会出现的数字比如-9999。如果已经裁出了黑边可以用条件函数Con把NoData重新赋值为周围高程的均值但尽量少用这种后悔药因为它会引入假高程。最靠谱的办法还是回到裁剪源头环境设置里把“捕捉栅格”设为原始DEM让掩膜提取按原始网格计算黑边基本能消除。5.3 范围对不上DEM和shp的bbox不一致导致白板现象跑完裁剪输出结果只有一个灰黑矩形完全没有高程纹理或者东西南北范围明显短一截。属性表里看范围发现DEM的覆盖范围根本包不住上海市界某一侧直接被切掉。原因数据拼盘时没有复核数据范围。这种情况常发生在用旧版的上海市边界shp配新版的DEM上或者数据本身只覆盖了部分行政区比如只有浦东、青浦等几个区的合集拿来当全市范围用。解决加载数据后先做一个快速范围目视检查。比较稳妥的做法是打开shp属性表查看几何范围字段与DEM的“源”选项卡里的范围对比。发现缺口时别急着裁剪先补数据源把缺失区域的DEM额外下载再镶嵌或者找一个覆盖范围完整的上海市界shp。这里有个热词场景经常出现——有人到处找DEM数据其实是没意识到自己手上的shp超出了DEM覆盖范围。只要把shp和DEM叠起来看一眼几秒钟就能发现问题很多白板输出都源于这一步马虎了。5.4 属性表丢失或字段乱码shp文件组不完整或编码问题现象加载shp后图层能显示图形但打开属性表发现是空的或者字段名变成“锟斤拷”之类的乱码。还有一种常见情况Open属性表时提示“未找到字段”工具批量裁剪时where子句中字段名怎么都匹配不上。原因shp文件的dbf属性表缺失或损坏是最直接的原因。或者dbf没问题但编码格式是GBK而ArcMap和QGIS的新版本默认按UTF-8读取中文编码错乱用ArcPy脚本时读不到属性值选不出区县名来。解决先补文件。检查shp同一目录下有没有同名的dbf文件没有就说明数据包不完整去找原始来源补一份。有dbf但读不到试着在QGIS里加载用“文件—属性—编码”切换为GBK或GB2312重新读取能正常看到中文后另存一份UTF-8的shp。如果只是要裁剪其实不依赖属性表也能干活把每个区县面手动选出来逐次裁剪或者用FID字段做循环这样绕开字段名和编码问题。这个技巧在紧急出图时非常管用。5.5 批处理循环里输出互相覆盖in_memory数据的生命周期问题现象ArcPy脚本批量裁剪多个区县跑完后发现输出文件只有最后一个区县的前面的全被覆盖或压根没生成。或者脚本执行到一半报错“已存在”手工删了重跑还是这样。原因循环里用了相同名称的临时要素或输出路径且没有清理中间数据。arcpy.Select_analysis把结果存到in_memory时如果不及时删除下一次循环同名写入会冲突。另一方面部分环境和后续脚本对地理数据库里的同名栅格管控更严格没有处理overwriteOutput状态的设置就出现覆盖。解决脚本开头显式设置arcpy.env.overwriteOutput True每一次循环里把in_memory要素存储为带循环编号的名字比如select_fea_{}用完再Delete_management。输出的tif文件名确保唯一建议把区县拼音或代码拼进文件名。另外一个实用习惯是给脚本加print打印当前处理的区县名出问题时能直接定位不用猜黑匣子一样看结果。6. 高程分析前的最后一道工序填洼、坡向计算与可视化验证拿到裁剪并统一投影后的12.5m DEM接下来最常做的事情是把它加工成可直接支撑决策的产品。我自己的固定流程是先做填洼再算坡度坡向。DEM里难免有坑洼——可能是真实地貌也可能是原始数据生成时的噪声。用填洼工具处理掉这些虚假洼地后坡度、流向计算才不会把水流堵在莫名其妙的地方。填洼时要注意填洼阈值设得越大地形越平滑但真实的小尺度地形会被抹掉。12.5m数据比30m数据保留了更多细小的冲沟和土坎阈值我一般控制在10m以内保留微地形。坡度和坡向两个工具直接基于像元邻域计算12.5m分辨率下算出的坡度比30m数据更锐利也更容易出现噪点。计算完可以做一个一致性检查随意取几个山谷线位置对比坡度图和原始DEM的等高线如果坡度高值都落在陡坎处说明计算没有异常。最后是可视化验证。把算好的坡度图叠加到影像底图或者用Global Mapper打开DEM做山体阴影渲染能快速发现数据里的异常值。山体阴影是最直观的黑匣子探测法任何突兀的高程跳变都会以黑白斑点的方式显现。建议再用“剖面图”工具拉几条穿过城市中心和河道的剖面线对照沿线的高程值是否合理比如黄浦江沿线应当是一条连续的低洼带而不是阶梯状断裂。这些检查全部通过这份12.5m DEM才算真正能扛事。如果数据源实际来自DSM也就是未去除地表建筑的高度做洪水淹没和通视分析之前还得先做建筑物高度剔除。最基本的操作是用形态学开运算做一轮滤波再跟人工修正的平坦区域对比。这一步不是必须但做城市分析时值得记得。每回处理新到手的DEM我都习惯先花小半天时间把坐标系、范围、基准、NoData全部查清楚再往下走省得后面返工。这算是干这行养成的老毛病了但对这份12.5m数据来说投入的时间绝对值得。希望这篇能让你少走一段弯路。本文还有配套的精品资源点击获取
返回列表