简介江苏省地质岩性栅格数据面向地理、地质、遥感及环境科学等领域的研究者提供一套按地表出露岩性划分的栅格数据集。数据涵盖中性深成岩、中性火山岩、变质岩、基性深成岩、火山碎屑岩等十四种岩性类别并依据岩石成因机理归并为三大岩类对第四系松散堆积物亦有区分分类方法参考化学风化CO2消耗量评估相关文献制作有助于区域地质制图、风化过程建模及碳循环研究。数据包压缩后约246KB内含7个文件以TIF栅格为主体辅以坐标配准TFW、属性表DBF、元数据XML等支持文件可较好地在GIS软件中直接加载使用空间参考为WGS84精度250米。目前已有182人学习下载适合需要快速获取江苏省级中低分辨率岩性底图进行教学或预研究的用户。1. 250m岩性栅格和普通tif不一样文件组与分类才是核心在拿到江苏省岩性栅格数据的压缩包时多数人的第一反应是解压后直接双击那个lithology_江苏省.tif看到一片灰白影像就当成普通遥感图处理。真正用 GDAL 打开过之后会发现这张栅格每个像元存的是一个整数类型编号对应的是出露地表的岩石成因类别中性深成岩、中性火山岩、冰川、变质岩、基性深成岩、火山碎屑岩等十余种。数据精度 250m坐标系为 WGS84 经纬度制作逻辑来自化学风化过程中 CO2 消耗量评估的公开方法在很多英文文献里都能看到同源的分区思路。江苏面积不大250m 网格足够跑省级尺度的碳循环、水文地质和土壤侵蚀粗算但前提是把 tif 旁边那几个.tfw、.aux.xml、.vat.dbf文件当成一等公民而不是无用的附带物。下面从文件结构讲起一直做到 CO2 消耗量栅格输出。2. rar 解压后的岩性栅格数据tif、tfw、aux.xml、vat.dbf 各管一段2.1 栅格主文件之外的影子文件少了哪个都会出问题地理栅格数据在普通图片查看器里看起来就是一张图但在 GIS 软件里它由坐标、投影、属性表、像素尺寸共同决定。这个压缩包里真正承担空间定位和属性解释的是 tif 主文件旁边的四个伴随文件。它们的职责和缺失后果如下文件作用缺失后的典型问题lithology_江苏省.tif主栅格单波段整数类型像元值为岩性编码无此文件则整包不可用lithology_江苏省.tfwworld file记录像素尺寸、旋转参数、左上角坐标栅格失去地理参考叠加矢量时整体偏移lithology_江苏省.tif.aux.xml栅格统计信息与坐标系统描述GDAL 自动生成工具链需要重新扫描统计某些软件直接报错lithology_江苏省.tif.vat.dbfValue Attribute Table记录像元值与岩性名称的对应关系只看得到编码数字无法知道每个数字代表什么岩性lithology_江苏省.tif.vat.cpgdbf 文件的编码声明一般指向 UTF-8属性表中文乱码岩性名称变成符号我一般会在解压后先把这几个文件的修改时间比对一遍如果vat.dbf和vat.cppg的时间戳不一致说明文件在传输或二次压缩时可能被重新写过属性表的可靠性要先打个问号。.tfw只有六行文本每行一个数字分别对应 X 方向像素尺寸、旋转项、再一个旋转项、Y 方向像素尺寸、左上角 X、左上角 Y用手写配置文件甚至十六进制编辑器都可以核对这是栅格能被 GIS 正确放置的关键。2.2 用 gdalinfo 快速确认像元大小和投影状态处理这份数据前第一步永远是用gdalinfo检查它到底是不是 GeoTIFF以及坐标参考是否完整。命令行如下cd 江苏省岩性 gdalinfo lithology_江苏省.tif逻辑上这条命令只做只读操作不会改动原文件。参数含义很简单gdalinfo是 GDAL 自带的元数据查看工具对 tif 输出影像尺寸、波段数、像素尺寸、坐标参考、Origin左上角坐标和Pixel Size。对于这份江苏省岩性栅格重点是确认三个信息。第一坐标系名称是否为 WGS 84如果是输出里会出现GEOGCRS[WGS 84]或DATUM[WGS_1984]第二像素尺寸单位是度不要习惯性当成米WGS84 经纬度栅格在华东地区每度约对应 111km 到 103km跨度明显第三Band 1的NoData Value是否声明。如果 NoData 没有声明后续统计像元数会把值 0 或 255 误算进岩性面积里这是最容易被忽略的误差来源。看完gdalinfo输出后再动手做重分类才是安全的。2.3 用 ogrinfo 打开 vat.dbf先把编号和岩性对应关系读出来vat.dbf是 ArcGIS 体系里的栅格值属性表本质上是一个 dBASE 格式的关系表每行对应一个岩性类别。可以用 ogrinfo 直接读ogrinfo -al lithology_江苏省.tif.vat.dbf-al表示列出全部要素so如果加上则只显示概要信息。输出结果里通常会有Value、Count和岩性名称字段字段名可能是NAME、CLASS或中文名取决于原始生产时用的 ArcGIS 中文版还是英文版。这里要强调一个容易出错的地方.vat.dbf只存属性表不存储坐标系统信息。如果主 tif 重命名或移动.vat.dbf不会自动更新关联某些环境下打开 tif 后属性表是空的显示层里是一片未知编号。遇到这种情况不要把责任推到数据上应该检查 tif 文件名是否和 dbf 前缀一致。这份数据的前缀是lithology_江苏省.tif伴随文件必须保持完整前缀名改动任意一个字符都会导致属性表无法挂接。3. 14 种岩性值的编码和三大成因类重映射先看 vat.dbf再写代码3.1 岩性编码表与表达含义每个像元值是整数类型一般从 1 开始递增最高到 14正好对应压缩包说明中的十四种岩性。下面这张表列出的是一份常见排列方式实际使用时应该以你解压后得到的vat.dbf为准因为不同生产流程对编号顺序的定义会不一样。像元值岩性名称归入大类说明1冲积物沉积岩第四系松散堆积未固结成岩2湖相沉积沉积岩湖泊环境沉积粉砂黏土为主3海相沉积沉积岩沿海地区海相层4冰川堆积沉积岩冰川搬运堆积物5碳酸盐岩沉积岩石灰岩、白云岩类化学风化敏感度高6砂岩沉积岩碎屑岩类7泥岩和页岩沉积岩细粒碎屑岩8基性深成岩火成岩辉长岩、辉绿岩类9中性深成岩火成岩闪长岩、正长岩类10酸性深成岩火成岩花岗岩类11基性火山岩火成岩玄武岩类12中性火山岩火成岩安山岩类13酸性火山岩与火山碎屑岩火成岩流纹岩、凝灰岩类14变质岩变质岩片麻岩、片岩、石英岩类这张表的核心意义在于它不是按地层时代分类而是按岩石成因机理分类。原因在于化学风化驱动的 CO2 消耗速率对岩性矿物组成更敏感而不是对沉积时代更敏感。同一个地层单元里如果既有碳酸盐岩又有碎屑岩用年代图做参数赋值会产生明显偏差用成因图则可以直接对应不同风化反应路径。3.2 用 Python 做三大类重映射拿到 1 到 14 的编码后下一步是把 14 类压缩成三大类便于在水文模型中做概化。用 rasterio 读取栅格并执行映射import numpy as np import rasterio with rasterio.open(lithology_江苏省.tif) as src: data src.read(1).astype(np.uint8) nodata src.nodata profile src.profile # 按三大岩类重映射1沉积岩2火成岩3变质岩 reclass_map { 1: 1, 2: 1, 3: 1, 4: 1, 5: 1, 6: 1, 7: 1, 8: 2, 9: 2, 10: 2, 11: 2, 12: 2, 13: 2, 14: 3 } reclassed np.zeros_like(data, dtypenp.uint8) for old, new in reclass_map.items(): reclassed[data old] new # 无值区域保持 0 if nodata is not None: reclassed[data nodata] 0 profile.update(count1, dtypeuint8, nodata0) with rasterio.open(lithology_jiangsu_3class.tif, w, **profile) as dst: dst.write(reclassed, 1)这段代码的逻辑是先读入原始单波段栅格然后构造一个字典把 1 到 7 映射到沉积岩8 到 13 映射到火成岩14 映射到变质岩。for old, new in reclass_map.items()逐条处理像元值比用np.where嵌套多层条件更容易阅读和维护。profile.update是为了复用原始栅格的变换参数和投影信息同时保证输出文件与输入文件在空间位置上完全重叠这一点在后续做面积统计时非常重要。一个常见误区是直接把编号当成连续数值做回归比如把 2 和 14 取平均当作 8这在地质语义上没有意义。岩性数据是典型的名义分类变量如果要作为回归模型的输入特征正确的做法是转为独热编码而不是保留原始数值。3.3 NoData 的处理优先级很多人在这一步直接使用np.bincount(data)统计像元数量但 bincount 会把 0 值也数进去。对这种岩性栅格0 通常不代表真实岩性要么是背景要么是掩膜外区域。正确处理方式valid data ! 0 counts np.bincount(data[valid], minlength15)这里的valid布尔数组用于过滤无效像元minlength15保证输出数组长度覆盖 0 到 14 全部取值后续counts[1]到counts[14]与 14 个岩性一一对应。如果原始文件在 tif 头信息里定义了 NoData比如 255 或 0那么以src.nodata为准跟输出文件写入时设置的nodata0保持一致否则下一次重新打开输出文件时坐标参考与无效值的解释会对不上。4. 岩性栅格上的化学风化 CO2 消耗量评估面积统计与系数加权4.1 为什么这张栅格适合做 CO2 消耗量评估学术意义上地表岩石与大气 CO2 和水发生化学风化反应时硅酸盐岩与碳酸盐岩的消耗路径不同。碳酸盐岩溶解速度快对短时间尺度碳循环有明显响应玄武岩、辉长岩这类基性火成岩含有丰富易风化的钙镁硅酸盐矿物化学风化也可以成为净碳汇花岗岩等高硅酸性岩风化速率低消耗量相应偏小。所以做江苏省尺度的 CO2 消耗量估算核心不是拿整张影像求平均值而是先获得每一类岩性的面积权重再乘以对应的风化速率参考系数。这份数据在制作时已经按暴露地表的岩性划分并且叠加了化学风化 CO2 消耗量的评估逻辑相当于帮我们省掉了野外填图的步骤。但需要明确一点栅格只能给出空间分布不能直接给出消耗通量你必须把岩性类型映射到合理系数。4.2 按岩性统计面积并计算加权消耗量面积计算放在地理坐标栅格上常见做法有两种一是先把栅格重投影到等面积投影再按像元个数乘以固定像元面积二是在经纬度栅格上直接统计像元数后用像元中心纬度做每行平均距离校正。第二种方法写起来麻烦我建议直接把栅格重投影到 CGCS2000 分带投影后统计。下面是统计重投影后的像元数和面积并且计算 CO2 消耗量相对指标的示例import numpy as np import rasterio with rasterio.open(lithology_jiangsu_3class.tif) as src: class_data src.read(1) profile src.profile nodata src.nodata # CGCS2000 分带投影后像元尺寸为 250m x 250m pixel_area 250 * 250 / 1_000_000 # 0.0625 km2 valid class_data ! nodata counts np.bincount(class_data[valid], minlength15) # 各岩性的相对消耗参考系数数值仅代表相对高低不构成实测通量 co2_coef np.array([ 0.0, # 0 背景 0.2, # 1 冲积物 0.2, # 2 湖相沉积 0.3, # 3 海相沉积 0.2, # 4 冰川堆积 2.0, # 5 碳酸盐岩 0.6, # 6 砂岩 0.7, # 7 泥岩和页岩 1.5, # 8 基性深成岩 1.1, # 9 中性深成岩 0.5, # 10 酸性深成岩 1.6, # 11 基性火山岩 1.0, # 12 中性火山岩 0.6, # 13 酸性火山岩与火山碎屑岩 0.4 # 14 变质岩 ], dtypenp.float32) area_km2 counts[1:15] * pixel_area co2_relative area_km2 * co2_coef[1:15] for i, name in enumerate( [冲积物, 湖相沉积, 海相沉积, 冰川堆积, 碳酸盐岩, 砂岩, 泥页岩, 基性深成岩, 中性深成岩, 酸性深成岩, 基性火山岩, 中性火山岩, 酸性火山岩, 变质岩] ): print(f{name}: 面积 {area_km2[i]:.0f} km2, 相对消耗指标 {co2_relative[i]:.0f})参数说明pixel_area是固定像元面积在分带投影后每像元是 250m 乘 250m 的正方形co2_coef数组下标对应栅格类别编号长度 15 是为了让下标 1 到 14 直接对齐计算公式是面积乘以相对系数最后得到的是“相对消耗指标”不是真实吨数。如果将来要换算真实通量需要从文献中找到江苏省当地气候校正因子用湿润指数和径流深进一步修正。4.3 输出逐像元的消耗量栅格面积加权只能得到总量如果想看 CO2 消耗量的空间分布需要把系数直接映射到每个像元上。操作上不复杂只需把刚才的系数数组做成索引映射co2_grid co2_coef[class_data] * valid co2_grid np.where(valid, co2_grid, 0).astype(np.float32) profile.update(dtypefloat32, count1, nodata0) with rasterio.open(lithology_jiangsu_co2_relative.tif, w, **profile) as dst: dst.write(co2_grid, 1)这里的核心技巧是co2_coef[class_data]numpy 会把栅格中的类别值当作数组的下标一次性生成与输入栅格尺寸相同的浮点数组速度快且避免循环。valid布尔数组保证背景像元不被赋予系数值np.where再把无效区域置 0。输出文件继续沿用原 profile 的坐标变换保证结果与江苏省岩性原始栅格像元完全对齐后面在 QGIS 里做专题图裁切或者与流域矢量叠加时不会出现半个像元的偏移。打开后如果看到苏南沿江一带碳酸盐岩分布区有明显高值说明系数赋值逻辑基本正确。5. 出图前的三次检查重投影、vat.dbf 编码与栅格对齐5.1 检查一WGS84 经纬度栅格要不要转投影WGS84 是地理坐标系单位是度直接做距离和面积计算会有系统性误差。江苏地处中纬度250m 分辨率在纬度方向约 0.0025 度经度方向同样的度数对应距离明显更短。如果只是出图WGS84 没问题如果要算面积必须转到 CGCS2000 分带投影。常见做法是执行重投影gdalwarp -t_srs EPSG:4549 -tr 250 250 -r near \ lithology_jiangsu_3class.tif lithology_jiangsu_3class_cgcs2000.tif-tr 250 250定义输出分辨率-r near表示最邻近法重采样这样能保留原始整数编码不变避免双线性插值生成小数岩性值。执行后要再看一次gdalinfo确认 Pixel Size 已经是 250250。5.2 检查二vat.dbf 文件名和中文编码属性表不显示时先核对伴随文件名与主 tif 是否完全一致。如果重命名了 tif需要同步把.vat.dbf和.vat.cpg改掉。中文乱码大多数情况是 dbf 是 UTF-8 而软件默认按 GBK 读取用文本编辑器打开.vat.cpg确认内容是否为UTF-8不是就手动修正后重新加载。5.3 检查三与其他国家级栅格对齐栅格网当江苏省岩性栅格要和全国城市形态栅格数据集或其他 250m 格网叠加分析时两个栅格的投影、原点、像元对齐方式可能不同。建议统一用gdalwarp把岩性栅格变换到目标栅格的范围和分辨率再用掩模提取江苏省边界这样可以消除因边界浮点计算导致的相邻像元错位。最后把重投影结果和遥感影像叠加检查长江河道与图层边缘是否平行如果出现相同距离的固定偏移基本就是忽略了.tfw坐标参考应该回头检查投影定义。本文还有配套的精品资源点击获取