企业官网建设流程全解析

简介：这份资源提供广东省全域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 = "+proj=aea +lat_1=20 +lat_2=26 +lat_0=23 +lon_0=113 +datum=WGS84 +units=m" 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( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=dst_crs, resampling=Resampling.nearest) # 分类栅格必须最近邻

这段代码先读取原始WGS84栅格的元信息，确认坐标系、像元尺寸和NoData值，再投影到以广东为中心的Albers等面积投影。参数说明：lat_1、lat_2是双标准纬线，取20和26覆盖广东南北；lon_0=113是中央经线；units=m保证输出单位是米。重采样方法锁定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_counts=True) 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, on="code", how="left") df["area_km2"] = df["pixels"] * 0.25 * 0.25 # 投影后每像元0.0625 km² print(df.sort_values("area_km2", ascending=False).head(10))

逻辑说明：np.unique同时返回唯一值和计数，比循环快得多。NoData必须剔除，否则会混进统计。面积计算用的是投影后250米×250米的假设，如果还在WGS84下，这个面积只是近似值。关联属性表用merge，how="left"保证即使有码值在属性表里缺失也能发现——缺失的会显示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", categorical=True, # 分类栅格必须开这个 nodata=-9999, geojson_out=False) # stats 是列表，每个元素是 {岩性码: 像元数} for i, s in enumerate(stats): total = sum(s.values()) top3 = sorted(s.items(), key=lambda x: -x[1])[:3] print(f"流域{i}: 主导岩性 {top3[0][0]}, 占比 {top3[0][1]/total:.1%}")

参数说明：categorical=True是关键，它让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, size=3, mode="nearest")

参数说明：generic_filter的size=3对应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值设置正确，categorical=True时NoData不会出现在结果字典里，但分母要自己用有效像元数算。边界歧义可以通过all_touched=False（默认）只算中心点落在多边形内的像元来规避，代价是面积略有偏差。如果精度要求高，先把多边形栅格化再统计。

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 COMPRESS=LZW gd_litho_wgs84.tif gd_litho_byte.tif gdaladdo -r nearest gd_litho_byte.tif 2 4 8 16

-ot Byte把数据类型压到8位无符号整数，前提是岩性码不超过255。-co COMPRESS=LZW做无损压缩。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(dtype="int32", 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个点用在线地质资料核对。这三步花不了半小时，但能避免后面几天的分析建立在错误数据上。岩性栅格不是“下载了就能用”的数据，它的价值取决于你对编码、坐标系和出露含义的理解程度。希望帮到你。

本文还有配套的精品资源，点击获取