企业官网建设流程全解析

简介：本资源为2020年河南省10米精度土地覆盖与土地利用数据包，面向地理信息、遥感、城市规划及生态环境研究等领域的科研人员与学生，可解决省级尺度高分辨率地表覆被分析、土地利用变化监测等场景下的数据获取与预处理问题。压缩包共126个文件，约65.7MB，以18个tif栅格影像为核心，配套tfw坐标文件、dbf属性表、cpg编码文件、xml元数据、xlsx统计表与png预览图，覆盖河南省各地级市，便于直接加载与属性查询。数据基于10米哨兵影像，采用深度学习方法制作，分为耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪冰等类别，并已由墨卡托投影转为WGS84地理坐标系，按最新省市级行政边界裁剪，省去了自行重投影与裁剪的繁琐步骤。目前已有353人学习下载，适合需要快速开展河南土地利用制图、空间统计与变化分析的研究者使用。

1. 拿到一份 2020 年 10m 分辨率河南省土地覆盖数据，先别急着打开

你手上如果有一份名为「2020年10m精度河南省土地覆盖土地利用.rar」的压缩包，大概率是从公开地理数据平台、科研数据共享站点或者同行手里流转过来的。它本质上是一套栅格数据产品：覆盖河南省全域，空间分辨率 10 米，时间截面是 2020 年，每个像元记录一个土地覆盖/土地利用类别码。10 米意味着什么？意味着一个像元在地面上是 10m×10m，一平方公里有 1 万个像元，河南约 16.7 万平方公里，粗算下来是十几亿个像元。这个量级决定了它不可能用 Excel 打开，也不适合直接扔进 QGIS 里靠鼠标点。

这类数据真正解决的问题是：你需要知道 2020 年河南每一块地是耕地、林地、草地、水体、建设用地还是未利用地，而且要精确到地块级别。做国土空间规划、耕地非农化监测、城市扩张分析、生态红线评估、遥感变化检测的人，都会盯上这种 10m 级别的省级产品。它比 30m 的全球产品细，比 1m 的商业影像便宜甚至免费，是性价比很高的中间档。但拿到 .rar 只是起点，真正的坑在于：解压之后你看到的可能是一堆 .tif、一个 .prj、一个说明文档，甚至还有 .img 或 .dat，坐标系、类别编码、无效值全都不一定写在显眼处。这篇就按我实际处理这类数据的顺序，把解压、校验、裁剪、统计、踩坑一条线讲清楚。

2. 解压之后先做三件事：看结构、验坐标系、读类别码

2.1 用命令行快速摸清压缩包里的文件构成

很多人习惯双击解压，然后被一堆同名不同后缀的文件搞晕。我一般先在 Linux 或 WSL 下用unzip -l看一眼列表，再决定解压到哪个目录。Windows 上如果没有 unzip，用 7-Zip 的命令行版本7z l也一样。

# 先列出压缩包内容，不实际解压，确认文件数量和类型 unzip -l "2020年10m精度河南省土地覆盖土地利用.rar" # 如果确认是分幅或分县存放，解压到独立目录，避免污染当前工作区 mkdir -p henan_lulc_2020 unzip "2020年10m精度河南省土地覆盖土地利用.rar" -d henan_lulc_2020 # 进入目录后按扩展名统计，快速判断数据组织方式 cd henan_lulc_2020 find . -type f | sed 's/.*\.//' | sort | uniq -c | sort -nr

这段命令的逻辑很直接：第一步只读列表，防止压缩包里有大量小文件直接铺满当前目录；第二步指定-d解压到独立文件夹；第三步用find加sed提取扩展名并计数。如果输出里.tif占绝大多数，说明是标准栅格分幅；如果出现.img，那是 ERDAS IMAGINE 格式，GDAL 也能读，但要注意可能附带.rrd金字塔文件；如果出现.dat加.hdr，那是 ENVI 格式，需要成对保留。参数上唯一要留意的是-d后面不要带中文空格，路径里如果有中文，建议先改成英文目录再操作，否则某些老版本 GDAL 会报编码错误。

2.2 用 gdalinfo 确认坐标系、分辨率和无效值

