简介江苏省地质岩性栅格数据包面向地理信息系统、地质与生态环境研究从业者提供按地表出露岩性精细划分的十四类参数涵盖中性深成岩、中性火山岩、基性深成岩、火山碎屑岩、变质岩、冰川沉积等类别同时依据化学风化引起的CO2消耗量评估思路进行岩性分类便于开展区域地质分析、土壤侵蚀与碳循环模拟。数据采用WGS84坐标系空间精度为250米适合中尺度制图与模型输入。压缩包共含7个文件以GeoTIFF栅格文件为核心配套tfw坐标参考、xml元数据、dbf属性表及cpg编码说明另有使用说明与缩略图结构清晰解压即可加载。目前已有182人学习下载适合需要快速获取江苏省岩性底图的研究者或GIS学习者可直接用于制图、空间统计或作为模型输入参数省去自行矢量化与重分类的繁琐步骤。1. 江苏省地质岩性栅格数据一份压缩包里装着的分类图项目组拿到省域地质数据是常事但真正把一份“江苏省地质岩性栅格数据.rar”处理好让它能参与坡度分析、地类统计和工程选址往往比想象中多花两三天。岩性栅格不是一张普通图片它每个像素的取值是地质分类代码比如“1”代表第四系松散沉积物“3”代表灰岩数值本身没有连续意义。你无法对它做插值也不该用双线性重采样去缩放它否则会出现“不存在的岩性”。这份数据处理的核心是理解它作为分类栅格的语义、坐标系细节和编码规则然后用GDAL、rasterio这类工具完成解压、检视、重投影、重分类和统计分析。本文按拿到压缩包后的实际操作顺序展开读完后你至少能独立完成从.rar到可用分析图层的完整链路。2. 解压前端与识别判断压缩包完整性和栅格格式2.1 先验包rar完整性检查和加密判断拿到“江苏省地质岩性栅格数据.rar”这类命名规范的文件不要急着双击解压。省域地质数据通常体积在数百MB到几个GB如果是从同事、网盘或数据共享平台辗转拷贝压缩包损坏的概率远高于普通软件安装包。我一般先做完整性校验而不是让解压器中途报错打断。Linux/macOS环境下用unrar或7z命令行Windows环境用WinRAR或7-Zip图形界面但校验逻辑相同# 检查压缩包完整性不执行解压 unrar t Jiangsu_Lithology.rar # 或使用7z 7z t Jiangsu_Lithology.rarunrar t中的t是test模式只读取压缩包并逐文件校验CRC不落盘。如果输出中每一行都是OK或All OK说明文件完整。若出现CRC Failed或Unexpected end of archive说明包已损坏继续解压会得到不完整的GeoTIFF加载时可能出现全黑或异常值。此外用unrar l Jiangsu_Lithology.rar可以列出包内文件清单而不解压这样能在解压前确认里面是不是GeoTIFF、IMG还是ASCII Grid以及有没有配套的说明文档.docx/.pdf或图例文件。需要注意加密问题有些地质数据分发时会加密压缩包或设置只读密码。unrar l的输出里文件名后有星号*表示文件有密码保护。这时需要向数据提供方索要密码不要尝试暴力破解——一方面是时间成本不划算另一方面省域地质数据往往涉及保密审查要求流程合规比技术手段更重要。2.2 解压后的栅格格式识别GeoTIFF、IMG还是ASCII Grid解压完成得到一批文件后你要先搞清楚栅格数据的封装格式。江苏地质岩性数据常见三种形态GeoTIFF.tif/.tiff最主流包含地理参考信息GDAL原生支持推荐优先采用ERDAS IMAGINE.img老项目遗留格式GDAL同样支持ESRI ASCII Grid.asc文本格式可直接用文本编辑器打开看值但文件巨大不要用文件扩展名作为唯一判断依据尤其当系统隐藏了扩展名时。用Linux的file命令识别或者直接让GDAL去读# file命令查看文件真实类型 file Jiangsu_Lithology.tif # 直接让GDAL读取并打印驱动信息 gdalinfo Jiangsu_Lithology.tif | head -20gdalinfo这个命令在本文会反复出现它是GDAL的核心工具专门输出栅格数据元数据。第一个参数是文件路径。| head -20是Unix管道只取前20行输出避免字段太多刷屏。Windows的cmd用户直接用gdalinfo Jiangsu_Lithology.tif输出会打在屏幕上可用 info.txt重定向到文件再打开。2.3 元数据体检从gdalinfo输出读懂数据底细拿到一份岩性栅格最忌直接丢进ArcGIS/QGIS看渲染结果因为默认渲染很可能把全图显示成一片黑或一片白。先用gdalinfo把底细摸清gdalinfo Jiangsu_Lithology.tif关键看以下几行输出字段含义判定要点Driver: GTiff/GeoTIFF格式驱动GTiff就是GeoTIFFSize is 43620, 28660栅格宽高像素数省域1比20万数据常见为几万乘几万像素Origin (116.3, 35.2)左上角坐标坐标值带小数是地理坐标系带6到7位整数是投影坐标系Pixel Size (0.0001, -0.0001)像元尺寸单位是度则分辨率约11米单位是米则直接读NoData Value -9999无效值标记统计时需剔除GEOGCRS[CGCS2000]坐标系定义江苏数据最常见的国产坐标系举个实际输出片段Driver: GTiff/GeoTIFF Files: Jiangsu_Lithology.tif Size is 43620, 28660 Coordinate System is: GEOGCRS[China Geodetic Coordinate System 2000, DATUM[China 2000, ELLIPSOID[CGCS2000,6378137,298.257222101]], PRIMEM[Greenwich,0], ORIGIN[116.3,35.2]] Pixel Size (0.0001,-0.0001) Image Structure Metadata: COMPRESSIONDEFLATE INTERLEAVEBAND NoData Value-9999看到Pixel Size (0.0001, -0.0001)时要注意这是以度为单位的像元大小换算成米约等于11米在江苏约北纬32度位置1经度≈94公里0.0001度≈9.4米1纬度≈111公里0.0001度≈11.1米。如果你之后要和30米分辨率的DEM做叠置分析这种基础分辨率远超分析需求4亿像素覆盖全省单波段内存就要1.6GB建议先降采样再参与计算。3. 坐标系判定与重投影分类栅格必须用最近邻3.1 江苏地质数据里的坐标系CGCS2000是主角江苏全省范围横跨约东经116.3度到121.9度按高斯-克吕格投影分带位于3度带的第38带到第40带中央经线分别为114度、117度和120度。在实际数据中江苏地质岩性栅格的坐标系常见三种形态地理坐标系CGCS2000EPSG:4490以经纬度表示位置单位是度投影坐标系CGCS2000 / 3-degree Gauss-Kruger CM 120EEPSG:4547江苏中东部使用单位是米西安1980或北京1954坐标系老图件转栅格时的遗留坐标值和CGCS2000有几十到上百米的偏移判断方法就是看gdalinfo输出的GEOGCRS或PROJCRS段。GEOGCRS表示地理坐标系没有投影PROJCRS表示投影坐标系里面会有PROJECTION[Gauss_Kruger]和PARAMETER[central_meridian,120]这类描述。如果输出只有LOCAL_CS甚至Unknown说明坐标参考信息在栅格转出时丢失了这需要我们按第3.3节手动指定。3.2 投影与地理坐标系混用会怎样岩性栅格最常见的使用场景是和行政区划矢量、DEM、土地利用数据叠加。一个反复踩坑的情况是岩性栅格是经纬度GEOGCRS而你在ArcGIS里加载的县界矢量是高斯-克吕格投影PROJCRS软件会做动态投影看起来也能叠上但任何涉及缓冲区、面积计算和坡度提取的操作都会得到错误结果。ArcGIS的“动态投影”是为显示服务的不是为分析服务的。因此拿到数据的第一步操作是让岩性栅格和分析项目中其他数据的坐标系完全一致。我通常以项目中的基础底图为基准统一转成CGCS2000 / 3-degree Gauss-Kruger投影EPSG:4547覆盖江苏中东部但江苏东西跨度大最稳妥的方案是转成Albers等积投影或者直接看底图的坐标系并保持一致。3.3 gdalwarp重投影近邻采样是唯一正确选项重投影的命令如下gdalwarp -t_srs EPSG:4547 -r near -tr 30 30 -co COMPRESSDEFLATE \ -overwrite Jiangsu_Lithology.tif Jiangsu_Lithology_4547_30m.tif参数拆解-t_srs EPSG:4547目标坐标系EPSG:4547对应CGCS2000 / 3-degree Gauss-Kruger CM 120E。如果整个项目最终统一到其他坐标系改成对应EPSG码即可-r near重采样方法指定为最近邻nearest neighbor。这是分类栅格唯一正确的采样方式。双线性bilinear或三次卷积cubic会计算周围像素的加权平均得到89.7、134.5这类整数值。对地质岩性而言89.7这个值不指代任何岩性整层数据被污染-tr 30 30目标像元大小30米乘30米。源数据约11米分辨率降到30米可以减轻后续计算的IO压力。若不指定gdalwarp默认保持原分辨率省域数据会非常大-co COMPRESSDEFLATE输出GeoTIFF时做无损压缩。岩性栅格有大量相邻同类像素DEFLATE压缩率通常很可观11米分辨率的全省数据压缩后可能只剩原来的三分之一执行后再次用gdalinfo检查输出。如果看到Pixel Size (30,-30)和PROJCRS开头说明重投影成功。3.4 坐标信息丢失时的手动补救有些历史数据的GeoTIFF里没有坐标系字段gdalinfo显示Coordinate System is: Unknown。这时先别急着用-a_srs硬赋坐标——你需要先确认数据真实的坐标系是什么。判断依据包括文件附带的说明文档、同批次矢量数据的坐标系、或者用QGIS加载后和已知地理底图对照。如果确认是CGCS2000地理坐标系但GeoTIFF头部信息丢失使用gdal_translate赋坐标gdal_translate -a_srs EPSG:4490 -of GTiff \ source_nocrs.tif source_crs_set.tif-a_srs是assign SRS直接向输出文件写入坐标系而不重采样。注意这个操作不会改变像素位置只修改元数据所以只适用于“坐标数值正确但缺少坐标系标识”的场景。如果数据本身是投影坐标值单位是米你却赋成了地理坐标系单位是度整个图层的位置会发生灾难性偏移务必先用QGIS对照检查。4. 岩性编码解读与值域清理分类栅格不能当数值栅格用4.1 地质图GB编码与栅格值的关系岩性栅格的每个像元值对应的是一个地质分类编码。国内省域地质图件在栅格化时通常遵循两类编码体系地质矿产行业标准编码如侵入岩按年代岩性组合编号火山岩按喷发旋回和岩性细分沉积岩按地层年代和岩石类型区分项目自定义简编码如1第四系冲积层2白垩系砂岩3侏罗系火山岩4燕山期花岗岩——这种在工程应用中最常见因为用户不需要理解地质学的复杂分类层级拿到数据后第一步不是统计最大最小值而是先找编码说明。压缩包内如果有图例.doc或说明.pdf优先阅读。如果没有打开一张分辨率较低的预览图看色彩对应的类别注释。最后一个办法是抽取几个典型像元值结合江苏省地质图公开资料反推——比如宁镇山脉一带读到的值是灰岩类编码苏北平原读到的值是第四系松散沉积物编码基本可以确认编码体系对应的是沉积岩、侵入岩、变质岩大类下的细分。4.2 用gdalinfo -hist看分布用Python做唯一值统计分类栅格的核心统计不是均值方差而是每个编码出现的频次。gdalinfo -hist可以输出直方图但对分类栅格来说直方图分箱是连续的不适合精确统计。更实用的方式是直接用rasterio和numpy读入统计唯一值import numpy as np import rasterio from collections import Counter with rasterio.open(Jiangsu_Lithology_4547_30m.tif) as src: # 读取第一波段 data src.read(1) # 提取NoData掩膜假设NoData-9999 nodata src.nodata valid_mask data ! nodata # 统计所有有效像元的编码分布 valid_values data[valid_mask] counter Counter(valid_values.tolist()) # 打印前20个主要岩性编码及面积30m×30m900㎡/像素 print(f有效像元总数: {len(valid_values)}) for code, count in counter.most_common(20): area_km2 count * 900 / 1e6 print(f编码 {code}: {count} 像元, 约 {area_km2:.1f} 平方千米)这段代码做了三层处理。第一层是只读第一波段岩性栅格是单波段分类数据不需要读多波段。第二层用src.nodata读取元数据中的无效值标识配合布尔掩膜排除NoData对统计的干扰。第三层用Counter统计每个编码出现的次数并乘像元面积换算成平方公里。代码中的900来自30米乘30米的像元面积若你用了11米分辨率数据应改为121。输出示例编码 4: 1862354 像元, 约 167.6 平方千米 编码 1: 12686822 像元, 约 1141.8 平方千米 编码 2: 8353607 像元, 约 751.8 平方千米这种统计的价值在于快速判断数据的完整性。如果江苏全省面积约10.72万平方公里但各编码面积加起来仅有6万平方公里说明数据有大量区域是NoData或白值值为0需要回看原始数据是否只覆盖了部分区域或者压缩包内存在分幅文件而你没有合并。4.3 必须处理的NoData与白值陷阱岩性栅格最常见的问题是不设置NoData而是用值为0或者255表示背景区域。gdalinfo输出中如果没有NoData Value字段或者值为none就要警惕了。用以下方式检查背景值# 查看栅格值域范围 python -c import rasterio with rasterio.open(Jiangsu_Lithology_4547_30m.tif) as src: print(src.read(1).min(), src.read(1).max()) 如果最小值是0而你的编码从1开始0就是背景值。正确处理方式是重定义NoData而不是删掉这些像元——删除会改变栅格矩阵尺寸之后和DEM叠置时会错位。重定义用gdal_translategdal_translate -a_nodata 0 -of GTiff \ Jiangsu_Lithology_4547_30m.tif Jiangsu_Lithology_nodata0.tif-a_nodata 0把值0标记为NoData。这样后续所有涉及统计、分析和可视化的步骤都会自动排除背景区不需要每次写代码时手动做掩膜。5. 从整省到目标区裁剪、重分类与批量处理5.1 按地市边界裁剪cutline和crop_to_cutline的配合做工程应用时整省范围往往不直接可用你需要针对某一个市或某一个流域做分析。用gdalwarp配合矢量边界裁剪是标准做法gdalwarp -cutline Nanjing_City.shp -crop_to_cutline \ -t_srs EPSG:4547 -r near -tr 30 30 \ Jiangsu_Lithology_nodata0.tif Nanjing_Lithology.tif-cutline指定裁剪边界矢量文件-crop_to_cutline让输出栅格的范围收紧到边界本身的外接矩形而不是保留整省范围再把边界外裁掉——不写后者会得到一个仍然是全省大小的文件只是边界外变成NoData文件尺寸一点没小。注意裁剪矢量的坐标系gdalwarp会自动将cutline重投影到源栅格的坐标系但为了透明建议先确认Nanjing_City.shp的坐标系比并做好预处理。5.2 岩性重分类gdal_calc.py把细类归并为大类地质部门提供的编码往往细分到岩石类型比如灰岩分“石灰岩”“白云质灰岩”“鲕状灰岩”对工程建设来说通常只需要归并成三大类松散沉积物、坚硬岩、软弱岩。重分类有多种方式最灵活的是gdal_calc.pygdal_calc.py -A Nanjing_Lithology.tif \ --outfileNaNjing_Lithology_reclass.tif \ --calc((A1)*(A5))*10 ((A6)*(A12))*20 ((A13)*(A20))*30 \ --NoDataValue0 --typeByte这条命令按编码范围做了三段归并编码1到5的松散沉积物对应新值10编码6到12的沉积岩对应新值20编码13到20的岩浆岩对应新值30。--calc参数中的表达式基于numpy语法((A1)*(A5))返回布尔数组与数值相乘时True转为1、False转为0于是代码落在对应段位就得到新值否则为0。--NoDataValue0保证背景区继续被排除--typeByte把输出数据类型压缩为8位无符号整数节省存储空间。执行后必须验证在QGIS中加载新文件随机取点对照原始栅格的编码确认重分类映射没有错位。5.3 分幅栅格的合并策略buildvrt永远比直接merge快省域地质栅格数据常见分幅存储比如按1比20万图幅切成了几十个文件。此时不要直接对几十个GeoTIFF执行gdal_merge.py那个工具会真实读写所有像素合并过程耗时很长且占用大量临时磁盘。正确做法是先构建VRT虚拟栅格再按需转译# 先构建VRT不复制像素数据 gdalbuildvrt -resolution highest -r near \ -input_file_list list.txt Jiangsu_Lithology_all.vrt # 按需输出为单个GeoTIFF gdal_translate -co COMPRESSDEFLATE -co TILEDYES \ Jiangsu_Lithology_all.vrt Jiangsu_Lithology_merged.tif-resolution highest让VRT统一所有分幅图的分辨率以最高的为准避免一些图幅11米、另一些图幅30米导致后续处理错乱。-r near同样是分类栅格的生命线。VRT文件可能只有几十KB但查看和读取它时GDAL会自动拼接虚拟的底层文件后续gdalinfo、gdalwarp都可以直接操作VRT甚至不必物理合并。5.4 与DEM叠置岩性-坡度联合统计分析投入实际生产时岩性栅格很少单独使用。一个高频分析场景是“岩性坡度”对地质灾害风险的影响。把岩性栅格和坡度栅格叠置统计每类岩性所在区域的坡度分布需要先将二者对齐import numpy as np import rasterio from rasterio.warp import reproject, Resampling # 打开岩性栅格和坡度栅格 with rasterio.open(Nanjing_Lithology_reclass.tif) as litho_src, \ rasterio.open(Nanjing_Slope_30m.tif) as slope_src: # 读取岩性数据并获取其元数据目标网格 litho litho_src.read(1) litho_profile litho_src.profile.copy() # 将坡度数据重投影/重采样到岩性栅格的网格 slope_resampled np.zeros_like(litho, dtypenp.float32) reproject( sourcerasterio.band(slope_src, 1), destinationslope_resampled, src_transformslope_src.transform, src_crsslope_src.crs, dst_transformlitho_src.transform, dst_crslitho_src.crs, resamplingResampling.bilinear ) # 按岩性编码统计平均坡度 import pandas as pd records [] for code in np.unique(litho): if code 0: continue mask litho code if mask.sum() 1000: # 只统计面积足够的岩性 records.append((code, slope_resampled[mask].mean(), slope_resampled[mask].max())) df pd.DataFrame(records, columns[litho_code, mean_slope, max_slope]) print(df)这里注意一个分类栅格与连续栅格叠置时的规范问题岩性栅格用最近邻保持原值坡度栅格用双线性插值计算平滑坡度两个栅格即便分辨率相同重采样逻辑也必须分开对待。6. 高频排错与验证技巧栅格看起来对结果却不对地质岩性栅格数据处理的投入产出比往往卡在最后几处细节。这里给出三个我反复用到的高频排错手段。第一直方图尾部异常值的识别。用QGIS加载后查看直方图如果发现某个编码值占比异常地高比如编码127占了85%的像元先别急着判定数据合理。在空值未正确设置时GDAL有时会用127填充背景区。此时回到元数据检查NoData设置。如果源数据本身没有NoData元数据需要重新审视4.3节的白值定义策略背景区域通常占据省域数据的20%到40%不清理干净后续任何统计都会严重偏差。第二坐标偏移的快速验证。完成重投影和裁剪后用QGIS“识别”工具点击一个已知地质点位。举例南京栖霞山的栖霞寺附近应是二叠系栖霞组灰岩如果栅格值对应灰岩编码则配准正确。即使只有三五个控制点这种“地质底图人肉验证”也远比只看边界形状可靠因为边界在低分辨率下错位几十米是肉眼难辨的。第三栅格读取的IO瓶颈处理。处理整省30米分辨率的GeoTIFF像素量约1.2亿Python循环读值完全不可行。正确的效率做法是使用rasterio窗口读取src.read(1, windowWindow(...))分块处理或者用gdal_calc.py这类C语言底层实现完成逐像元运算。遇到“内存不足”报错不要加大内存改用--co TILEDYES输出瓦片化GeoTIFF并用numpy.memmap做阵列映射性能会提升一个量级。最后当你对重分类结果有疑问时用gdalinfo -stats重新计算全图层直方图并和原始数据对比分类栅格正确处理后的统计特征应该是“各编码面积在归并方向上守恒”否则追查重分类表达式里漏掉的范围分支。本文还有配套的精品资源点击获取 地质岩性栅格数据处理实战：从RAR解压到GIS分析全流程
简介江苏省地质岩性栅格数据包面向地理信息系统、地质与生态环境研究从业者提供按地表出露岩性精细划分的十四类参数涵盖中性深成岩、中性火山岩、基性深成岩、火山碎屑岩、变质岩、冰川沉积等类别同时依据化学风化引起的CO2消耗量评估思路进行岩性分类便于开展区域地质分析、土壤侵蚀与碳循环模拟。数据采用WGS84坐标系空间精度为250米适合中尺度制图与模型输入。压缩包共含7个文件以GeoTIFF栅格文件为核心配套tfw坐标参考、xml元数据、dbf属性表及cpg编码说明另有使用说明与缩略图结构清晰解压即可加载。目前已有182人学习下载适合需要快速获取江苏省岩性底图的研究者或GIS学习者可直接用于制图、空间统计或作为模型输入参数省去自行矢量化与重分类的繁琐步骤。1. 江苏省地质岩性栅格数据一份压缩包里装着的分类图项目组拿到省域地质数据是常事但真正把一份“江苏省地质岩性栅格数据.rar”处理好让它能参与坡度分析、地类统计和工程选址往往比想象中多花两三天。岩性栅格不是一张普通图片它每个像素的取值是地质分类代码比如“1”代表第四系松散沉积物“3”代表灰岩数值本身没有连续意义。你无法对它做插值也不该用双线性重采样去缩放它否则会出现“不存在的岩性”。这份数据处理的核心是理解它作为分类栅格的语义、坐标系细节和编码规则然后用GDAL、rasterio这类工具完成解压、检视、重投影、重分类和统计分析。本文按拿到压缩包后的实际操作顺序展开读完后你至少能独立完成从.rar到可用分析图层的完整链路。2. 解压前端与识别判断压缩包完整性和栅格格式2.1 先验包rar完整性检查和加密判断拿到“江苏省地质岩性栅格数据.rar”这类命名规范的文件不要急着双击解压。省域地质数据通常体积在数百MB到几个GB如果是从同事、网盘或数据共享平台辗转拷贝压缩包损坏的概率远高于普通软件安装包。我一般先做完整性校验而不是让解压器中途报错打断。Linux/macOS环境下用unrar或7z命令行Windows环境用WinRAR或7-Zip图形界面但校验逻辑相同# 检查压缩包完整性不执行解压 unrar t Jiangsu_Lithology.rar # 或使用7z 7z t Jiangsu_Lithology.rarunrar t中的t是test模式只读取压缩包并逐文件校验CRC不落盘。如果输出中每一行都是OK或All OK说明文件完整。若出现CRC Failed或Unexpected end of archive说明包已损坏继续解压会得到不完整的GeoTIFF加载时可能出现全黑或异常值。此外用unrar l Jiangsu_Lithology.rar可以列出包内文件清单而不解压这样能在解压前确认里面是不是GeoTIFF、IMG还是ASCII Grid以及有没有配套的说明文档.docx/.pdf或图例文件。需要注意加密问题有些地质数据分发时会加密压缩包或设置只读密码。unrar l的输出里文件名后有星号*表示文件有密码保护。这时需要向数据提供方索要密码不要尝试暴力破解——一方面是时间成本不划算另一方面省域地质数据往往涉及保密审查要求流程合规比技术手段更重要。2.2 解压后的栅格格式识别GeoTIFF、IMG还是ASCII Grid解压完成得到一批文件后你要先搞清楚栅格数据的封装格式。江苏地质岩性数据常见三种形态GeoTIFF.tif/.tiff最主流包含地理参考信息GDAL原生支持推荐优先采用ERDAS IMAGINE.img老项目遗留格式GDAL同样支持ESRI ASCII Grid.asc文本格式可直接用文本编辑器打开看值但文件巨大不要用文件扩展名作为唯一判断依据尤其当系统隐藏了扩展名时。用Linux的file命令识别或者直接让GDAL去读# file命令查看文件真实类型 file Jiangsu_Lithology.tif # 直接让GDAL读取并打印驱动信息 gdalinfo Jiangsu_Lithology.tif | head -20gdalinfo这个命令在本文会反复出现它是GDAL的核心工具专门输出栅格数据元数据。第一个参数是文件路径。| head -20是Unix管道只取前20行输出避免字段太多刷屏。Windows的cmd用户直接用gdalinfo Jiangsu_Lithology.tif输出会打在屏幕上可用 info.txt重定向到文件再打开。2.3 元数据体检从gdalinfo输出读懂数据底细拿到一份岩性栅格最忌直接丢进ArcGIS/QGIS看渲染结果因为默认渲染很可能把全图显示成一片黑或一片白。先用gdalinfo把底细摸清gdalinfo Jiangsu_Lithology.tif关键看以下几行输出字段含义判定要点Driver: GTiff/GeoTIFF格式驱动GTiff就是GeoTIFFSize is 43620, 28660栅格宽高像素数省域1比20万数据常见为几万乘几万像素Origin (116.3, 35.2)左上角坐标坐标值带小数是地理坐标系带6到7位整数是投影坐标系Pixel Size (0.0001, -0.0001)像元尺寸单位是度则分辨率约11米单位是米则直接读NoData Value -9999无效值标记统计时需剔除GEOGCRS[CGCS2000]坐标系定义江苏数据最常见的国产坐标系举个实际输出片段Driver: GTiff/GeoTIFF Files: Jiangsu_Lithology.tif Size is 43620, 28660 Coordinate System is: GEOGCRS[China Geodetic Coordinate System 2000, DATUM[China 2000, ELLIPSOID[CGCS2000,6378137,298.257222101]], PRIMEM[Greenwich,0], ORIGIN[116.3,35.2]] Pixel Size (0.0001,-0.0001) Image Structure Metadata: COMPRESSIONDEFLATE INTERLEAVEBAND NoData Value-9999看到Pixel Size (0.0001, -0.0001)时要注意这是以度为单位的像元大小换算成米约等于11米在江苏约北纬32度位置1经度≈94公里0.0001度≈9.4米1纬度≈111公里0.0001度≈11.1米。如果你之后要和30米分辨率的DEM做叠置分析这种基础分辨率远超分析需求4亿像素覆盖全省单波段内存就要1.6GB建议先降采样再参与计算。3. 坐标系判定与重投影分类栅格必须用最近邻3.1 江苏地质数据里的坐标系CGCS2000是主角江苏全省范围横跨约东经116.3度到121.9度按高斯-克吕格投影分带位于3度带的第38带到第40带中央经线分别为114度、117度和120度。在实际数据中江苏地质岩性栅格的坐标系常见三种形态地理坐标系CGCS2000EPSG:4490以经纬度表示位置单位是度投影坐标系CGCS2000 / 3-degree Gauss-Kruger CM 120EEPSG:4547江苏中东部使用单位是米西安1980或北京1954坐标系老图件转栅格时的遗留坐标值和CGCS2000有几十到上百米的偏移判断方法就是看gdalinfo输出的GEOGCRS或PROJCRS段。GEOGCRS表示地理坐标系没有投影PROJCRS表示投影坐标系里面会有PROJECTION[Gauss_Kruger]和PARAMETER[central_meridian,120]这类描述。如果输出只有LOCAL_CS甚至Unknown说明坐标参考信息在栅格转出时丢失了这需要我们按第3.3节手动指定。3.2 投影与地理坐标系混用会怎样岩性栅格最常见的使用场景是和行政区划矢量、DEM、土地利用数据叠加。一个反复踩坑的情况是岩性栅格是经纬度GEOGCRS而你在ArcGIS里加载的县界矢量是高斯-克吕格投影PROJCRS软件会做动态投影看起来也能叠上但任何涉及缓冲区、面积计算和坡度提取的操作都会得到错误结果。ArcGIS的“动态投影”是为显示服务的不是为分析服务的。因此拿到数据的第一步操作是让岩性栅格和分析项目中其他数据的坐标系完全一致。我通常以项目中的基础底图为基准统一转成CGCS2000 / 3-degree Gauss-Kruger投影EPSG:4547覆盖江苏中东部但江苏东西跨度大最稳妥的方案是转成Albers等积投影或者直接看底图的坐标系并保持一致。3.3 gdalwarp重投影近邻采样是唯一正确选项重投影的命令如下gdalwarp -t_srs EPSG:4547 -r near -tr 30 30 -co COMPRESSDEFLATE \ -overwrite Jiangsu_Lithology.tif Jiangsu_Lithology_4547_30m.tif参数拆解-t_srs EPSG:4547目标坐标系EPSG:4547对应CGCS2000 / 3-degree Gauss-Kruger CM 120E。如果整个项目最终统一到其他坐标系改成对应EPSG码即可-r near重采样方法指定为最近邻nearest neighbor。这是分类栅格唯一正确的采样方式。双线性bilinear或三次卷积cubic会计算周围像素的加权平均得到89.7、134.5这类整数值。对地质岩性而言89.7这个值不指代任何岩性整层数据被污染-tr 30 30目标像元大小30米乘30米。源数据约11米分辨率降到30米可以减轻后续计算的IO压力。若不指定gdalwarp默认保持原分辨率省域数据会非常大-co COMPRESSDEFLATE输出GeoTIFF时做无损压缩。岩性栅格有大量相邻同类像素DEFLATE压缩率通常很可观11米分辨率的全省数据压缩后可能只剩原来的三分之一执行后再次用gdalinfo检查输出。如果看到Pixel Size (30,-30)和PROJCRS开头说明重投影成功。3.4 坐标信息丢失时的手动补救有些历史数据的GeoTIFF里没有坐标系字段gdalinfo显示Coordinate System is: Unknown。这时先别急着用-a_srs硬赋坐标——你需要先确认数据真实的坐标系是什么。判断依据包括文件附带的说明文档、同批次矢量数据的坐标系、或者用QGIS加载后和已知地理底图对照。如果确认是CGCS2000地理坐标系但GeoTIFF头部信息丢失使用gdal_translate赋坐标gdal_translate -a_srs EPSG:4490 -of GTiff \ source_nocrs.tif source_crs_set.tif-a_srs是assign SRS直接向输出文件写入坐标系而不重采样。注意这个操作不会改变像素位置只修改元数据所以只适用于“坐标数值正确但缺少坐标系标识”的场景。如果数据本身是投影坐标值单位是米你却赋成了地理坐标系单位是度整个图层的位置会发生灾难性偏移务必先用QGIS对照检查。4. 岩性编码解读与值域清理分类栅格不能当数值栅格用4.1 地质图GB编码与栅格值的关系岩性栅格的每个像元值对应的是一个地质分类编码。国内省域地质图件在栅格化时通常遵循两类编码体系地质矿产行业标准编码如侵入岩按年代岩性组合编号火山岩按喷发旋回和岩性细分沉积岩按地层年代和岩石类型区分项目自定义简编码如1第四系冲积层2白垩系砂岩3侏罗系火山岩4燕山期花岗岩——这种在工程应用中最常见因为用户不需要理解地质学的复杂分类层级拿到数据后第一步不是统计最大最小值而是先找编码说明。压缩包内如果有图例.doc或说明.pdf优先阅读。如果没有打开一张分辨率较低的预览图看色彩对应的类别注释。最后一个办法是抽取几个典型像元值结合江苏省地质图公开资料反推——比如宁镇山脉一带读到的值是灰岩类编码苏北平原读到的值是第四系松散沉积物编码基本可以确认编码体系对应的是沉积岩、侵入岩、变质岩大类下的细分。4.2 用gdalinfo -hist看分布用Python做唯一值统计分类栅格的核心统计不是均值方差而是每个编码出现的频次。gdalinfo -hist可以输出直方图但对分类栅格来说直方图分箱是连续的不适合精确统计。更实用的方式是直接用rasterio和numpy读入统计唯一值import numpy as np import rasterio from collections import Counter with rasterio.open(Jiangsu_Lithology_4547_30m.tif) as src: # 读取第一波段 data src.read(1) # 提取NoData掩膜假设NoData-9999 nodata src.nodata valid_mask data ! nodata # 统计所有有效像元的编码分布 valid_values data[valid_mask] counter Counter(valid_values.tolist()) # 打印前20个主要岩性编码及面积30m×30m900㎡/像素 print(f有效像元总数: {len(valid_values)}) for code, count in counter.most_common(20): area_km2 count * 900 / 1e6 print(f编码 {code}: {count} 像元, 约 {area_km2:.1f} 平方千米)这段代码做了三层处理。第一层是只读第一波段岩性栅格是单波段分类数据不需要读多波段。第二层用src.nodata读取元数据中的无效值标识配合布尔掩膜排除NoData对统计的干扰。第三层用Counter统计每个编码出现的次数并乘像元面积换算成平方公里。代码中的900来自30米乘30米的像元面积若你用了11米分辨率数据应改为121。输出示例编码 4: 1862354 像元, 约 167.6 平方千米 编码 1: 12686822 像元, 约 1141.8 平方千米 编码 2: 8353607 像元, 约 751.8 平方千米这种统计的价值在于快速判断数据的完整性。如果江苏全省面积约10.72万平方公里但各编码面积加起来仅有6万平方公里说明数据有大量区域是NoData或白值值为0需要回看原始数据是否只覆盖了部分区域或者压缩包内存在分幅文件而你没有合并。4.3 必须处理的NoData与白值陷阱岩性栅格最常见的问题是不设置NoData而是用值为0或者255表示背景区域。gdalinfo输出中如果没有NoData Value字段或者值为none就要警惕了。用以下方式检查背景值# 查看栅格值域范围 python -c import rasterio with rasterio.open(Jiangsu_Lithology_4547_30m.tif) as src: print(src.read(1).min(), src.read(1).max()) 如果最小值是0而你的编码从1开始0就是背景值。正确处理方式是重定义NoData而不是删掉这些像元——删除会改变栅格矩阵尺寸之后和DEM叠置时会错位。重定义用gdal_translategdal_translate -a_nodata 0 -of GTiff \ Jiangsu_Lithology_4547_30m.tif Jiangsu_Lithology_nodata0.tif-a_nodata 0把值0标记为NoData。这样后续所有涉及统计、分析和可视化的步骤都会自动排除背景区不需要每次写代码时手动做掩膜。5. 从整省到目标区裁剪、重分类与批量处理5.1 按地市边界裁剪cutline和crop_to_cutline的配合做工程应用时整省范围往往不直接可用你需要针对某一个市或某一个流域做分析。用gdalwarp配合矢量边界裁剪是标准做法gdalwarp -cutline Nanjing_City.shp -crop_to_cutline \ -t_srs EPSG:4547 -r near -tr 30 30 \ Jiangsu_Lithology_nodata0.tif Nanjing_Lithology.tif-cutline指定裁剪边界矢量文件-crop_to_cutline让输出栅格的范围收紧到边界本身的外接矩形而不是保留整省范围再把边界外裁掉——不写后者会得到一个仍然是全省大小的文件只是边界外变成NoData文件尺寸一点没小。注意裁剪矢量的坐标系gdalwarp会自动将cutline重投影到源栅格的坐标系但为了透明建议先确认Nanjing_City.shp的坐标系比并做好预处理。5.2 岩性重分类gdal_calc.py把细类归并为大类地质部门提供的编码往往细分到岩石类型比如灰岩分“石灰岩”“白云质灰岩”“鲕状灰岩”对工程建设来说通常只需要归并成三大类松散沉积物、坚硬岩、软弱岩。重分类有多种方式最灵活的是gdal_calc.pygdal_calc.py -A Nanjing_Lithology.tif \ --outfileNaNjing_Lithology_reclass.tif \ --calc((A1)*(A5))*10 ((A6)*(A12))*20 ((A13)*(A20))*30 \ --NoDataValue0 --typeByte这条命令按编码范围做了三段归并编码1到5的松散沉积物对应新值10编码6到12的沉积岩对应新值20编码13到20的岩浆岩对应新值30。--calc参数中的表达式基于numpy语法((A1)*(A5))返回布尔数组与数值相乘时True转为1、False转为0于是代码落在对应段位就得到新值否则为0。--NoDataValue0保证背景区继续被排除--typeByte把输出数据类型压缩为8位无符号整数节省存储空间。执行后必须验证在QGIS中加载新文件随机取点对照原始栅格的编码确认重分类映射没有错位。5.3 分幅栅格的合并策略buildvrt永远比直接merge快省域地质栅格数据常见分幅存储比如按1比20万图幅切成了几十个文件。此时不要直接对几十个GeoTIFF执行gdal_merge.py那个工具会真实读写所有像素合并过程耗时很长且占用大量临时磁盘。正确做法是先构建VRT虚拟栅格再按需转译# 先构建VRT不复制像素数据 gdalbuildvrt -resolution highest -r near \ -input_file_list list.txt Jiangsu_Lithology_all.vrt # 按需输出为单个GeoTIFF gdal_translate -co COMPRESSDEFLATE -co TILEDYES \ Jiangsu_Lithology_all.vrt Jiangsu_Lithology_merged.tif-resolution highest让VRT统一所有分幅图的分辨率以最高的为准避免一些图幅11米、另一些图幅30米导致后续处理错乱。-r near同样是分类栅格的生命线。VRT文件可能只有几十KB但查看和读取它时GDAL会自动拼接虚拟的底层文件后续gdalinfo、gdalwarp都可以直接操作VRT甚至不必物理合并。5.4 与DEM叠置岩性-坡度联合统计分析投入实际生产时岩性栅格很少单独使用。一个高频分析场景是“岩性坡度”对地质灾害风险的影响。把岩性栅格和坡度栅格叠置统计每类岩性所在区域的坡度分布需要先将二者对齐import numpy as np import rasterio from rasterio.warp import reproject, Resampling # 打开岩性栅格和坡度栅格 with rasterio.open(Nanjing_Lithology_reclass.tif) as litho_src, \ rasterio.open(Nanjing_Slope_30m.tif) as slope_src: # 读取岩性数据并获取其元数据目标网格 litho litho_src.read(1) litho_profile litho_src.profile.copy() # 将坡度数据重投影/重采样到岩性栅格的网格 slope_resampled np.zeros_like(litho, dtypenp.float32) reproject( sourcerasterio.band(slope_src, 1), destinationslope_resampled, src_transformslope_src.transform, src_crsslope_src.crs, dst_transformlitho_src.transform, dst_crslitho_src.crs, resamplingResampling.bilinear ) # 按岩性编码统计平均坡度 import pandas as pd records [] for code in np.unique(litho): if code 0: continue mask litho code if mask.sum() 1000: # 只统计面积足够的岩性 records.append((code, slope_resampled[mask].mean(), slope_resampled[mask].max())) df pd.DataFrame(records, columns[litho_code, mean_slope, max_slope]) print(df)这里注意一个分类栅格与连续栅格叠置时的规范问题岩性栅格用最近邻保持原值坡度栅格用双线性插值计算平滑坡度两个栅格即便分辨率相同重采样逻辑也必须分开对待。6. 高频排错与验证技巧栅格看起来对结果却不对地质岩性栅格数据处理的投入产出比往往卡在最后几处细节。这里给出三个我反复用到的高频排错手段。第一直方图尾部异常值的识别。用QGIS加载后查看直方图如果发现某个编码值占比异常地高比如编码127占了85%的像元先别急着判定数据合理。在空值未正确设置时GDAL有时会用127填充背景区。此时回到元数据检查NoData设置。如果源数据本身没有NoData元数据需要重新审视4.3节的白值定义策略背景区域通常占据省域数据的20%到40%不清理干净后续任何统计都会严重偏差。第二坐标偏移的快速验证。完成重投影和裁剪后用QGIS“识别”工具点击一个已知地质点位。举例南京栖霞山的栖霞寺附近应是二叠系栖霞组灰岩如果栅格值对应灰岩编码则配准正确。即使只有三五个控制点这种“地质底图人肉验证”也远比只看边界形状可靠因为边界在低分辨率下错位几十米是肉眼难辨的。第三栅格读取的IO瓶颈处理。处理整省30米分辨率的GeoTIFF像素量约1.2亿Python循环读值完全不可行。正确的效率做法是使用rasterio窗口读取src.read(1, windowWindow(...))分块处理或者用gdal_calc.py这类C语言底层实现完成逐像元运算。遇到“内存不足”报错不要加大内存改用--co TILEDYES输出瓦片化GeoTIFF并用numpy.memmap做阵列映射性能会提升一个量级。最后当你对重分类结果有疑问时用gdalinfo -stats重新计算全图层直方图并和原始数据对比分类栅格正确处理后的统计特征应该是“各编码面积在归并方向上守恒”否则追查重分类表达式里漏掉的范围分支。本文还有配套的精品资源点击获取