广东250米岩性栅格数据：GIS分析、分区统计与机器学习特征工程实战

简介这份资源提供广东省全域250米分辨率的地表出露岩性分布栅格数据采用WGS84地理坐标系共划分14类岩性类型涵盖中性深成岩、中性火山岩、冰川沉积物、各类变质岩、基性深成岩及火山碎屑岩等并按岩石成因归入火成岩、沉积岩、变质岩三大类其中“土”类指未固结的第四纪松散堆积物。分类依据为地表可见岩性不涉及地下深部岩层推断适合从事地质环境评价、水文模拟、土壤侵蚀分析及CO2化学风化消耗量估算等方向的研究人员与高年级学生使用。资源包共13个文件约831KB以tif主栅格为核心辅以dbf属性表、tfw坐标信息、cpg编码文件、aux.xml元数据并附Python脚本、说明文档与预览图便于直接加载与二次处理。目前已有17人学习下载可帮助读者快速获取广东岩性空间分布底图支撑区域地质与地表过程建模。1. 广东全省250米分辨率地表出露岩性分布栅格数据一份能直接进GIS的底图做区域地质背景分析的人大概都经历过这样的场景手头有一份地球化学采样点数据想叠加岩性做分区统计翻遍硬盘只找到几张扫描的纸质地质图配准误差几十米矢量化到一半就放弃了。广东全省250米分辨率地表出露岩性分布栅格数据解决的正是这个卡点——它把全省地表出露的岩性类型按250米×250米的格网编码成一张栅格图坐标系是WGS84可以直接拖进QGIS、ArcGIS或者用Python的rasterio读取和DEM、降水、NDVI、采样点做逐像元叠加。这份数据的价值不在“好看”而在“能算”。矢量地质图做区域统计时边界处的面积权重很难处理一个多边形跨了几十个像元按属性连接会把整个多边形算给一个类别栅格化之后每个像元只有一个岩性码分区统计、缓冲区分析、机器学习特征工程都能直接跑。250米这个尺度对省级尺度的成土母质分析、地下水富水性评价、生态地质背景划分来说精度和体量平衡得比较合适——全省约18万平方公里栅格行列数在可接受范围内普通笔记本也能处理。适合谁用做土壤、水文、生态、地质灾害背景分析的从业者需要给模型加地质先验特征的算法工程师以及带学生做区域实习、不想在数据配准上耗掉半学期的导师。下面从数据本身怎么理解、怎么落地、怎么避坑一路讲到进阶用法。2. 岩性栅格数据的编码逻辑与坐标系统先搞懂像元值代表什么2.1 岩性码不是连续值别拿它做插值拿到栅格第一件事是看像元值分布。岩性数据是分类栅格像元值通常是整数编码比如1代表花岗岩、2代表砂岩、3代表碳酸盐岩具体编码表随数据附带的属性表或说明文件。很多人习惯性打开栅格就做重采样、做插值这是典型的翻车起点——分类栅格用双线性插值重采样会在两类岩性交界处生成不存在的“混合岩性码”比如花岗岩码10和砂岩码20之间插出15而15可能恰好对应另一个真实岩性结果就是凭空多出一片错误分区。正确做法是重采样一律用最近邻nearest neighbor。在QGIS里导出不同分辨率时重采样方法选“最近邻”在Python里用rasterio的Resampling.nearest。如果确实需要平滑也应该先转成各岩性的独热one-hot栅格再做滤波而不是直接对类别码做算术。另一个常见误解是拿岩性码当数值大小比较。码值大小只代表类别编号不代表硬度、年龄或任何物理量的高低。做统计分析时岩性应该作为分类变量处理用众数、类别频率而不是均值。2.2 WGS84坐标系下的像元对齐问题WGS84是一个地理坐标系单位是度不是米。250米分辨率在地理坐标系下并不是固定的度数——赤道处1度约111公里北纬23度附近广东大致纬度范围1度经度约102公里1度纬度约111公里。所以“250米分辨率”在WGS84下经度和纬度方向的像元跨度是不一样的大约是0.00245度经度、0.00225度纬度具体以数据实际像元尺寸为准。这带来两个实操问题。第一和投影坐标系的数据叠加时比如你有一份UTM投影的DEM直接叠加会有轻微错位因为WGS84的像元是“度”的等间隔不是“米”的等间隔。第二做面积统计时不能简单用像元数乘以250×250因为每个像元实际面积随纬度变化。省级尺度上这个误差在百分之几量级如果只是做相对分区还凑合要做绝对面积统计建议先投影到等面积投影如Albers再统计。import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling # 打开WGS84岩性栅格查看基本信息 with rasterio.open(gd_litho_wgs84.tif) as src: print(CRS:, src.crs) # 应为 EPSG:4326 print(尺寸:, src.width, x, src.height) print(像元尺寸:, src.res) # 地理坐标系下单位为度 print(数据类型:, src.dtypes[0]) # 通常为 uint8 或 int16 print(NoData:, src.nodata) # 投影到Albers等面积投影便于面积统计 dst_crs projaea lat_120 lat_226 lat_023 lon_0113 datumWGS84 unitsm with rasterio.open(gd_litho_wgs84.tif) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs src.meta.copy() kwargs.update({crs: dst_crs, transform: transform, width: width, height: height}) with rasterio.open(gd_litho_albers.tif, w, **kwargs) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest) # 分类栅格必须最近邻这段代码先读取原始WGS84栅格的元信息确认坐标系、像元尺寸和NoData值再投影到以广东为中心的Albers等面积投影。参数说明lat_1、lat_2是双标准纬线取20和26覆盖广东南北lon_0113是中央经线unitsm保证输出单位是米。重采样方法锁定Resampling.nearest这是分类栅格投影变换的铁律。投影后像元尺寸会变成接近250米的整数面积统计就准了。2.3 属性表怎么读、怎么关联岩性栅格的像元值本身只是编号真正的岩性名称、大类、成因类型通常在配套的属性表里可能是CSV、DBF或者栅格的颜色表。常见做法是先读栅格的唯一值列表再和属性表做键值关联生成一个“码值—岩性名—大类”的查找字典。import numpy as np import rasterio import pandas as pd with rasterio.open(gd_litho_wgs84.tif) as src: arr src.read(1) nodata src.nodata # 统计各岩性码的像元数 unique, counts np.unique(arr, return_countsTrue) df pd.DataFrame({code: unique, pixels: counts}) if nodata is not None: df df[df[code] ! nodata] # 关联属性表假设有 litho_lookup.csv含 code 和 litho_name 两列 lookup pd.read_csv(litho_lookup.csv) df df.merge(lookup, oncode, howleft) df[area_km2] df[pixels] * 0.25 * 0.25 # 投影后每像元0.0625 km² print(df.sort_values(area_km2, ascendingFalse).head(10))逻辑说明np.unique同时返回唯一值和计数比循环快得多。NoData必须剔除否则会混进统计。面积计算用的是投影后250米×250米的假设如果还在WGS84下这个面积只是近似值。关联属性表用mergehowleft保证即使有码值在属性表里缺失也能发现——缺失的会显示NaN这时候要回去查属性表是不是不全。3. 把岩性栅格接进分析流程分区统计、采样与特征工程3.1 按流域或行政区做岩性分区统计最常见的需求是给定一批多边形流域、行政区、采样缓冲区统计每个多边形内各岩性的面积占比。用rasterstats的zonal_stats可以一行搞定但分类栅格要指定统计函数为majority众数或自定义的频率统计。from rasterstats import zonal_stats import geopandas as gpd zones gpd.read_file(basins.shp) # 流域多边形 # 众数岩性 各类别像元计数 stats zonal_stats( zones, gd_litho_albers.tif, categoricalTrue, # 分类栅格必须开这个 nodata-9999, geojson_outFalse) # stats 是列表每个元素是 {岩性码: 像元数} for i, s in enumerate(stats): total sum(s.values()) top3 sorted(s.items(), keylambda x: -x[1])[:3] print(f流域{i}: 主导岩性 {top3[0][0]}, 占比 {top3[0][1]/total:.1%})参数说明categoricalTrue是关键它让zonal_stats按类别计数而不是求均值。nodata要和栅格一致否则NoData会被当成一个类别。返回的字典键是岩性码需要再用前面的查找表翻译成名称。如果多边形很小、内部像元很少统计结果波动会很大这时候要么合并小多边形要么在报告里标注样本量。3.2 采样点提取岩性属性做地球化学或生态调查时经常需要给每个采样点提取所在位置的岩性。点提取比分区统计简单但要注意坐标系统一和点位落在NoData上的处理。import geopandas as gpd import rasterio points gpd.read_file(samples.shp).to_crs(EPSG:4326) with rasterio.open(gd_litho_wgs84.tif) as src: coords [(p.x, p.y) for p in points.geometry] sampled [x[0] for x in src.sample(coords)] points[litho_code] sampled # 标记落在NoData或边界外的点 points[valid] points[litho_code] ! src.nodata print(points[valid].value_counts())逻辑说明src.sample接受坐标列表返回每个点的像元值。注意采样点必须先转到和栅格相同的坐标系这里栅格是WGS84所以点也转EPSG:4326。如果点落在栅格边界外sample会返回0或NoData所以要显式检查。实际项目中采样点坐标常有偏移如果发现大量点落在水体或城区这些在岩性图上可能是NoData或“人工填土”类要回头核对点位精度。3.3 给机器学习模型加岩性先验特征把岩性作为分类特征喂给模型时直接输入码值会让模型误以为码值有大小关系。正确做法是独热编码或者按岩性大类做目标编码。250米栅格做特征时还要考虑邻域信息——一个点的岩性往往不如“周围1公里内岩性组合”有信息量。import numpy as np import rasterio from scipy.ndimage import generic_filter with rasterio.open(gd_litho_albers.tif) as src: arr src.read(1).astype(np.float32) nodata src.nodata arr[arr nodata] np.nan # 用众数滤波提取3x3邻域主导岩性 def mode_filter(values): vals values[~np.isnan(values)] if len(vals) 0: return np.nan return np.bincount(vals.astype(int)).argmax() dominant generic_filter(arr, mode_filter, size3, modenearest)参数说明generic_filter的size3对应3×3像元在250米分辨率下约覆盖750米×750米。mode_filter里先剔除NaN再取众数避免NoData污染。如果要做多尺度特征可以分别算3×3、5×5、7×7的众数作为不同尺度的岩性背景特征。注意generic_filter比较慢全省范围建议分块处理或者用numba加速。4. 避坑与排查岩性栅格处理中最容易翻车的五件事4.1 现象重采样后岩性类别数变多出现不存在的码值原因用了双线性或三次卷积插值。分类栅格的码值是离散标签插值会在类别之间生成中间值这些中间值可能恰好落在另一个合法码值上。解决所有重采样、投影变换、缩放操作一律用最近邻。在QGIS导出时检查“重采样方法”下拉框在GDAL命令里加-r near在Python里显式指定Resampling.nearest。如果已经生成了错误栅格只能从原始数据重做没有后悔药。4.2 现象分区统计结果里某个流域的岩性占比加起来不是100%原因NoData像元被计入分母或者多边形边界处的像元归属有歧义。zonal_stats默认只统计多边形内部像元但边界像元可能被算给相邻多边形。解决统计前先确认NoData值设置正确categoricalTrue时NoData不会出现在结果字典里但分母要自己用有效像元数算。边界歧义可以通过all_touchedFalse默认只算中心点落在多边形内的像元来规避代价是面积略有偏差。如果精度要求高先把多边形栅格化再统计。4.3 现象采样点提取的岩性和实地不符大量点落在“水体”或“第四系”原因地表出露岩性图反映的是基岩出露情况水体、城区、厚层第四系覆盖区往往被归为单独类别或NoData。采样点如果落在这些区域提取到的不是基岩岩性。解决这是数据本身的属性不是bug。分析时要么剔除这些点要么把“第四系覆盖”作为一个有意义的类别保留——它本身就说明该点没有基岩出露。如果研究目标是基岩应该在采样设计阶段就避开覆盖区或者用地质图交叉验证。4.4 现象WGS84栅格和投影数据叠加时边界对不上偏移几百米原因WGS84是地理坐标系ArcGIS和QGIS在显示时会做动态投影但不同软件的基准面转换参数可能不同导致视觉偏移。另外如果栅格本身没有正确的坐标定义比如只是写了4326但实际是别的基准偏移会更大。解决先确认栅格的CRS定义是否准确用gdalinfo查看。如果要做精确叠加统一投影到同一个投影坐标系再操作不要依赖软件的动态投影。基准面转换参数在省级尺度上影响通常在米级如果偏移达到几百米多半是CRS定义错了。4.5 现象全省栅格太大打开慢、处理卡原因250米全省栅格虽然不算巨大但如果用float32存储、没有金字塔、还带着全范围NoData读写效率会很低。解决把数据类型压到uint8或int16岩性码通常不超过255个类别用gdaladdo建金字塔裁剪到研究区范围再处理用窗口读取rasterio.windows代替全图读取。如果只是做点采样根本不需要读全图src.sample本身就是按需读取。# 压缩数据类型并建金字塔 gdal_translate -ot Byte -co COMPRESSLZW gd_litho_wgs84.tif gd_litho_byte.tif gdaladdo -r nearest gd_litho_byte.tif 2 4 8 16-ot Byte把数据类型压到8位无符号整数前提是岩性码不超过255。-co COMPRESSLZW做无损压缩。gdaladdo建金字塔-r nearest保证金字塔也是最近邻不会在缩放显示时产生假类别。5. 进阶用法岩性栅格与多源数据的耦合分析5.1 岩性—地形—气候的联合分区单一岩性图的信息量有限真正有解释力的是岩性和地形、气候的组合。一个实用技巧是把岩性码、高程带、坡度带做交叉制表生成“岩性—地形组合单元”再统计每个单元的地球化学背景值。这种做法在成土母质研究和地下水富水性评价里很常见。具体操作先把DEM重采样到250米并分带比如200米、200-500米、500米再把坡度分带5°、5-15°、15°然后用栅格计算器把三个分类栅格组合成一个整数码combo litho * 100 elev_band * 10 slope_band。这样每个组合单元有唯一码值后续分区统计就按这个combo做。import rasterio import numpy as np with rasterio.open(gd_litho_albers.tif) as src: litho src.read(1) profile src.profile with rasterio.open(dem_250m.tif) as src: dem src.read(1) with rasterio.open(slope_250m.tif) as src: slope src.read(1) # 分带 elev_band np.digitize(dem, [200, 500]) # 0,1,2 slope_band np.digitize(slope, [5, 15]) # 0,1,2 combo litho.astype(np.int32) * 100 elev_band * 10 slope_band profile.update(dtypeint32, nodata-9999) with rasterio.open(litho_terrain_combo.tif, w, **profile) as dst: dst.write(combo, 1)参数说明np.digitize把连续值分成区间[200, 500]产生0、1、2三个带。组合码用litho*100 elev*10 slope前提是岩性码小于100、分带数小于10否则要调整位权。生成的combo栅格可以像普通分类栅格一样做分区统计但解读时要记得码值的构成规则。5.2 用岩性栅格做采样设计的分层依据如果下一步要布设野外采样点岩性图可以直接做分层抽样。思路是按岩性类别分层每层按面积比例分配采样数再在层内随机或按可达性选点。这样比均匀网格采样更能覆盖稀有岩性。一个容易忽略的点是空间自相关——同一岩性大片连续分布时随机撒点会浪费在相似位置上。改进做法是在每个岩性层内先做空间聚类每个聚类至少采一个点。用sklearn的KMeans对像元坐标聚类或者简单用规则网格在每个岩性层内抽点。5.3 精度验证怎么知道这份岩性图靠不靠谱任何栅格数据都有误差岩性图也不例外。验证方法分两步一是和独立来源比对比如用更高精度的局部地质图做交叉验证看类别一致性二是实地抽查选若干点核对出露岩性。如果拿不到独立验证数据至少做内部一致性检查看岩性边界和地形边界是否合理——比如碳酸盐岩区往往对应喀斯特地貌花岗岩区往往对应圆缓山丘。如果岩性边界和地形完全无关可能是配准出了问题。另一个检查是看类别面积分布是否符合区域地质常识比如广东花岗岩大面积出露如果统计结果里花岗岩占比极低就要怀疑编码或分类有问题。我自己做这类分析时习惯在拿到数据的第一天就做三件事画一张全省岩性图缩略图看空间格局是否合理、统计各类别面积占比是否和已知区域地质背景吻合、随机抽20个点用在线地质资料核对。这三步花不了半小时但能避免后面几天的分析建立在错误数据上。岩性栅格不是“下载了就能用”的数据它的价值取决于你对编码、坐标系和出露含义的理解程度。希望帮到你。本文还有配套的精品资源点击获取