企业官网建设流程全解析

简介：这份资源提供广东省全域250米分辨率的地表出露岩性分布栅格数据，采用WGS84地理坐标系，共划分14类岩性类型，按成因归为火成岩、沉积岩、变质岩三大类，其中“土”类指未固结的第四纪松散堆积物，分类依据为地表可见岩性，不涉及地下深部推断。数据可服务于地质环境评价、水文模拟、土壤侵蚀分析及CO2化学风化消耗量估算等研究场景，适合地质、地理、环境与水文方向的研究人员及学生使用。资源包共13个文件，约831KB，以tif主栅格为核心，辅以dbf属性表、tfw坐标信息、cpg编码文件与aux.xml元数据，另含png预览图、xls说明表、py脚本及使用说明压缩包，便于快速读取与二次处理。目前已有17人学习下载。借助该数据，读者可直接获得覆盖全省的岩性空间底图，结合属性表与说明文档完成分类统计、制图与建模，为区域地质与地表过程研究提供基础支撑。

1. 广东全省250米分辨率地表出露岩性分布栅格数据：从拿到文件到真正用起来

做区域地质、水文模拟或者工程选址的同行，大概率都遇到过同一个尴尬：手头有行政区划、有DEM、有降雨栅格，唯独缺一张能直接叠加分析的岩性底图。纸质地质图翻起来费劲，矢量化又是个无底洞。广东全省250米分辨率地表出露岩性分布栅格数据（WGS84坐标系）解决的正是这个痛点——它把全省地表出露的岩性类型，按250米网格重采样成一张带属性编码的栅格，坐标系是WGS84，能直接和大多数公开地理数据对齐。适合谁用？做滑坡易发性评价的、做地下水补给分区的、做土壤侵蚀建模的，以及需要快速出区域地质背景图的规划人员。它不替代1:5万精度的野外填图，但作为区域尺度的背景底图，省下的矢量化时间够你多跑两轮模型。

2. 岩性栅格到底存了什么：编码体系与250米网格的取舍

2.1 岩性不是连续变量，栅格值本质是分类码

很多人第一次打开这类数据会愣一下：像素值不是0到1的反射率，也不是高程，而是一串整数。这些整数是岩性类型的编码，比如某个值代表花岗岩，另一个值代表砂岩。关键在于，编码本身没有大小意义，21不代表比15“更硬”或“更老”。如果你直接把它当连续变量丢进回归模型，结果一定是错的。常见做法是把它当分类栅格处理，需要参与数值计算时，先转成独热编码或者按工程需求做硬度、渗透性分级。

我一般会先做一件事：拿到数据后立刻查属性表。没有属性表的岩性栅格就是黑匣子。属性表里通常有编码、岩性名称、可能还有年代地层。如果属性表缺失，就得靠区域地质资料反推，这一步偷懒后面全是坑。

2.2 为什么是250米，而不是30米或1公里

250米这个分辨率不是随便定的。广东全省面积约18万平方公里，30米分辨率下全省栅格约2亿个像元，普通工作站打开都费劲，做全省迭代模拟基本不现实。1公里又太粗，珠江三角洲冲积平原和周边丘陵的岩性边界会被抹掉，影响水文响应单元的划分。250米大约对应1:25万地质图的成图精度，既能保留主要岩性单元的边界形态，又把数据量压到单机可处理的范围内。

从数据量估算：广东东西跨度约800公里，南北约600公里，250米网格下大约是3200列乘2400行，不到800万个像元。用16位整型存储，文件大小约15MB左右，压缩后更小。这个体量，笔记本16GB内存就能流畅做全省分析。

2.3 WGS84坐标系意味着什么

WGS84是地理坐标系，单位是度，不是米。这一点直接决定了你不能拿它直接算面积、算距离。比如你想统计某类岩性的分布面积，直接按像元数乘250米乘250米会出错，因为经纬度网格在高纬度地区实际面积会收缩。广东位于北纬20度到25度之间，这个误差虽然不像高纬度那么夸张，但在全省尺度上仍不可忽略。

正确做法是先投影到等面积投影或等距投影，再算面积。常用的是Albers等面积投影，中央经线取东经113度左右，双标准纬线取20度和25度。投影之后再重采样到250米，或者直接在投影坐标系下做面积统计。如果你只是做叠加显示，WGS84可以直接用；一旦涉及定量计算，先投影。

3. 把岩性栅格接进工作流：从裁剪到重分类的完整操作

3.1 用Python读取并检查数据的基本信息

拿到文件后，第一步不是急着分析，而是确认数据本身没问题。下面这段代码用rasterio读取栅格，打印基本信息，并检查有没有异常值。