解压完别急着写代码，先用gdalinfo把第一个 .tif 的元数据打出来。这一步能直接告诉你三件关键事：投影是什么、像元大小是不是真的 10m、NoData 值是多少。

# 查看单个栅格的完整元数据，重点关注 Coordinate System 和 NoData Value gdalinfo henan_lulc_2020/henan_2020_lulc.tif # 如果文件很多，批量提取坐标系和分辨率，输出成表格方便比对 for f in henan_lulc_2020/*.tif; do echo "=== $f ===" gdalinfo "$f" | grep -E "Pixel Size|Coordinate System|NoData Value" done

gdalinfo的输出里，Pixel Size如果显示(10, 10)，说明分辨率达标；如果显示(0.0001, 0.0001)且坐标系是地理坐标（WGS84），那实际地面分辨率会随纬度变化，河南大约在 10m 上下浮动，需要重投影才能得到严格等面积统计。Coordinate System常见的是WGS 84 / UTM zone 49N或50N，河南跨了 49 和 50 两个带，如果数据是分幅的，很可能一部分在 49N 一部分在 50N，合并前必须统一。NoData Value如果是0，而类别码里 0 又代表「无数据」，那统计时就要显式排除，否则会把背景值算成某一类。我见过最坑的情况是 NoData 写成-128，但说明文档里没提，结果统计出来多出一个莫名其妙的类别。

2.3 类别码对照表：没有它，栅格只是一堆数字

10m 土地覆盖产品的类别码通常沿用几套主流体系：一套是 6 大类（耕地、林地、草地、水体、建设用地、未利用地），一套是 8 到 10 类的细分体系，还有的会带二级码。压缩包里一般会有一个.txt或.pdf说明，但有时候被漏掉。如果确实没有，可以按像元值分布反推，但更稳妥的是找同源产品的公开文档。

拿到对照表后，建议在代码里写成一个字典，后续统计直接映射，不要每次靠记忆。如果说明文档里的类别码和gdalinfo看到的实际值对不上，以实际值分布为准，先做一次直方图统计再下结论。

3. 用 Python 把整省数据裁成分析单元并做面积统计

3.1 环境准备与读取：rasterio 比 gdal 更顺手

命令行看元数据可以，但真正做裁剪和统计，我习惯用 Python 的rasterio加numpy。它比直接调 GDAL 的 Python 绑定干净，读进来就是数组，配合geopandas做矢量边界裁剪很顺。

import rasterio from rasterio.mask import mask import geopandas as gpd import numpy as np # 打开整省栅格，先看基本信息 src_path = "henan_lulc_2020/henan_2020_lulc.tif" with rasterio.open(src_path) as src: print("CRS:", src.crs) print("分辨率:", src.res) print("波段数:", src.count) print("NoData:", src.nodata) # 读取第一波段为数组，方便后续统计 data = src.read(1) print("数组形状:", data.shape) print("唯一值:", np.unique(data)[:20])

这段代码先打印元数据，再把第一波段读成 numpy 数组。src.read(1)里的1是波段索引，土地覆盖产品通常只有一个波段，所以固定读 1。np.unique只打印前 20 个唯一值，防止类别太多刷屏。如果src.nodata是None，说明文件里没写无效值，需要手动根据直方图判断，比如 0 或 255 出现频率异常高，那基本就是背景值。参数上唯一要注意的是内存：整省 10m 数据读成数组后可能占几个 GB，如果机器内存不够，不要一次性read，改用分块读取，后面会讲。

3.2 按行政区裁剪：用矢量边界切出你要的市或县

整省数据直接统计只能得到全省总面积，但实际分析往往要落到市、县甚至乡镇。这时候需要一份行政边界矢量，用rasterio.mask按几何裁剪。

# 读取行政边界矢量，假设字段名是 NAME，要提取郑州市 admin = gpd.read_file("henan_admin.shp") zhengzhou = admin[admin["NAME"] == "郑州市"] # 用郑州市边界裁剪栅格 with rasterio.open(src_path) as src: # 确保矢量坐标系和栅格一致，不一致先转 if zhengzhou.crs != src.crs: zhengzhou = zhengzhou.to_crs(src.crs) # mask 返回裁剪后的数组和变换参数 out_image, out_transform = mask(src, zhengzhou.geometry, 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("zhengzhou_lulc_2020.tif", "w", **out_meta) as dest: dest.write(out_image)

