★ FEATURED ARTICLE 简介这份资源提供广东省全域250米分辨率的地表出露岩性分布栅格数据采用WGS84地理坐标系共划分14类岩性类型涵盖中性深成岩、中性火山岩、冰川沉积物、各类变质岩、基性深成岩及火山碎屑岩等并按岩石成因归入火成岩、沉积岩、变质岩三大类其中“土”类指未固结的第四纪松散堆积物。分类依据为地表可见岩性不涉及地下深部岩层推断适合从事地质环境评价、水文模拟、土壤侵蚀分析及CO2化学风化消耗量估算等方向的研究人员与高年级学生使用。资源包共13个文件约831KB以tif主栅格为核心辅以dbf属性表、tfw坐标信息、cpg编码文件与aux.xml元数据另含png预览图、xls说明表、py脚本及使用说明压缩包便于直接加载与二次处理。目前已有17人学习下载可帮助读者快速获取全省岩性底图支撑空间分析与建模工作。1. 广东全省250米分辨率地表出露岩性分布栅格数据一份能直接进GIS的底图资源做华南地区水文建模、地质灾害易发性评价或者土壤类型推演的同行大概率都经历过同一个卡脖子环节手上有了DEM、降雨、NDVI唯独缺一张能跟其他栅格对齐的岩性底图。公开的全国地质图要么是1:250万的小比例尺矢量要么是分幅扫描件配准和矢量化能耗掉一整周。这份广东全省250米分辨率地表出露岩性分布栅格数据解决的正是这个痛点——它把全省地表出露的岩性类型直接栅格化到250米格网坐标系是WGS84地理坐标拿到手就能和SRTM、MODIS、Landsat系列产品做像元级叠加。适合做区域尺度空间分析的研究生、做地质灾害普查的技术员以及需要岩性因子参与权重计算的环境建模从业者。它不是扫描件不是需要二次矢量化的图片而是一张已经编码好的分类栅格。2. 岩性栅格的数据组织像元值、配色表与坐标对齐2.1 像元值编码与岩性分类体系拿到一份岩性栅格第一件事不是急着出图而是搞清楚每个像元值到底代表什么岩性。这份数据采用整型编码像元值从1开始顺序编号对应一套岩性分类表。常见做法是随数据附一个.clr配色文件或者.csv对照表里面至少包含三列像元值、岩性名称、岩性大类。岩性大类一般归并到沉积岩、岩浆岩、变质岩三大类再往下细分到具体岩性比如花岗岩、砂岩、灰岩、片麻岩等。如果你拿到的只有栅格没有对照表别慌可以用GDAL先把唯一值拉出来再结合区域地质资料反查。下面这段Python代码就是干这个的from osgeo import gdal import numpy as np # 打开岩性栅格 ds gdal.Open(guangdong_lithology_250m.tif) band ds.GetRasterBand(1) arr band.ReadAsArray() # 统计唯一像元值及其频数 unique, counts np.unique(arr, return_countsTrue) for u, c in zip(unique, counts): print(f像元值 {u}: {c} 个像元, 占比 {c/arr.size*100:.2f}%) # 顺便看一下NoData值是多少 nodata band.GetNoDataValue() print(fNoData 值: {nodata})这段代码的逻辑很直接ReadAsArray()把整幅栅格读进内存np.unique同时返回唯一值和出现次数。参数上唯一需要注意的是NoData值如果数据里用0或者-9999表示无数据统计前最好先掩掉否则会把无数据区当成一种岩性。占比那一列能帮你快速判断分类是否合理——比如广东大面积出露花岗岩如果花岗岩占比过低可能编码表对不上。2.2 WGS84坐标下的像元对齐与重采样WGS84地理坐标系意味着栅格的像元大小是度不是米。250米分辨率在赤道附近约等于0.00225度但广东在北纬20°到25°之间经度方向的实际地面距离会随纬度变化。这一点在做多源数据叠加时特别容易翻车你的DEM如果是UTM投影直接叠加会出现像元错位。常见做法是统一到同一投影再做分析。如果坚持用WGS84那就确保所有栅格的GeoTransform一致。用GDAL查看和设置地理变换参数的代码如下from osgeo import gdal ds gdal.Open(guangdong_lithology_250m.tif) gt ds.GetGeoTransform() print(f左上角坐标: ({gt[0]}, {gt[3]})) print(f像元大小: {gt[1]} 度 x {gt[5]} 度) print(f行列数: {ds.RasterXSize} x {ds.RasterYSize}) # 如果需要对其他栅格做对齐可以用gdal.Warp # 假设要把一个DEM重采样到岩性栅格的网格上 gdal.Warp(dem_aligned.tif, dem_utm.tif, dstSRSEPSG:4326, xResgt[1], yResabs(gt[5]), outputBounds(gt[0], gt[3] gt[5]*ds.RasterYSize, gt[0] gt[1]*ds.RasterXSize, gt[3]), resampleAlgnear) # 分类数据必须用最近邻这里的关键参数是resampleAlg。岩性是分类数据重采样只能用near最近邻用双线性或者立方卷积会把类别值插成小数直接毁掉数据。outputBounds四个值分别是左上角x、左下角y、右下角x、右上角y顺序别搞反。我一般会先把岩性栅格的边界算出来再拿这个边界去裁其他数据保证像元一一对应。2.3 属性表的挂接与出图配色栅格本身不带属性表但做专题图的时候需要把像元值和岩性名称对应起来。QGIS里可以直接加载.clr文件或者手动在图层属性里设置分类渲染。如果要在Python里出图用matplotlib配合ListedColormap更可控import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap, BoundaryNorm import numpy as np from osgeo import gdal ds gdal.Open(guangdong_lithology_250m.tif) arr ds.GetRasterBand(1).ReadAsArray() # 假设有8类岩性自定义配色 colors [#a6cee3, #1f78b4, #b2df8a, #33a02c, #fb9a99, #e31a1c, #fdbf6f, #ff7f00] cmap ListedColormap(colors) bounds np.arange(0.5, 9.5, 1) norm BoundaryNorm(bounds, cmap.N) fig, ax plt.subplots(figsize(10, 8)) im ax.imshow(arr, cmapcmap, normnorm) cbar plt.colorbar(im, axax, ticksrange(1, 9)) cbar.set_label(岩性编码) plt.title(广东地表出露岩性分布250m) plt.savefig(lithology_map.png, dpi300, bbox_inchestight)BoundaryNorm的作用是让每个整数编码落在色带正中间避免出现颜色渐变。ticks设成1到8对应你的分类数。出图之后如果发现某类颜色太浅看不清直接改colors列表里的十六进制值就行不用动数据本身。3. 把岩性栅格接进空间分析流程从重分类到因子提取3.1 岩性因子的重分类与赋值做地质灾害易发性评价或者土壤侵蚀模型时岩性很少直接以原始类别参与计算通常要重分类成工程岩组或者硬度等级。比如把花岗岩、闪长岩归为“坚硬岩组”页岩、泥岩归为“软弱岩组”灰岩单独作为“可溶岩组”。这一步用numpy的映射操作最快import numpy as np from osgeo import gdal ds gdal.Open(guangdong_lithology_250m.tif) arr ds.GetRasterBand(1).ReadAsArray() # 假设原始编码1-花岗岩, 2-砂岩, 3-页岩, 4-灰岩, 5-片麻岩 # 重分类1-1(坚硬), 2-2(较硬), 3-3(软弱), 4-4(可溶), 5-1(坚硬) reclass_map {1: 1, 2: 2, 3: 3, 4: 4, 5: 1} reclassed np.zeros_like(arr) for old, new in reclass_map.items(): reclassed[arr old] new # 写出重分类结果 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(lithology_reclass.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(reclassed) out_ds.GetRasterBand(1).SetNoDataValue(0) out_ds Nonereclass_map字典的键是原始编码值是重分类后的编码。用循环逐类赋值虽然不如np.select快但胜在直观类别少的时候完全够用。写出的数据类型选Byte因为重分类后类别数一般不超过10没必要用Int16浪费空间。SetNoDataValue(0)把背景值标成无数据后续做掩膜分析时不会干扰统计。3.2 与DEM叠加计算坡度-岩性组合单元岩性对边坡稳定性的影响往往要通过和坡度组合才能体现。常见做法是把坡度分成若干等级和岩性大类做交叉制表生成组合单元图。下面这段代码演示如何用rasterio读取两个栅格并做交叉import rasterio import numpy as np with rasterio.open(lithology_reclass.tif) as src_lith: lith src_lith.read(1) profile src_lith.profile with rasterio.open(slope_250m.tif) as src_slope: slope src_slope.read(1) # 坡度分级0-5, 5-15, 15-25, 25-35, 35 slope_class np.digitize(slope, bins[5, 15, 25, 35]) # 组合编码岩性*10 坡度等级 combined lith * 10 slope_class combined[lith 0] 0 # 保持无数据 # 写出组合栅格 profile.update(dtyperasterio.int16, nodata0) with rasterio.open(lith_slope_combo.tif, w, **profile) as dst: dst.write(combined.astype(rasterio.int16), 1)np.digitize返回的是每个像元落在哪个区间bins设了四个边界所以坡度等级是0到4。组合编码用岩性*10 坡度等级解码的时候十位数是岩性个位数是坡度等级。这种编码方式在后续做统计时特别方便用整除和取余就能拆开。注意profile是从岩性栅格继承的保证了输出栅格的地理参考和岩性完全一致。3.3 分区统计按行政区汇总岩性面积如果要做区域地质报告经常需要按市或者县统计各类岩性的面积。用rasterstats可以一行代码搞定但前提是矢量边界和栅格的坐标系一致from rasterstats import zonal_stats import geopandas as gpd # 读取行政区矢量 gdf gpd.read_file(guangdong_counties.shp) # 确保坐标系是WGS84 gdf gdf.to_crs(EPSG:4326) # 分区统计每个县各类岩性的像元数 stats zonal_stats(gdf, guangdong_lithology_250m.tif, categoricalTrue, nodata0) # 把结果挂回GeoDataFrame gdf[lith_counts] [s for s in stats] print(gdf[[NAME, lith_counts]].head())categoricalTrue告诉rasterstats这是分类数据返回的是每个类别出现的像元数而不是均值或者总和。nodata0把背景排除掉。拿到像元数之后乘以每个像元的实际面积250米×250米但要注意纬度变化就能换算成平方公里。我一般会把这个结果导出成CSV再用pandas做透视表按岩性大类汇总。4. 避坑与排查岩性栅格处理中的五个血泪教训4.1 像元值被重采样插成小数现象重采样之后原本整型的岩性编码变成了2.5、3.7这种小数分类统计全乱套。原因用了双线性或者立方卷积插值。这两种算法是为连续数据设计的对分类数据是灾难。解决所有涉及岩性栅格的重采样、投影变换、裁剪操作resampleAlg一律设成near。在QGIS里导出栅格时重采样方法选“最近邻”。4.2 WGS84下直接算面积导致偏差现象按像元数乘以250×250算出来的面积和行政区实际面积对不上越往北偏差越大。原因WGS84是地理坐标系像元的经度跨度对应的地面距离随纬度增加而减小。广东南北跨了5个纬度面积误差能到百分之几。解决做面积统计前先投影到等面积投影比如Albers等积投影参数用中央经线114°E、双标准纬线20°N和25°N。投影之后再统计像元数结果才靠谱。4.3 NoData值参与运算污染结果现象做岩性重分类或者叠加分析时背景区域被赋成了某个岩性类别统计多出一大片。原因没有正确设置或识别NoData值。有些数据用0表示无数据有些用-9999还有些压根没设。解决拿到数据先用gdalinfo或者band.GetNoDataValue()查一下。如果没设根据数据分布判断——岩性编码一般从1开始0大概率是背景。处理前统一用掩膜把NoData盖掉。4.4 分类配色表与像元值错位现象出图之后发现颜色和岩性对不上花岗岩显示成了灰岩的颜色。原因配色文件的类别顺序和栅格像元值编码顺序不一致或者中间漏了某个编码。解决出图前先用唯一值统计确认实际有哪些编码再按编码顺序逐一指定颜色。别直接套用网上找的通用地质配色方案编码体系不一样。4.5 大区域处理时内存溢出现象用ReadAsArray()读全省栅格时程序崩溃提示内存不足。原因250米分辨率全省栅格虽然不算特别大但如果同时读多个图层做运算内存占用会成倍增长。解决用分块读取。GDAL的ReadAsArray支持指定窗口每次只读一块处理完写出再读下一块。或者用rasterio的block_windows迭代。实在不行就换64位Python32位版本单进程最多用2GB内存。5. 进阶技巧用岩性栅格做滑坡易发性快速制图5.1 信息量模型的计算流程滑坡易发性评价里信息量模型是入门门槛最低、解释性又不错的方法。核心思路是统计每个因子类别内滑坡点的分布密度和全区平均密度做对比取对数得到信息量值。岩性作为其中一个因子计算步骤如下import numpy as np import pandas as pd from osgeo import gdal # 读取岩性栅格和滑坡点标签栅格 lith_ds gdal.Open(guangdong_lithology_250m.tif) lith lith_ds.GetRasterBand(1).ReadAsArray() label_ds gdal.Open(landslide_label.tif) # 1为滑坡0为背景 label label_ds.GetRasterBand(1).ReadAsArray() # 统计各类岩性的像元总数和滑坡像元数 classes np.unique(lith) classes classes[classes ! 0] # 去掉NoData total_pixels np.sum(lith ! 0) total_landslides np.sum(label 1) results [] for c in classes: mask (lith c) n_pixels np.sum(mask) n_landslides np.sum(label[mask] 1) # 信息量计算 if n_landslides 0 and n_pixels 0: I np.log((n_landslides / total_landslides) / (n_pixels / total_pixels)) else: I 0 results.append({岩性编码: c, 像元数: n_pixels, 滑坡数: n_landslides, 信息量: round(I, 4)}) df pd.DataFrame(results) print(df.sort_values(信息量, ascendingFalse))这段代码的输出是一张表信息量越大说明该岩性单元越容易发生滑坡。np.log里面是滑坡比例和面积比例的比值如果某个岩性里滑坡占比高于面积占比信息量为正。参数上唯一要注意的是滑坡标签栅格必须和岩性栅格严格对齐否则统计出来的数没有意义。5.2 把信息量值回写到栅格拿到各类岩性的信息量之后下一步是生成一张连续的信息量栅格方便和其他因子的信息量图层叠加# 构建映射字典 info_dict dict(zip(df[岩性编码], df[信息量])) # 回写 info_raster np.zeros_like(lith, dtypenp.float32) for c, v in info_dict.items(): info_raster[lith c] v # 写出 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(lith_info.tif, lith_ds.RasterXSize, lith_ds.RasterYSize, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(lith_ds.GetGeoTransform()) out_ds.SetProjection(lith_ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(info_raster) out_ds.GetRasterBand(1).SetNoDataValue(0) out_ds Noneinfo_raster用float32是因为信息量有小数Byte存不下。SetNoDataValue(0)把背景排除后续做叠加时不会把0当成有效信息量。生成的这张栅格可以直接扔进加权叠加模型和其他因子的信息量图层做栅格加法。5.3 一个容易忽略的验证习惯信息量模型跑完之后别急着出图。我一般会做一件事把滑坡点随机分成训练集和验证集用训练集算信息量用验证集画ROC曲线。如果AUC低于0.7说明岩性这个因子单独区分能力有限得考虑和其他因子组合。这个习惯是从一次翻车经历来的——当时直接用全部滑坡点算信息量出图漂亮得很结果审稿人一问验证指标就露馅了。从那以后我每次做易发性制图都强制走一遍训练/验证拆分哪怕只是快速评估。希望这份岩性栅格的处理思路能帮到你少走几个我踩过的坑。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