import rasterio import numpy as np # 打开岩性栅格文件，假设文件名为gd_litho_250m.tif with rasterio.open('gd_litho_250m.tif') as src: # 打印元数据：坐标系、尺寸、波段数、数据类型 print('CRS:', src.crs) print('尺寸:', src.width, 'x', src.height) print('波段数:', src.count) print('数据类型:', src.dtypes[0]) print('无效值:', src.nodata) # 读取第一波段数据 litho = src.read(1) # 统计唯一值和出现频次 unique, counts = np.unique(litho, return_counts=True) for u, c in zip(unique, counts): print(f'编码 {u}: {c} 个像元')

逻辑说明：src.crs确认是不是WGS84；src.nodata告诉你无效值是多少，常见是-9999或0；np.unique统计每个岩性编码有多少像元，如果某个编码像元数极少，可能是边界噪声或分类错误。参数方面，如果nodata是None，说明文件没定义无效值，后续计算要手动屏蔽0或负值。

3.2 按行政区或流域裁剪：mask与clip的区别

全省数据直接跑模型太慢，通常先裁剪到研究区。这里有两个概念容易混：mask是按掩膜提取，保留掩膜范围内的像元，范围外变无效值；clip是按矩形范围裁剪，输出尺寸变小。如果你有研究区的矢量边界，用mask更合适。

import geopandas as gpd from rasterio.mask import mask # 读取研究区矢量边界，假设是shapefile gdf = gpd.read_file('study_area.shp') # 确保矢量坐标系和栅格一致，都是WGS84 if gdf.crs != 'EPSG:4326': gdf = gdf.to_crs('EPSG:4326') # 提取几何对象 geoms = gdf.geometry.values with rasterio.open('gd_litho_250m.tif') as src: # mask函数返回裁剪后的数据和新的transform out_image, out_transform = mask(src, geoms, crop=True) out_meta = src.meta.copy() # 更新元数据 out_meta.update({ 'height': out_image.shape[1], 'width': out_image.shape[2], 'transform': out_transform }) # 写出裁剪结果 with rasterio.open('litho_clipped.tif', 'w', **out_meta) as dest: dest.write(out_image)

逻辑说明：mask的crop=True表示裁剪到掩膜的最小外接矩形，减少输出尺寸；out_meta继承原数据元数据，只改尺寸和transform。参数上，如果矢量边界很大，crop=False会保持原图尺寸，只把范围外设为nodata，适合后续还要按原图对齐的场景。

3.3 岩性重分类：把几十类合并成工程需要的几类

原始岩性编码可能有三四十类，做水文模型时通常只关心渗透性等级。这时候需要重分类。下面用numpy做映射，把编码映射到1到5的渗透性等级。

# 假设原始编码：1-10为松散沉积物，11-20为砂岩页岩，21-30为碳酸盐岩，31-40为火成岩，41以上为变质岩 # 定义映射规则：高渗透、中渗透、低渗透、极低渗透、隔水 reclass_map = {} for i in range(1, 11): reclass_map[i] = 1 # 松散沉积物，高渗透 for i in range(11, 21): reclass_map[i] = 2 # 砂岩页岩，中渗透 for i in range(21, 31): reclass_map[i] = 3 # 碳酸盐岩，中低渗透 for i in range(31, 41): reclass_map[i] = 4 # 火成岩，低渗透 for i in range(41, 100): reclass_map[i] = 5 # 变质岩，极低渗透 # 应用映射 reclassed = np.zeros_like(litho, dtype=np.uint8) for old, new in reclass_map.items(): reclassed[litho == old] = new # 保持无效值 reclassed[litho == src.nodata] = 0

逻辑说明：映射规则必须基于区域地质认识，不能随便定。比如广东的碳酸盐岩主要集中在粤北，岩溶发育区渗透性可能比砂岩还高，这时候分类就要调整。参数上，dtype选uint8足够，因为类别数很少。重分类后建议用颜色表可视化检查一遍，确认边界合理。

4. 避坑指南：岩性栅格处理中最容易翻车的五个地方

4.1 现象：面积统计结果偏大或偏小；原因：直接用经纬度算面积；解决：先投影再统计

这是最常见的翻车点。WGS84下每个像元代表的实际面积随纬度变化，广东南北跨5个纬度，误差能到百分之几。正确做法是用pyproj或gdal的warp先转到Albers等面积投影，再统计像元数乘像元面积。如果只是做相对比较，不投影也能凑合，但一旦要出面积数字，必须投影。

