GIS遥感数据工程【免费下载链接】gdalGDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.项目地址https://gitcode.com/gh_mirrors/gd/gdal点击查看免费下载导读本文围绕 GDAL 栅格驱动体系中用于读取USGS ASCII DEM以及加拿大CDED数据产品的USGSDEM驱动展开系统讲解该驱动的格式识别、坐标系与地理配准解析、高程单位与 nodata 处理、多版本格式变体兼容等核心机制。读完本文你将掌握如何用gdalinfo、gdal_translate以及 Python 绑定高效读取 USGS DEM 文件并理解底层实现如何依据固定偏移解析文件头、如何判定数据类型与投影以及驱动以单一大瓦片方式载入数据的性能特征。驱动概览与历史背景USGSDEM是 GDAL 内置built-in by default的栅格只读驱动对应官方文档 drivers/raster/usgsdem.rst实现代码位于 frmts/usgsdem/usgsdemdataset.cpp。USGS ASCII DEM 是美国地质调查局USGS在被 SDTS 格式取代之前使用的传统栅格高程交换格式同时它也是加拿大CDEDCanadian Digital Elevation Data产品的承载格式。GDAL 对该格式的支持覆盖了市面上绝大多数流行变体包括正确识别文件中声明的坐标系地理坐标 / UTM / 州平面投影依据四角点坐标与像元间距完成地理配准定位正确标记边缘缺失数据区域的 nodata 值正确区分米m与英尺ft两种高程单位。代码渊源文档明确指出usgsdemdataset.cpp的读取逻辑部分源自 Ben Discoe 的VTPVirtual Terrain Project软件中的 USGS DEM 导入器importer源码头部注释亦注明了这一出处见 usgsdemdataset.cpp#L7-L8。此外文档还记录了该驱动的导出写入能力曾由加拿大 Yukon 省环境部门提供资金支持而开发不过从当前仓库源码结构看USGSDEMCreateCopy在 usgsdemdataset.cpp#L36 仅有前置声明、并未实现在本文件内且GDALRegister_USGSDEM()只注册了pfnOpen与pfnIdentify、未注册pfnCreateCopy因此可以推断当前版本的该驱动以只读为主更新模式下打开会被拒绝Open()中通过ReportUpdateNotSupportedByDriver(USGSDEM)明确拒绝GA_Update见 usgsdemdataset.cpp#L940-L945。驱动能力清单根据文档中的supports_georeferencing与supports_virtualio标记并结合驱动注册代码usgsdemdataset.cpp#L972-L993USGSDEM驱动的能力可归纳如下能力项是否支持依据栅格数据读取✅ 支持注册时设置GDAL_DCAP_RASTER YES地理配准GetGeoTransform✅ 支持GetGeoTransform()返回基于四角点与间距计算的仿射变换空间参考识别GetSpatialRef✅ 支持GetSpatialRef()返回 UTM / 州平面 / 地理坐标系虚拟文件系统/vsimem/ 等✅ 支持注册时设置GDAL_DCAP_VIRTUALIO YES更新写回❌ 不支持GA_Update打开被拒绝仅注册只读回调创建 / 复制输出❌ 当前源码未实现未注册pfnCreateCopy/pfnCreate默认扩展名.demGDAL_DMD_EXTENSION dem驱动长名USGS Optional ASCII DEM (and CDED)GDAL_DMD_LONGNAMEPAM 辅助信息与概览✅ 支持继承GDALPamDataset打开时TryLoadXML()并初始化oOvManager驱动的 CMake 构建定义位于 frmts/usgsdem/CMakeLists.txt采用add_gdal_driver(... PLUGIN_CAPABLE ...)方式注册既可静态内置于 GDAL也可编译为独立插件其在全局注册表中的调用点为 frmts/gdalallregister.cpp#L659GDALRegister_USGSDEM()。USGS DEM 文件格式基础头部结构与字段偏移USGS ASCII DEM 是定宽 ASCII 文本格式文件头按字节偏移组织。仓库中的 frmts/usgsdem/CDED.notes 提供了一份完整的 CDED 字段映射表与源码中LoadFromFile()usgsdemdataset.cpp#L584使用的偏移完全对应。关键字段如下偏移字节长度字段含义说明040文件名由输出文件名生成4060数据生产方模板或空白10926西南角地理坐标生成1351处理代码Process Code如 8ANUDEM、ATopoGrid1446DEM 级别代码CDED50K 产品为 11506高程模式常规模式为 11566水平参考系统0地理坐标Geographic1UTM2州平面1626UTM/SP 带号地理坐标时为 0168360投影参数通常为 0.05286水平单位1英尺ft、3角秒arc seconds等5346垂直单位1英尺、2米546192四角坐标SW、NW、NE、SE每角 24 字节x、y 各 12 字节73824最小高程生成76224最大高程生成81636空间分辨率X/Y 间距生成Z 通常为 1.08526剖面行数ProfilesCDED 硬编码为 18586剖面列数Columns of Profiles即 DEM 的列数像素宽度8764数据来源日期模板或空白8882垂直基准模板或 1MSL8902水平基准datum模板或 4NAD838924版本 / 规格版本CDED50K 产品规格版本为 10208964空白void百分比生成9008边缘匹配标志模板或空白上述字段映射来自仓库内 CDED.notesField Map 一节其写入方向面向 CDED 生产但对理解任意 USGS DEM 文件头的偏移布局同样适用且与读取驱动的解析偏移完全吻合。源码中驱动对头的解析均通过VSIFSeekL()定位到上述固定偏移后读取例如偏移 156水平参考系统代码nCoordSystem偏移 162UTM 带号iUTMZoneusgsdemdataset.cpp#L653-L655偏移 528 / 534水平单位nGUnit与垂直单位nVUnitusgsdemdataset.cpp#L657-L659偏移 816 起连续三个 12 字节字段X 间距dxdelta、Y 间距dydelta、垂直分辨率fVResusgsdemdataset.cpp#L667-L672偏移 858剖面总数nProfiles即栅格宽度usgsdemdataset.cpp#L703-L704偏移 546 起 192 字节四个角点坐标SW、NW、NE、SE各 24 字节usgsdemdataset.cpp#L685-L691。数据类型Data Type与单位识别USGSDEMRasterBand的GetUnitType()依据垂直单位代码返回高程单位字符串垂直单位为英尺nVUnit 1时返回ft否则返回m见 usgsdemdataset.cpp#L661-L665 与 usgsdemdataset.cpp#L537-L542。与此同时源码根据单位与垂直分辨率决定输出波段的数据类型usgsdemdataset.cpp#L676-L680条件数据类型说明垂直单位为英尺nVUnit1GDT_Float32保留原始数据精度垂直分辨率 1.0fVRes1.0GDT_Float32亚米级/亚英尺级分辨率需小数表达其他米制且分辨率 ≥ 1.0GDT_Int16整数米高程可直接用短整型承载从源码结构看这一判定的设计意图是DEM 数据以米为单位时值通常存为 short若以英尺为单位则以 float 存储以保留原始数据精度见LoadFromFile()上方注释usgsdemdataset.cpp#L574-L583。nodata 与缺失数据7.5 分钟图幅的边缘空白文档特别指出7.5 分钟UTM 网格USGS DEM 图幅在边缘通常存在缺失数据区域这些区域会被正确标记为 nodata。源码中定义的 nodata 常量为constexpr int USGSDEM_NODATA -32767;见 usgsdemdataset.cpp#L34。USGSDEMRasterBand::GetNoDataValue()恒返回该值usgsdemdataset.cpp#L525-L532。在数据读取阶段IReadBlock()usgsdemdataset.cpp#L350缓冲区先整体初始化为 nodata随后逐剖面将有效高程填入// 将输出块预填充为 nodata GDALCopyWords(USGSDEM_NODATA, GDT_Int32, 0, pImage, GetRasterDataType(), ...); // 逐剖面读取时若某像元为 nodata则保持输出缓冲区中的 nodata 不变 else if (nElev USGSDEM_NODATA) /* leave in output buffer as nodata */; else fComputedElev (float)(nElev * poGDS-fVRes dfElevOffset);因此边缘缺失区域无需文件显式存储即可自然呈现为 nodata这与文档描述完全一致。高程值的换算与写入剖面数据块中每条剖面profile包含行号、列号、点数、列数、X 起点、Y 起点、高程偏移量、最小/最大 Z 等字段。实际高程在读取时按以下公式还原usgsdemdataset.cpp#L482-L484elevation raw_int_value * fVRes dfElevOffset即原始整数值 × 垂直分辨率 高程偏移。若为地理坐标系Y 起点还会从角秒换算为十进制度数除以 3600usgsdemdataset.cpp#L439-L440。写入GDT_Int16波段时还会做范围钳制-32768 ~ 32767防止溢出usgsdemdataset.cpp#L486-L493。坐标参考系统识别机制USGSDEMDataset::LoadFromFile()依据头文件字段构建OGRSpatialReferenceusgsdemdataset.cpp#L709-L793水平基准datum映射水平基准代码位于偏移 890新格式旧格式默认 NAD27代码基准源码处理1NAD27SetWellKnownGeogCS(NAD27)2WGS 72SetWellKnownGeogCS(WGS72)3WGS 84SetWellKnownGeogCS(WGS84)4NAD83SetWellKnownGeogCS(NAD83)-9未知不设置基准其他 / 旧格式默认 NAD27SetWellKnownGeogCS(NAD27)投影系统分支nCoordSystem偏移 156决定投影分支usgsdemdataset.cpp#L768-L7911 UTM若带号在 [-60, 60] 区间内调用SetUTM()南北半球由带号正负决定若水平单位为英尺nGUnit 1还会以 US survey foot 作为线性单位并重命名投影字符串2 State Plane州平面调用SetStatePlane(iUTMZone, bNAD83, ...)同样区分 NAD27/NAD83 与英尺/米其他含 -9999 未知走通用分支按角秒换算地理坐标。注意新格式头的偏移 876 还存放数据编译年份驱动会读取该字段虽不直接参与 SRS 构建但用于区分格式版本见 usgsdemdataset.cpp#L714-L729。栅格尺寸与地理配准计算投影坐标UTM / 州平面栅格列数取nProfiles行数按(extent_max.y - extent_min.y) / dydelta 1.5计算X 起点强制取第一条剖面的 X 起点且像素锚点会按像元间距取模对齐m_gt中yorig取北边界之上半个像元、yscale为负的dydeltausgsdemdataset.cpp#L801-L831地理坐标直接用四角点外扩半个像元后的范围并把角秒转换为十进制度数除以 3600xscale/yscale亦同步除以 3600usgsdemdataset.cpp#L835-L852。这一外扩半个像元 负行间距的处理保证像元中心落在网格点上从而得到与 USGS DEM 规范一致的亚像元配准精度。测试 autotest/gdrivers/usgsdem.py 中的check_gt即是对该配准结果如(-67.00041667, 0.00083333, 0.0, 50.000416667, 0.0, -0.00083333)的回归校验。格式识别与多版本变体兼容Identify()如何判定文件是 USGS DEM驱动通过Identify()usgsdemdataset.cpp#L890-L908做快速识别检查头部偏移处的定宽文本偏移 156 处必须是 0、 1、 2、 3或 -9999水平参考系统代码的合法取值偏移 150 处必须是 1或 4高程模式为常规或 Level-4 的 CDED 变体。两者同时满足即判定为 USGSDEM 文件。四种数据起始偏移header 布局变体由于历史上有多个版本的记录布局LoadFromFile()会尝试依次探测数据区起始偏移usgsdemdataset.cpp#L586-L651起始偏移对应格式探测方法1024新格式New Format偏移 864 处读取的行/列号不为 1或文件位置 ≥ 1024且偏移 1024 处读到i1、j1 或 0893未文档化格式Undocumented Format偏移 893 处读到i1、j1如 39109h1.dem918A 记录最新迭代偏移 918 处读到i1、j1如 fema06-140cm_2995441b.dem864旧格式Old Format上述均不满足时的默认如 4619old.dem此外源码还针对两类行尾变体做了兼容usgsdemdataset.cpp#L632-L651数据区采用1025 字节记录、以换行符LF结尾的文件nDataStartOffset会 1对应 OSGeo/gdal#5007 的 1025 字节记录问题数据区采用1026 字节记录、以 CRLF 结尾的文件nDataStartOffset会 2对应 #15050 的 CRLF 问题。读取剖面时若数据起始偏移为 1024驱动还会在每条剖面结束后对齐到下一个 1024 字节边界以跳过部分文件中声明的有效剖面之外的垃圾剖面值usgsdemdataset.cpp#L503-L514。测试对变体覆盖的验证autotest/gdrivers/usgsdem.py 针对上述变体提供了完整回归测试测试数据位于 autotest/gdrivers/data/usgsdem/022gdeme_truncated地理坐标系、NAD27验证check_gt114p01_0100_deme_truncated.dem地理坐标系39079G6_truncated.demWGS72 UTM 17 带39109h1_truncated.dem未文档化格式893 偏移4619old_truncated.dem旧格式864 偏移usgsdem_with_extra_values_at_end_of_profile.dem剖面末尾多余值对应 OSGeo/gdal#583usgsdem_with_spaces_after_byte_864.dem偏移 864 之后出现空格对应 Novato.dem 的 #4901 问题fema06-140cm_2995441b_truncated.dem918 字节头部NAD83 UTM 15 带record_1025_ending_with_linefeed.dem1025 字节 LF 记录crlf.dem1026 字节 CRLF 记录。这些测试共同印证了大多数流行变体均应支持的文档承诺。实际使用命令行与 Python 示例用 gdalinfo 查看文件信息gdalinfo 022gdeme_truncated输出中你会看到Driver: USGSDEM/USGS Optional ASCII DEM (and CDED)Size is nProfiles x nRows列数来自偏移 858行数由范围/间距推算Origin与Pixel Size构成的仿射地理配准Coordinate System is中体现的 NAD27 / UTM 等空间参考波段 1 的TypeInt16/Float32、NoData Value-32767、Unit Typem/ft。用 gdal_translate 转换为其他格式由于驱动以只读为主常见的落地方式是将其作为读取源转换到 GeoTIFF 等可写格式gdal_translate -of GTiff -co COMPRESSDEFLATE -co TILEDYES \ in.dem out.tif若需要导出为 USGSDEM能力本文档提及曾由 Yukon 环境部门资助开发但如前所述当前仓库源码中未见对应的CreateCopy实现请以实际构建版本的驱动能力为准可通过gdalinfo --formats查看该驱动的 create 支持标记。Python 绑定示例from osgeo import gdal ds gdal.Open(in.dem) band ds.GetRasterBand(1) print(ds.RasterXSize, ds.RasterYSize) print(ds.GetGeoTransform()) # 仿射配准参数 print(ds.GetProjection()) # 空间参考 WKT print(band.GetNoDataValue()) # 通常为 -32767 print(band.GetUnitType()) # m 或 ft data band.ReadAsArray() # 读取整幅高程数组 ds None读取后即可对data做统计、裁剪、坡度/坡向计算等后续处理也可配合/vsimem/虚拟文件系统直接打开内存中的 DEM 字节流驱动已声明支持 VIRTUALIO。性能注意事项单一大瓦片的数据读取模型文档给出了两条明确的性能指引需要在实际使用中重视整个文件被表示为一个大 tile。USGSDEMRasterBand构造时nBlockXSize RasterXSize、nBlockYSize RasterYSizeusgsdemdataset.cpp#L334-L344即整个栅格是一个巨型块。首次读取像元时整个文件会被一次性消化ingest。IReadBlock()会 seek 到数据区并逐剖面读取全部内容同时文件内容读入使用一个最大 32 KB 的滚动缓冲区Bufferusgsdemdataset.cpp#L80-L103。由此带来两个工程建议保持 GDAL tile cache 足够大如果GDAL_CACHEMAX设置的 tile cache 偏小单一大块与缓存机制配合不佳可能引发缓存抖动cache thrashing。建议在批量读取前调高缓存上限。首像素延迟是预期的第一次ReadAsArray()或ReadBlock()会包含整个文件的解析开销属于该格式的正常行为并非程序卡死若只需范围/元数据用gdalinfo或仅调用GetGeoTransform()/GetProjection()不会触发数据区读取。源码实现结构与阅读指引驱动实现集中于 frmts/usgsdem/usgsdemdataset.cpp核心模块包括模块位置作用Identify()usgsdemdataset.cpp#L890头部特征快速识别Open()usgsdemdataset.cpp#L914打开流程识别 → 读取 → 建带 → PAM/概览初始化LoadFromFile()usgsdemdataset.cpp#L584核心解析格式探测、SRS 构建、尺寸与配准计算IReadBlock()usgsdemdataset.cpp#L350一次性读入全图并做高程换算与 nodata 填充GetUnitType()/GetNoDataValue()usgsdemdataset.cpp#L537、usgsdemdataset.cpp#L525单位与 nodata 元数据GDALRegister_USGSDEM()usgsdemdataset.cpp#L972驱动注册与能力声明配套资源文件头字段规范备注frmts/usgsdem/CDED.notes驱动构建配置frmts/usgsdem/CMakeLists.txt回归测试autotest/gdrivers/usgsdem.py 及测试数据目录 autotest/gdrivers/data/usgsdem/驱动全局注册frmts/gdalallregister.cpp#L659。局限性与注意事项小结只读驱动当前源码未实现写入/创建更新模式打开会被拒绝单一大块读取全图载入模型决定了内存与 I/O 的固定开销不适合对超大 DEM 做随机小块频繁读取单位与类型自动判定英尺或亚米分辨率数据会以Float32输出米制整米数据为Int16下游处理需注意类型差异格式变体依赖偏移探测对于极个别非标准头如 893/918 字节未文档化布局、1025/1026 字节记录驱动已做兼容但自定义实验性变体仍可能存在识别失败风险坐标基准默认值旧格式缺少 datum 字段时默认按 NAD27 处理与实际数据基准不一致时需在转换阶段通过-a_srs显式修正。综合来看USGSDEM驱动是 GDAL 读取 USGS 传统 ASCII DEM 与加拿大 CDED 数据的标准入口其固定偏移解析 多版本探测 自动配准的实现方式保证了历史数据与现代 GIS 工作流的无缝衔接是处理这类老式高程档案时可靠且开箱即用的解决方案。赞分享GIS遥感数据工程【免费下载链接】gdalGDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.项目地址https://gitcode.com/gh_mirrors/gd/gdal点击查看免费下载相关推荐GDAL OGR 之 Arc/Info E00 (ASCII) Coverage 读取驱动 AVCE00 完全指南GDAL OGR 之 Arc/Info E00 ASCII Coverage 读取驱动 AVCE00 完全指南 导读 Arc/Info E00ASCIICoGIS遥感数据工程GDAL XYZ 栅格驱动完全指南读写 ASCII Gridded XYZ 规则格网数据GDAL XYZ 栅格驱动完全指南读写 ASCII Gridded XYZ 规则格网数据 GDAL 内置的 XYZ 驱动 ASCII Gridded XYGIS遥感数据工程GDAL ISIS2 驱动完全指南USGS 行星影像 Cube v2 格式的读取、标签解析与地理参考处理GDAL ISIS2 驱动完全指南USGS 行星影像 Cube v2 格式的读取、标签解析与地理参考处理 导读 ISIS2 是 USGS美国地质调查局行星GIS遥感数据工程上一篇终极指南Android-BLE连接队列如何实现高效稳定的蓝牙通信下一篇LAVIS安装指南从PyPI到源码编译的完整流程创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考