简介这份资源为2020年广东省10米精度土地覆盖土地利用数据面向地理信息、遥感、城市规划及生态环境研究等方向的学习者与从业者。数据基于哨兵影像与深度学习方法制作经墨卡托投影转WGS84地理坐标系并按最新省市级行政边界裁剪覆盖广东省各地级市共分耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪冰等十类。压缩包共147个文件约75.8MB包含21个tif栅格主数据、21个xlsx属性表、21个dbf与cpg等配套文件、21个tfw坐标信息、21个xml元数据及21个png预览图便于直接加载与属性查询。已有329人学习下载。读者可获得统一坐标系、边界规范的广东省市级土地利用数据用于空间分析、制图与变化监测省去投影转换与裁剪步骤。1. 2020年10m精度广东省土地覆盖土地利用数据这份栅格包到底能干什么如果你手头正好缺一份能直接扔进 GIS 里做分析的广东省地表覆盖底图又不想花几天时间去拼 Landsat 或者等哨兵影像慢慢下载那这份「2020年10m精度广东省土地覆盖土地利用.rar」大概率就是你要找的东西。它本质上是一套已经分类好的栅格数据产品空间分辨率 10 米覆盖广东全省时间截面锁定在 2020 年。10 米这个精度意味着什么一个像素代表地面上 10 米 × 10 米的方块广州珠江新城一个标准地铁站出入口的面积大概也就几十个像素用来做市域、县域尺度的国土空间分析、生态评估、城市扩张监测完全够用。它适合做自然资源调查的技术人员、做城市研究的规划师、写论文需要土地利用底图的研究生以及做遥感应用开发需要一份现成分类结果做验证的工程师。但先别急着解压这份数据能不能用、怎么用、坑在哪得先把它拆开看明白。2. 10米分辨率土地覆盖栅格从分类体系到广东地类的对应关系2.1 为什么是10米而不是30米或1米土地覆盖产品的分辨率选择从来不是拍脑袋定的。30 米级别的产品比如早期基于 Landsat 的全球地表覆盖在省级尺度上做宏观趋势没问题但一旦落到珠三角这种城市蔓延剧烈、地块破碎度高的区域30 米像素会把一个城中村和旁边的工业园混成一个类别边界模糊得让人抓狂。1 米级别的高分影像当然精细但广东省全域覆盖的高分数据获取成本高、处理量大而且分类精度受阴影和建筑纹理干扰反而可能下降。10 米恰好卡在一个甜点上哨兵二号Sentinel-2的多光谱波段原生就是 10 米、20 米和 60 米其中可见光和近红外是 10 米用这些波段做分类既保留了城市内部地块的区分度又不至于让数据量爆炸。这份 2020 年的产品大概率就是基于哨兵二号时序影像生产的这也是目前省级 10 米土地覆盖产品的主流技术路线。2.2 分类体系怎么读一级类和二级类的映射拿到栅格后第一个要搞清楚的是它的分类编码。国内常见的土地覆盖产品通常参考《土地利用现状分类》国标或者类似 FROM-GLC、CLCD 等产品的分类体系。这份数据大概率包含耕地、林地、草地、灌木地、湿地、水体、不透水面、裸地这几大类。你需要找到随数据附带的说明文件或者属性表确认每个像元值对应什么类别。常见做法是值 1 代表耕地2 代表林地3 代表草地4 代表灌木5 代表湿地6 代表水体7 代表不透水面8 代表裸地。但不同生产单位编码顺序可能不同所以务必以实际说明为准。如果压缩包里没有说明文件那就得靠目视比对打开珠三角区域看哪个值大面积出现在城市核心区那个值大概率就是不透水面。2.3 数据格式与坐标系统确认解压后你大概率会看到一个或多个 .tif 文件可能按地级市分幅也可能是一整幅广东全省的镶嵌影像。坐标系统通常是 WGS84 地理坐标系EPSG:4326或者 CGCS2000 投影坐标系。如果是 4326单位是度做面积统计时需要先投影到合适的投影坐标系比如 CGCS2000 3 度带高斯克吕格投影否则算出来的面积是平方度毫无意义。检查方法很简单在 QGIS 或 ArcGIS 里加载后看图层属性里的 CRS 信息。如果是分幅的还需要先镶嵌成一整幅再裁剪到研究区边界。# 用 gdalinfo 查看栅格基本信息确认坐标系、行列数、像元大小和 NoData 值 gdalinfo gd_landcover_2020.tif # 关键输出解读 # Size is 120000, 85000 - 全省幅面行列数很大 # Pixel Size (0.0000898, -0.0000898) - 约10米分辨率地理坐标系 # Coordinate System is WGS 84 - EPSG:4326 # NoData Value 0 - 背景值统计时要排除上面这段gdalinfo是拿到任何栅格数据后的第一步。重点看四个东西行列数判断数据范围是否符合预期像元大小确认是不是真的 10 米坐标系决定后续要不要投影转换NoData 值决定统计时怎么设掩膜。如果 NoData 是 0 而你的分类编码里 0 恰好代表某种地类那就得小心了需要跟数据说明核对清楚。2.4 用 Python 快速做面积统计确认完基本信息后最常做的操作就是统计各地类面积。下面这段代码用 rasterio 和 numpy 实现逻辑清晰适合直接抄。import rasterio import numpy as np # 打开栅格文件 with rasterio.open(gd_landcover_2020.tif) as src: data src.read(1) # 读第一个波段 transform src.transform nodata src.nodata # 计算单个像元面积平方米 # transform[0] 是 x 方向像元大小度transform[4] 是 y 方向像元大小 # 地理坐标系下需要换算成米粗略按纬度30度附近1度约111km和96km pixel_width_m transform[0] * 111000 * np.cos(np.radians(23.5)) # 广东平均纬度约23.5 pixel_height_m abs(transform[4]) * 111000 pixel_area_m2 pixel_width_m * pixel_height_m # 排除 NoData valid data ! nodata unique, counts np.unique(data[valid], return_countsTrue) # 输出各地类面积平方公里 for cls, cnt in zip(unique, counts): area_km2 cnt * pixel_area_m2 / 1e6 print(f类别 {cls}: {area_km2:.2f} 平方公里)这段代码的核心逻辑是先算出一个像元代表多少平方米再统计每个类别有多少个像元相乘得到总面积。参数方面transform[0]和transform[4]是仿射变换里的像元尺寸np.cos(np.radians(23.5))是因为地理坐标系下经度方向的实际距离随纬度变化广东跨纬度约 20.5 到 25.5 度取中间值 23.5 度做近似。如果你追求更精确的面积建议先用gdalwarp投影到等面积投影再做统计那样每个像元面积恒定不用做余弦校正。3. 把栅格用起来裁剪、重分类与叠加分析的操作链3.1 按行政区裁剪用掩膜提取研究区全省数据直接分析往往没必要你通常只关心某个市或某个流域。裁剪有两种常见做法用矢量边界做掩膜提取或者按坐标范围做窗口读取。前者更精准后者更快。下面用gdalwarp做矢量裁剪。# 用广州市行政边界裁剪土地覆盖栅格 # -cutline 指定矢量边界文件 # -crop_to_cutline 让输出范围严格贴合边界 # -dstnodata 设置输出背景值 gdalwarp -cutline guangzhou_boundary.shp \ -crop_to_cutline \ -dstnodata 0 \ -tr 0.0000898 0.0000898 \ gd_landcover_2020.tif \ gz_landcover_2020.tif-cutline后面跟矢量文件路径支持 shp、geojson 等格式。-crop_to_cutline是关键参数不加的话输出还是全省范围只是边界外被设为 NoData文件大小没变。-tr指定输出分辨率这里保持和原数据一致避免重采样引入误差。裁剪完之后再用前面的 Python 脚本统计广州各地类面积就方便多了。3.2 重分类把二级类合并成你需要的一级类原始数据的分类可能很细比如把耕地分成水田和旱地把林地分成有林地和灌木林。但你的分析可能只需要「耕地、林地、水体、建设用地」这几大类。这时候就需要重分类。用gdal_calc或者 Python 的 numpy 都能做。import rasterio import numpy as np with rasterio.open(gz_landcover_2020.tif) as src: data src.read(1) profile src.profile # 假设原始编码1水田 2旱地 3有林地 4灌木 5水体 6建设用地 7裸地 # 重分类规则1,2-1(耕地) 3,4-2(林地) 5-3(水体) 6-4(建设用地) 7-5(其他) reclass_map {1:1, 2:1, 3:2, 4:2, 5:3, 6:4, 7:5} reclassified np.zeros_like(data) for old, new in reclass_map.items(): reclassified[data old] new # 保持 NoData reclassified[data 0] 0 profile.update(dtyperasterio.uint8, nodata0) with rasterio.open(gz_landcover_reclass.tif, w, **profile) as dst: dst.write(reclassified, 1)重分类的逻辑就是建立一个旧值到新值的映射字典然后遍历每个旧值把对应像元赋成新值。注意np.zeros_like(data)创建的全零数组如果原始数据里有 0 值代表某种地类这里就会混淆所以务必确认 NoData 和有效类别不重叠。profile.update里把数据类型改成 uint8因为重分类后类别数变少了用不着原来的 uint16 或 float能省不少存储空间。3.3 叠加分析土地覆盖与人口格网的交叉统计土地覆盖数据很少单独用通常要和其他空间数据叠加。比如你想知道广州市每个公里格网里建设用地占比和人口密度的关系就需要把土地覆盖栅格和人口格网做分区统计。常见做法是先用gdalwarp把两套数据重采样到同一分辨率和坐标系然后用 Python 做逐像元计算或者用 QGIS 的 Zonal Statistics 插件。import rasterio import numpy as np import pandas as pd # 读取土地覆盖和人口栅格确保两者已经对齐 with rasterio.open(gz_landcover_reclass.tif) as src: landcover src.read(1) with rasterio.open(gz_population_2020.tif) as src: population src.read(1) # 创建 1km 格网编号假设分辨率约10米100个像元约1km block_size 100 rows, cols landcover.shape results [] for i in range(0, rows, block_size): for j in range(0, cols, block_size): lc_block landcover[i:iblock_size, j:jblock_size] pop_block population[i:iblock_size, j:jblock_size] if lc_block.size 0: continue # 计算建设用地类别4占比 built_ratio np.sum(lc_block 4) / lc_block.size # 计算平均人口 mean_pop np.mean(pop_block) results.append({built_ratio: built_ratio, mean_pop: mean_pop}) df pd.DataFrame(results) print(df.corr()) # 输出相关系数矩阵这段代码把影像切成 100×100 像元的块每块约 1 公里见方然后统计每块里建设用地的比例和平均人口最后算相关系数。block_size可以根据你的研究尺度调整想粗一点就设 200 或 500。注意人口栅格如果是 WorldPop 或 GPW 这类产品单位可能是每像元人数直接取平均没问题如果是密度值那就要乘以面积换算成人数再统计。3.4 变化检测的前置准备和 2010 年数据对齐如果你手头还有 2010 年的同类型数据想做十年变化检测那第一步是确保两期数据坐标系、分辨率、范围完全一致。常见做法是以 2020 年数据为基准把 2010 年数据重采样和裁剪到相同网格。# 将2010年数据对齐到2020年的网格 gdalwarp -t_srs EPSG:4326 \ -tr 0.0000898 0.0000898 \ -te 109.5 20.5 117.5 25.5 \ -r near \ gd_landcover_2010.tif \ gd_landcover_2010_aligned.tif-t_srs指定目标坐标系-tr指定目标分辨率-te指定目标范围xmin ymin xmax ymax-r near表示用最近邻重采样这对分类数据是必须的不能用双线性或立方卷积否则会出现不存在的类别值。对齐之后两期数据逐像元相减非零的地方就是变化区域。4. 避坑与排查10米土地覆盖数据常见的五个翻车点4.1 面积算出来偏大或偏小现象用地理坐标系直接统计面积发现各地类加总跟官方公布的行政区面积对不上偏差能到百分之十几。原因地理坐标系下像元面积随纬度变化高纬度地区像元实际面积比低纬度小用固定值乘会引入系统误差。解决先投影到等面积投影或高斯克吕格投影再做统计。广东常用 CGCS2000 3 度带中央经线根据研究区经度选择比如 114°E 对应 EPSG:4547 左右。投影后像元面积恒定统计结果才可靠。4.2 分类编码对不上说明书现象按说明书里的编码去提取水体结果提取出来的是大片山区。原因不同生产批次或不同来源的数据编码可能调整过说明书没同步更新。解决不要盲信文档先做目视验证。在 QGIS 里加载栅格叠加天地图或影像底图找几个已知地物点比如广州塔附近应该是建设用地珠江应该是水体查看像元值反推编码含义。确认后再批量处理。4.3 分幅数据接边处出现裂缝现象把各地级市分幅数据镶嵌后发现市界处有一条明显的类别突变线。原因分幅生产时各幅独立分类接边处分类结果不一致或者镶嵌时没有做羽化处理。解决如果数据是分幅的优先找全省整幅版本。如果没有镶嵌后用众数滤波做一下后处理或者在接受范围内忽略接边处的少量误差。做变化检测时尤其要注意接边处的假变化会干扰结果。4.4 NoData 值被当成有效类别统计现象统计结果里多出一个面积巨大的未知类别或者面积加总远超全省面积。原因NoData 值通常是 0 或 255没有排除被当成一个地类参与了统计。解决统计前先确认 NoData 值用data ! nodata做掩膜。如果 NoData 是 0 而 0 又恰好是某个地类的编码那就得回去找生产方确认或者用栅格边缘的 0 值区域判断——通常 NoData 会出现在影像边界外。4.5 重采样方法选错导致类别污染现象把 10 米数据重采样到 30 米后出现了一些原始数据里没有的类别值。原因用了双线性或立方卷积等连续型重采样方法把类别值当连续值插值了。解决分类栅格的重采样必须用最近邻法nearest neighbor。gdalwarp里加-r nearPython 里用rasterio的Resampling.nearest。如果确实需要降分辨率也可以先用众数滤波再抽样效果比直接最近邻更好。5. 进阶技巧用土地覆盖数据做生态质量评价的完整链路5.1 从土地覆盖到生态指数土地覆盖数据本身是基础底图真正体现价值的是用它算出衍生指标。一个经典应用是计算区域生态质量指数比如基于土地利用的景观生态风险指数或者生境质量模型。这里给一个简化但可复现的链路用土地覆盖计算香农多样性指数再结合不透水面比例做生态质量分级。import rasterio import numpy as np from scipy.stats import entropy with rasterio.open(gz_landcover_reclass.tif) as src: data src.read(1) profile src.profile # 滑动窗口计算香农多样性指数 window_size 50 # 500米窗口 rows, cols data.shape shannon np.zeros_like(data, dtypenp.float32) for i in range(0, rows, window_size): for j in range(0, cols, window_size): block data[i:iwindow_size, j:jwindow_size] if block.size 0: continue # 统计各类别出现频率 unique, counts np.unique(block, return_countsTrue) # 排除 NoData mask unique ! 0 if np.sum(mask) 0: continue probs counts[mask] / np.sum(counts[mask]) shannon[i:iwindow_size, j:jwindow_size] entropy(probs) profile.update(dtyperasterio.float32, nodata-9999) with rasterio.open(gz_shannon.tif, w, **profile) as dst: dst.write(shannon, 1)这段代码用 50×50 像元的滑动窗口计算香农多样性指数窗口约 500 米见方。entropy函数来自 scipy输入是各类别的概率分布。指数越高说明窗口内土地覆盖类型越丰富生态多样性越好。计算完后可以按自然断点法分成高、中、低几档再和不透水面比例做叠加识别出「高多样性但高不透水面」的冲突区域这些往往是城市扩张前沿值得重点关注。5.2 验证分类精度的一个土办法如果你手头没有独立的验证样本但又想对数据质量有个基本判断可以用一个土办法找几幅同时期的高分影像截图在土地覆盖栅格上随机撒点目视比对。具体操作是在 QGIS 里加载土地覆盖和影像底图创建 100 个随机点逐个查看点所在位置的影像地物和栅格类别是否一致。如果一致率低于 80%说明这份数据在你的研究区可能不太靠谱需要谨慎使用或者考虑自己重新分类。这个方法虽然粗糙但比盲目信任数据强。我一般会在正式分析前花半小时做这个抽查翻车过几次之后就养成了习惯。5.3 数据融合的边界最后说一个容易上头的地方有人会想把这份 10 米土地覆盖和夜间灯光、POI、路网数据全部融合起来做城市边界提取。思路没错但要注意尺度匹配。夜间灯光数据常见的是 500 米分辨率POI 是点数据路网是矢量线直接和 10 米栅格叠加会产生大量空值和尺度不一致问题。常见做法是先把所有数据统一到同一个格网比如 100 米或 250 米再做多要素加权。不要试图在 10 米尺度上融合所有数据那样计算量大不说精度提升也有限。从那以后我每次做多源融合之前都强制先跑一遍尺度一致性检查确认所有数据的空间支撑匹配了再往下走。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/3 18:29:46 · 阅读 1.2k 2020年10米精度广东省土地覆盖栅格数据实战指南
简介这份资源为2020年广东省10米精度土地覆盖土地利用数据面向地理信息、遥感、城市规划及生态环境研究等方向的学习者与从业者。数据基于哨兵影像与深度学习方法制作经墨卡托投影转WGS84地理坐标系并按最新省市级行政边界裁剪覆盖广东省各地级市共分耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪冰等十类。压缩包共147个文件约75.8MB包含21个tif栅格主数据、21个xlsx属性表、21个dbf与cpg等配套文件、21个tfw坐标信息、21个xml元数据及21个png预览图便于直接加载与属性查询。已有329人学习下载。读者可获得统一坐标系、边界规范的广东省市级土地利用数据用于空间分析、制图与变化监测省去投影转换与裁剪步骤。1. 2020年10m精度广东省土地覆盖土地利用数据这份栅格包到底能干什么如果你手头正好缺一份能直接扔进 GIS 里做分析的广东省地表覆盖底图又不想花几天时间去拼 Landsat 或者等哨兵影像慢慢下载那这份「2020年10m精度广东省土地覆盖土地利用.rar」大概率就是你要找的东西。它本质上是一套已经分类好的栅格数据产品空间分辨率 10 米覆盖广东全省时间截面锁定在 2020 年。10 米这个精度意味着什么一个像素代表地面上 10 米 × 10 米的方块广州珠江新城一个标准地铁站出入口的面积大概也就几十个像素用来做市域、县域尺度的国土空间分析、生态评估、城市扩张监测完全够用。它适合做自然资源调查的技术人员、做城市研究的规划师、写论文需要土地利用底图的研究生以及做遥感应用开发需要一份现成分类结果做验证的工程师。但先别急着解压这份数据能不能用、怎么用、坑在哪得先把它拆开看明白。2. 10米分辨率土地覆盖栅格从分类体系到广东地类的对应关系2.1 为什么是10米而不是30米或1米土地覆盖产品的分辨率选择从来不是拍脑袋定的。30 米级别的产品比如早期基于 Landsat 的全球地表覆盖在省级尺度上做宏观趋势没问题但一旦落到珠三角这种城市蔓延剧烈、地块破碎度高的区域30 米像素会把一个城中村和旁边的工业园混成一个类别边界模糊得让人抓狂。1 米级别的高分影像当然精细但广东省全域覆盖的高分数据获取成本高、处理量大而且分类精度受阴影和建筑纹理干扰反而可能下降。10 米恰好卡在一个甜点上哨兵二号Sentinel-2的多光谱波段原生就是 10 米、20 米和 60 米其中可见光和近红外是 10 米用这些波段做分类既保留了城市内部地块的区分度又不至于让数据量爆炸。这份 2020 年的产品大概率就是基于哨兵二号时序影像生产的这也是目前省级 10 米土地覆盖产品的主流技术路线。2.2 分类体系怎么读一级类和二级类的映射拿到栅格后第一个要搞清楚的是它的分类编码。国内常见的土地覆盖产品通常参考《土地利用现状分类》国标或者类似 FROM-GLC、CLCD 等产品的分类体系。这份数据大概率包含耕地、林地、草地、灌木地、湿地、水体、不透水面、裸地这几大类。你需要找到随数据附带的说明文件或者属性表确认每个像元值对应什么类别。常见做法是值 1 代表耕地2 代表林地3 代表草地4 代表灌木5 代表湿地6 代表水体7 代表不透水面8 代表裸地。但不同生产单位编码顺序可能不同所以务必以实际说明为准。如果压缩包里没有说明文件那就得靠目视比对打开珠三角区域看哪个值大面积出现在城市核心区那个值大概率就是不透水面。2.3 数据格式与坐标系统确认解压后你大概率会看到一个或多个 .tif 文件可能按地级市分幅也可能是一整幅广东全省的镶嵌影像。坐标系统通常是 WGS84 地理坐标系EPSG:4326或者 CGCS2000 投影坐标系。如果是 4326单位是度做面积统计时需要先投影到合适的投影坐标系比如 CGCS2000 3 度带高斯克吕格投影否则算出来的面积是平方度毫无意义。检查方法很简单在 QGIS 或 ArcGIS 里加载后看图层属性里的 CRS 信息。如果是分幅的还需要先镶嵌成一整幅再裁剪到研究区边界。# 用 gdalinfo 查看栅格基本信息确认坐标系、行列数、像元大小和 NoData 值 gdalinfo gd_landcover_2020.tif # 关键输出解读 # Size is 120000, 85000 - 全省幅面行列数很大 # Pixel Size (0.0000898, -0.0000898) - 约10米分辨率地理坐标系 # Coordinate System is WGS 84 - EPSG:4326 # NoData Value 0 - 背景值统计时要排除上面这段gdalinfo是拿到任何栅格数据后的第一步。重点看四个东西行列数判断数据范围是否符合预期像元大小确认是不是真的 10 米坐标系决定后续要不要投影转换NoData 值决定统计时怎么设掩膜。如果 NoData 是 0 而你的分类编码里 0 恰好代表某种地类那就得小心了需要跟数据说明核对清楚。2.4 用 Python 快速做面积统计确认完基本信息后最常做的操作就是统计各地类面积。下面这段代码用 rasterio 和 numpy 实现逻辑清晰适合直接抄。import rasterio import numpy as np # 打开栅格文件 with rasterio.open(gd_landcover_2020.tif) as src: data src.read(1) # 读第一个波段 transform src.transform nodata src.nodata # 计算单个像元面积平方米 # transform[0] 是 x 方向像元大小度transform[4] 是 y 方向像元大小 # 地理坐标系下需要换算成米粗略按纬度30度附近1度约111km和96km pixel_width_m transform[0] * 111000 * np.cos(np.radians(23.5)) # 广东平均纬度约23.5 pixel_height_m abs(transform[4]) * 111000 pixel_area_m2 pixel_width_m * pixel_height_m # 排除 NoData valid data ! nodata unique, counts np.unique(data[valid], return_countsTrue) # 输出各地类面积平方公里 for cls, cnt in zip(unique, counts): area_km2 cnt * pixel_area_m2 / 1e6 print(f类别 {cls}: {area_km2:.2f} 平方公里)这段代码的核心逻辑是先算出一个像元代表多少平方米再统计每个类别有多少个像元相乘得到总面积。参数方面transform[0]和transform[4]是仿射变换里的像元尺寸np.cos(np.radians(23.5))是因为地理坐标系下经度方向的实际距离随纬度变化广东跨纬度约 20.5 到 25.5 度取中间值 23.5 度做近似。如果你追求更精确的面积建议先用gdalwarp投影到等面积投影再做统计那样每个像元面积恒定不用做余弦校正。3. 把栅格用起来裁剪、重分类与叠加分析的操作链3.1 按行政区裁剪用掩膜提取研究区全省数据直接分析往往没必要你通常只关心某个市或某个流域。裁剪有两种常见做法用矢量边界做掩膜提取或者按坐标范围做窗口读取。前者更精准后者更快。下面用gdalwarp做矢量裁剪。# 用广州市行政边界裁剪土地覆盖栅格 # -cutline 指定矢量边界文件 # -crop_to_cutline 让输出范围严格贴合边界 # -dstnodata 设置输出背景值 gdalwarp -cutline guangzhou_boundary.shp \ -crop_to_cutline \ -dstnodata 0 \ -tr 0.0000898 0.0000898 \ gd_landcover_2020.tif \ gz_landcover_2020.tif-cutline后面跟矢量文件路径支持 shp、geojson 等格式。-crop_to_cutline是关键参数不加的话输出还是全省范围只是边界外被设为 NoData文件大小没变。-tr指定输出分辨率这里保持和原数据一致避免重采样引入误差。裁剪完之后再用前面的 Python 脚本统计广州各地类面积就方便多了。3.2 重分类把二级类合并成你需要的一级类原始数据的分类可能很细比如把耕地分成水田和旱地把林地分成有林地和灌木林。但你的分析可能只需要「耕地、林地、水体、建设用地」这几大类。这时候就需要重分类。用gdal_calc或者 Python 的 numpy 都能做。import rasterio import numpy as np with rasterio.open(gz_landcover_2020.tif) as src: data src.read(1) profile src.profile # 假设原始编码1水田 2旱地 3有林地 4灌木 5水体 6建设用地 7裸地 # 重分类规则1,2-1(耕地) 3,4-2(林地) 5-3(水体) 6-4(建设用地) 7-5(其他) reclass_map {1:1, 2:1, 3:2, 4:2, 5:3, 6:4, 7:5} reclassified np.zeros_like(data) for old, new in reclass_map.items(): reclassified[data old] new # 保持 NoData reclassified[data 0] 0 profile.update(dtyperasterio.uint8, nodata0) with rasterio.open(gz_landcover_reclass.tif, w, **profile) as dst: dst.write(reclassified, 1)重分类的逻辑就是建立一个旧值到新值的映射字典然后遍历每个旧值把对应像元赋成新值。注意np.zeros_like(data)创建的全零数组如果原始数据里有 0 值代表某种地类这里就会混淆所以务必确认 NoData 和有效类别不重叠。profile.update里把数据类型改成 uint8因为重分类后类别数变少了用不着原来的 uint16 或 float能省不少存储空间。3.3 叠加分析土地覆盖与人口格网的交叉统计土地覆盖数据很少单独用通常要和其他空间数据叠加。比如你想知道广州市每个公里格网里建设用地占比和人口密度的关系就需要把土地覆盖栅格和人口格网做分区统计。常见做法是先用gdalwarp把两套数据重采样到同一分辨率和坐标系然后用 Python 做逐像元计算或者用 QGIS 的 Zonal Statistics 插件。import rasterio import numpy as np import pandas as pd # 读取土地覆盖和人口栅格确保两者已经对齐 with rasterio.open(gz_landcover_reclass.tif) as src: landcover src.read(1) with rasterio.open(gz_population_2020.tif) as src: population src.read(1) # 创建 1km 格网编号假设分辨率约10米100个像元约1km block_size 100 rows, cols landcover.shape results [] for i in range(0, rows, block_size): for j in range(0, cols, block_size): lc_block landcover[i:iblock_size, j:jblock_size] pop_block population[i:iblock_size, j:jblock_size] if lc_block.size 0: continue # 计算建设用地类别4占比 built_ratio np.sum(lc_block 4) / lc_block.size # 计算平均人口 mean_pop np.mean(pop_block) results.append({built_ratio: built_ratio, mean_pop: mean_pop}) df pd.DataFrame(results) print(df.corr()) # 输出相关系数矩阵这段代码把影像切成 100×100 像元的块每块约 1 公里见方然后统计每块里建设用地的比例和平均人口最后算相关系数。block_size可以根据你的研究尺度调整想粗一点就设 200 或 500。注意人口栅格如果是 WorldPop 或 GPW 这类产品单位可能是每像元人数直接取平均没问题如果是密度值那就要乘以面积换算成人数再统计。3.4 变化检测的前置准备和 2010 年数据对齐如果你手头还有 2010 年的同类型数据想做十年变化检测那第一步是确保两期数据坐标系、分辨率、范围完全一致。常见做法是以 2020 年数据为基准把 2010 年数据重采样和裁剪到相同网格。# 将2010年数据对齐到2020年的网格 gdalwarp -t_srs EPSG:4326 \ -tr 0.0000898 0.0000898 \ -te 109.5 20.5 117.5 25.5 \ -r near \ gd_landcover_2010.tif \ gd_landcover_2010_aligned.tif-t_srs指定目标坐标系-tr指定目标分辨率-te指定目标范围xmin ymin xmax ymax-r near表示用最近邻重采样这对分类数据是必须的不能用双线性或立方卷积否则会出现不存在的类别值。对齐之后两期数据逐像元相减非零的地方就是变化区域。4. 避坑与排查10米土地覆盖数据常见的五个翻车点4.1 面积算出来偏大或偏小现象用地理坐标系直接统计面积发现各地类加总跟官方公布的行政区面积对不上偏差能到百分之十几。原因地理坐标系下像元面积随纬度变化高纬度地区像元实际面积比低纬度小用固定值乘会引入系统误差。解决先投影到等面积投影或高斯克吕格投影再做统计。广东常用 CGCS2000 3 度带中央经线根据研究区经度选择比如 114°E 对应 EPSG:4547 左右。投影后像元面积恒定统计结果才可靠。4.2 分类编码对不上说明书现象按说明书里的编码去提取水体结果提取出来的是大片山区。原因不同生产批次或不同来源的数据编码可能调整过说明书没同步更新。解决不要盲信文档先做目视验证。在 QGIS 里加载栅格叠加天地图或影像底图找几个已知地物点比如广州塔附近应该是建设用地珠江应该是水体查看像元值反推编码含义。确认后再批量处理。4.3 分幅数据接边处出现裂缝现象把各地级市分幅数据镶嵌后发现市界处有一条明显的类别突变线。原因分幅生产时各幅独立分类接边处分类结果不一致或者镶嵌时没有做羽化处理。解决如果数据是分幅的优先找全省整幅版本。如果没有镶嵌后用众数滤波做一下后处理或者在接受范围内忽略接边处的少量误差。做变化检测时尤其要注意接边处的假变化会干扰结果。4.4 NoData 值被当成有效类别统计现象统计结果里多出一个面积巨大的未知类别或者面积加总远超全省面积。原因NoData 值通常是 0 或 255没有排除被当成一个地类参与了统计。解决统计前先确认 NoData 值用data ! nodata做掩膜。如果 NoData 是 0 而 0 又恰好是某个地类的编码那就得回去找生产方确认或者用栅格边缘的 0 值区域判断——通常 NoData 会出现在影像边界外。4.5 重采样方法选错导致类别污染现象把 10 米数据重采样到 30 米后出现了一些原始数据里没有的类别值。原因用了双线性或立方卷积等连续型重采样方法把类别值当连续值插值了。解决分类栅格的重采样必须用最近邻法nearest neighbor。gdalwarp里加-r nearPython 里用rasterio的Resampling.nearest。如果确实需要降分辨率也可以先用众数滤波再抽样效果比直接最近邻更好。5. 进阶技巧用土地覆盖数据做生态质量评价的完整链路5.1 从土地覆盖到生态指数土地覆盖数据本身是基础底图真正体现价值的是用它算出衍生指标。一个经典应用是计算区域生态质量指数比如基于土地利用的景观生态风险指数或者生境质量模型。这里给一个简化但可复现的链路用土地覆盖计算香农多样性指数再结合不透水面比例做生态质量分级。import rasterio import numpy as np from scipy.stats import entropy with rasterio.open(gz_landcover_reclass.tif) as src: data src.read(1) profile src.profile # 滑动窗口计算香农多样性指数 window_size 50 # 500米窗口 rows, cols data.shape shannon np.zeros_like(data, dtypenp.float32) for i in range(0, rows, window_size): for j in range(0, cols, window_size): block data[i:iwindow_size, j:jwindow_size] if block.size 0: continue # 统计各类别出现频率 unique, counts np.unique(block, return_countsTrue) # 排除 NoData mask unique ! 0 if np.sum(mask) 0: continue probs counts[mask] / np.sum(counts[mask]) shannon[i:iwindow_size, j:jwindow_size] entropy(probs) profile.update(dtyperasterio.float32, nodata-9999) with rasterio.open(gz_shannon.tif, w, **profile) as dst: dst.write(shannon, 1)这段代码用 50×50 像元的滑动窗口计算香农多样性指数窗口约 500 米见方。entropy函数来自 scipy输入是各类别的概率分布。指数越高说明窗口内土地覆盖类型越丰富生态多样性越好。计算完后可以按自然断点法分成高、中、低几档再和不透水面比例做叠加识别出「高多样性但高不透水面」的冲突区域这些往往是城市扩张前沿值得重点关注。5.2 验证分类精度的一个土办法如果你手头没有独立的验证样本但又想对数据质量有个基本判断可以用一个土办法找几幅同时期的高分影像截图在土地覆盖栅格上随机撒点目视比对。具体操作是在 QGIS 里加载土地覆盖和影像底图创建 100 个随机点逐个查看点所在位置的影像地物和栅格类别是否一致。如果一致率低于 80%说明这份数据在你的研究区可能不太靠谱需要谨慎使用或者考虑自己重新分类。这个方法虽然粗糙但比盲目信任数据强。我一般会在正式分析前花半小时做这个抽查翻车过几次之后就养成了习惯。5.3 数据融合的边界最后说一个容易上头的地方有人会想把这份 10 米土地覆盖和夜间灯光、POI、路网数据全部融合起来做城市边界提取。思路没错但要注意尺度匹配。夜间灯光数据常见的是 500 米分辨率POI 是点数据路网是矢量线直接和 10 米栅格叠加会产生大量空值和尺度不一致问题。常见做法是先把所有数据统一到同一个格网比如 100 米或 250 米再做多要素加权。不要试图在 10 米尺度上融合所有数据那样计算量大不说精度提升也有限。从那以后我每次做多源融合之前都强制先跑一遍尺度一致性检查确认所有数据的空间支撑匹配了再往下走。希望帮到你。本文还有配套的精品资源点击获取