4.2 现象：重采样后岩性边界出现大量混合像元；原因：用双线性或三次卷积插值；解决：分类栅格必须用最近邻

岩性是分类变量，重采样时只能用最近邻法。双线性插值会算出2.5这种不存在的编码，后续映射全乱。在GDAL里用-r near，在rasterio里用Resampling.nearest。这个坑我见过不止一次，有人用默认参数重采样，结果属性表里多出一堆小数编码，排查半天。

4.3 现象：裁剪后像元值全变成0；原因：矢量边界和栅格坐标系不一致；解决：先统一CRS再裁剪

如果矢量是投影坐标系而栅格是WGS84，mask函数不会报错，但提取结果全是nodata。排查方法是打印两者的crs对比。统一到WGS84再操作，或者把栅格转到矢量坐标系。注意转换会引入重采样，分类栅格同样用最近邻。

4.4 现象：属性表丢失，编码无法解读；原因：格式转换时未保留RAT；解决：转格式前导出属性表为CSV备份

GeoTIFF可以内嵌栅格属性表，但转成某些格式或经过某些软件处理后，属性表会丢。我的习惯是拿到数据先导出属性表为CSV，单独保存。这样即使栅格属性丢了，还能靠编码对照表恢复。如果原始数据就没有属性表，只能查区域地质志或1:25万地质图图例。

4.5 现象：全省分析内存溢出；原因：一次性读入全分辨率数据做浮点运算；解决：分块读取或先降采样

800万像元不算大，但如果做多图层叠加、迭代运算，内存会迅速膨胀。rasterio支持窗口读取，可以分块处理。或者先用out_shape参数降采样读取，做快速预览，确认流程没问题再全分辨率跑。参数上，out_shape设为原尺寸的一半，内存占用降到四分之一。

5. 进阶用法：把岩性栅格变成模型能吃的因子

5.1 计算岩性多样性指数：一个被低估的区域地质指标

岩性多样性指数反映一定窗口内岩性类型的丰富程度，对滑坡易发性、地下水富集性评价有指示意义。用移动窗口统计唯一值数量即可。下面用scipy的generic_filter实现。

from scipy.ndimage import generic_filter def litho_diversity(window): # 排除无效值0 valid = window[window > 0] if len(valid) == 0: return 0 return len(np.unique(valid)) # 用5x5窗口计算多样性 diversity = generic_filter(reclassed, litho_diversity, size=5)

逻辑说明：窗口大小决定多样性指数的空间尺度，5x5在250米分辨率下对应1.25公里见方，适合区域评价。参数上，窗口越大指数越平滑，但会掩盖小范围岩性变化。建议做多个尺度对比，选与研究目标最相关的。

5.2 岩性接触带提取：两种岩性的边界往往是工程关注点

岩性接触带是风化、蚀变、地下水富集的高发区。提取方法是先做边缘检测，再筛选出两侧岩性不同的边界。用numpy的梯度函数可以快速实现。

# 计算水平和垂直梯度 grad_y, grad_x = np.gradient(reclassed.astype(float)) # 梯度非零处即为边界 edge = np.sqrt(grad_x**2 + grad_y**2) > 0 # 进一步筛选：只保留两侧岩性不同的边界 # 用roll做邻域比较 contact = np.zeros_like(edge, dtype=bool) for dy in [-1, 0, 1]: for dx in [-1, 0, 1]: if dy == 0 and dx == 0: continue neighbor = np.roll(np.roll(reclassed, dy, axis=0), dx, axis=1) contact |= edge & (neighbor != reclassed) & (neighbor > 0) & (reclassed > 0)

逻辑说明：np.gradient在分类栅格上会产生阶梯状梯度，所以用阈值判断边界。np.roll做邻域比较时要注意边界环绕问题，实际应用中可以先padding再计算。参数上，接触带宽度可以通过膨胀操作控制，一般膨胀1到2个像元即可。

5.3 验证方法：用已知地质点做混淆矩阵

岩性栅格准不准，不能只看图。如果有野外地质点或钻孔数据，可以做混淆矩阵。把点位坐标提取对应栅格编码，和实际岩性对比，算总体精度和Kappa系数。如果精度低于70%，要么是重采样引入误差，要么是原始数据分类体系和你用的不一致。这时候需要回到属性表核对，而不是硬着头皮往下做。

我自己的习惯是：拿到任何岩性栅格，先裁一小块做可视化，叠加已知地质点看吻合度。吻合度不行，后面所有分析都是沙上建塔。这个250米数据在平原区精度尚可，山区岩性突变处会有混合像元，用的时候心里要有数。希望帮到你。

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