mask函数的第一个参数是打开的栅格数据集，第二个是几何对象列表，crop=True表示裁剪后去掉外围空白，减小文件体积。out_meta复制原元数据后必须更新height、width和transform，否则写出的文件地理参考会错位。坐标系一致性检查不能省，我遇到过矢量是 CGCS2000 而栅格是 WGS84 UTM 的情况，不转直接裁，结果边界偏了几百米。如果裁剪多个县，建议写成循环，每个县输出一个文件，命名带上行政区代码，方便后续批量统计。

3.3 面积统计：从像元数到平方公里，别忘了几何校正

裁剪完就可以统计各类别面积了。核心逻辑是：统计每个类别码的像元数，乘以单个像元面积，再换算成平方公里。如果数据是投影坐标系且分辨率严格 10m，单个像元面积就是 100 平方米。

# 读取裁剪后的栅格 with rasterio.open("zhengzhou_lulc_2020.tif") as src: data = src.read(1) nodata = src.nodata # 排除无效值 if nodata is not None: valid = data[data != nodata] else: valid = data # 统计每个类别的像元数 classes, counts = np.unique(valid, return_counts=True) # 单个像元面积，投影坐标系下 10m 分辨率即 100 平方米 pixel_area = 100 # 平方米 # 换算成平方公里 for cls, cnt in zip(classes, counts): area_km2 = cnt * pixel_area / 1e6 print(f"类别 {cls}: {cnt} 像元, 面积 {area_km2:.2f} 平方公里")

np.unique配合return_counts=True直接返回类别和对应像元数，比循环快得多。pixel_area这里写死 100 是因为确认了投影坐标系下分辨率为 10m；如果数据是地理坐标系，这个值会随纬度变化，必须先用gdalwarp重投影到等面积投影或 UTM，再统计。1e6是平方米到平方公里的换算系数。如果要做多年对比，建议把结果存成 CSV，字段包括行政区、年份、类别、面积，后续直接读表做变化矩阵。我一般还会加一列百分比，方便快速看结构。

4. 避坑与排查：10m 省级土地覆盖数据最常见的五个翻车点

4.1 现象：统计出来总面积比河南实际面积大很多

原因通常有两个：一是没排除 NoData，背景值被算进了某一类；二是数据分幅之间有重叠，合并时没有去重。解决方法是先确认src.nodata，如果为None，用直方图看哪个值出现频率异常高，手动设为无效值；分幅数据用gdal_merge.py合并时加-n 0指定无效值，或者用gdalbuildvrt先建虚拟镶嵌再转 tif，避免重复像元被累加。

4.2 现象：裁剪后的栅格边界和矢量对不上，整体偏移

原因几乎都是坐标系不一致。栅格是 UTM 49N，矢量是 WGS84 地理坐标，直接裁就会偏。解决方法是统一坐标系，用to_crs把矢量转到栅格坐标系，或者用gdalwarp把栅格转到矢量坐标系。注意不要反复转，转一次就固定下来，反复转会引入重采样误差。

4.3 现象：类别码和说明文档对不上，多出几个没见过的值

原因可能是数据经过了二次处理，或者说明文档对应的是另一版产品。解决方法是先做全图唯一值统计，把实际出现的类别列出来，再对照公开文档反推。如果多出的值集中在边界区域，可能是混合像元或边缘效应，可以考虑用众数滤波平滑，但不要过度处理，否则会损失真实细节。

4.4 现象：整省数据读进内存直接爆掉

原因是一次性read了十几亿像元。解决方法是分块读取，用rasterio的block_windows或者window参数，每次只读一块，统计完累加。也可以先用gdal_translate做降采样，生成低分辨率概览用于快速预览，确认无误后再对目标区域做全分辨率处理。

4.5 现象：面积统计结果和官方公布数据差几个百分点

原因可能是分辨率导致的边界混合像元、类别定义差异、或者统计时包含了水域中的岛屿等细节。解决方法是对比时明确口径，比如官方耕地面积是否包含临时种植园，建设用地是否包含农村道路。如果差异在 5% 以内，通常属于正常范围；如果超过 10%，要检查投影和无效值处理。我一般会保留一份统计脚本和中间结果，方便回溯每一步。

5. 进阶技巧：用分块统计加变化矩阵，把 10m 数据用出 30m 没有的细节

5.1 分块统计：内存不够时的标准操作

整省 10m 数据在普通笔记本上确实吃力，但分块统计可以解决。核心思路是按行块或窗口遍历，每块统计类别像元数，最后累加。

from rasterio.windows import Window def block_statistics(tif_path, block_size=2048): """分块统计各类别像元数，适合大文件""" counts = {} with rasterio.open(tif_path) as src: nodata = src.nodata height, width = src.height, src.width # 按块遍历，block_size 控制每次读取的行数 for row in range(0, height, block_size): for col in range(0, width, block_size): # 计算当前窗口，防止越界 win = Window(col, row, min(block_size, width - col), min(block_size, height - row)) data = src.read(1, window=win) if nodata is not None: data = data[data != nodata] classes, cnts = np.unique(data, return_counts=True) for c, n in zip(classes, cnts): counts[c] = counts.get(c, 0) + int(n) return counts # 调用并换算面积 result = block_statistics("henan_lulc_2020/henan_2020_lulc.tif") for cls, cnt in sorted(result.items()): print(f"类别 {cls}: {cnt * 100 / 1e6:.2f} 平方公里")

block_size默认 2048，意味着每次读 2048 行，内存占用可控。Window的四个参数是列起点、行起点、宽度、高度，用min防止最后一块越界。统计结果累加到字典里，最后统一换算。这个方法比一次性读整图慢一点，但内存占用从几个 GB 降到几十 MB，普通办公机也能跑。如果要做多年变化矩阵，把每年的统计结果存成字典，再两两对比，就能得到耕地转建设用地、林地转耕地等转移面积。

5.2 变化矩阵：两期数据叠加，看谁在变

有了分块统计的基础，做变化矩阵就是多一步叠加。假设你有 2020 和 2015 两期同分辨率数据，先确保坐标系和范围一致，然后逐像元比较。

# 假设两期数据已经对齐，读取同一窗口 with rasterio.open("henan_2015_lulc.tif") as src15, \ rasterio.open("henan_2020_lulc.tif") as src20: # 读取同一区域，这里以整图为例，实际建议分块 old = src15.read(1) new = src20.read(1) # 排除无效值 mask_valid = (old != src15.nodata) & (new != src20.nodata) old_valid = old[mask_valid] new_valid = new[mask_valid] # 组合成转移对，统计频次 transition = np.zeros((10, 10), dtype=np.int64) for o, n in zip(old_valid, new_valid): transition[o, n] += 1 # 打印主要转移方向 for i in range(10): for j in range(10): if transition[i, j] > 0 and i != j: print(f"{i} -> {j}: {transition[i, j] * 100 / 1e6:.2f} 平方公里")

transition矩阵的行是旧类别，列是新类别，对角线是未变化部分。mask_valid确保两期都有效的像元才参与统计。循环部分如果数据量大可以用numpy的bincount加速，但为了逻辑清晰，这里用双重循环展示。实际跑的时候，建议先分块再叠加，否则两期整图同时读进内存容易爆。变化矩阵最有价值的地方是能看出「耕地转建设用地」这种关键转移，10m 分辨率下，连农村宅基地的扩张都能捕捉到，这是 30m 数据做不到的。

5.3 一个我常犯的错：忘了记录处理链

最后说个血泪经验。我早期处理这类数据时，裁剪、重投影、统计各写一个脚本，跑完就扔，结果三个月后要复现某个市的面积，发现忘了当时用的哪版矢量边界，也忘了有没有做众数滤波。后来我养成了一个习惯：每个项目目录下放一个README.md，记录数据来源、处理步骤、关键参数、每一步的输出文件名。哪怕只是临时分析，也写三行。这个习惯帮我省了无数次后悔药。10m 省级数据量不小，处理一次不容易，把处理链记清楚，比多跑几次脚本重要得多。希望帮到你。

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