简介2020年10m精度浙江省土地覆盖土地利用数据基于哨兵影像与深度学习制作共划分耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪/冰等十类地物。数据已由原始墨卡托坐标系重投影为WGS84地理坐标系并依据最新省市行政边界裁剪形成浙江省下辖11个市的独立TIF栅格文件。压缩包共77个文件总大小约40.16MB除TIF栅格外还包含匹配的tfw坐标参考、vat.dbf属性表、cpg编码文件、aux.xml元数据以及各市的土地覆盖类型xlsx说明表和png示意图片目录按城市分类便于按需调用。已有547人浏览学习。这套数据可直接用于ArcGIS、QGIS等软件的土地利用分类统计、城市扩张分析、生态环境评估或相关课程教学免去了自行拼接、投影和裁剪的预处理工作尤其适合需要高分辨率土地利用底图的GIS从业者与科研人员。1. 2020年10m浙江省土地覆盖数据包先分清“覆盖”和“利用”再动手拿到一个名为“2020年10m精度浙江省土地覆盖土地利用.rar”的压缩包第一步不是急着解压而是先确认一件事里面到底是“土地覆盖”还是“土地利用”。这两个词在遥感产品里经常被混用但语义和用途完全不同——“土地覆盖”回答的是地表是什么水体、林地、不透水面由卫星影像光谱解译得到“土地利用”回答的是人在上面做什么水田、果园、住宅用地需要结合社会经济数据和权属信息。10米分辨率这个精度意味着数据大概率来自哨兵-2影像的地表覆盖分类产品直接按土地利用去用会出大问题比如湖州的水田在覆盖分类里往往只被标成“作物”而非“水田”。这篇文章就把一个省级10米覆盖栅格从解压、检查、重分类、面积统计到变化检测的完整链路讲清楚适合要做国土空间分析、生态评估或数据加工的工程师。2. 数据包体检与预处理从.rar到可分析的GeoTIFF2.1 解压与文件组织先列目录再解压.rar是商业压缩格式Linux 服务器默认没有unrar命令Windows 下 WinRAR 虽常用但对自动化不友好。我一般先用 7-Zip 列目录确认包内结构再决定解压路径避免一次性解出几十个 GB 的文件后才发现组织混乱。# 列出压缩包内容不实际解压 7z l 2020年10m精度浙江省土地覆盖土地利用.rar # 解压到 data 目录-o 后不要加空格 7z x 2020年10m精度浙江省土地覆盖土地利用.rar -o./data用7z l先看包内是单幅完整栅格、分幅瓦片还是带shp的工程目录。分幅栅格后期需要合并gdal_merge.py或gdalbuildvrt带lyr/qml配色文件的包通常还包含样式信息对后续出图有用。解压后先留意文件名编码包内中文文件名在部分 Linux 环境解压会乱码这是因为压缩时使用了 GBK 编码而系统默认 UTF-8。7-Zip 在部分版本能自动转换乱码时用convmv -f GBK -t UTF-8批量修正文件名不要手动一个个改。2.2 用 GDAL/Rasterio 核对元数据CRS、分辨率、nodata拿到tif后第一件事是读元数据而不是直接扔进 ArcGIS 看图层。重点核对四件事投影坐标系、像素分辨率、nodata 值、波段数。用 Python 的rasterio几行就能完成。import rasterio with rasterio.open(data/zhejiang_landcover_2020_10m.tif) as src: print(CRS:, src.crs) # 坐标系 print(分辨率:, src.res) # 单位米或度 print(范围:, src.bounds) # 左上右下坐标 print(nodata:, src.nodata) # 无效值 print(波段数:, src.count) # 分类栅格通常为1波段注意src.res返回两个值(x, y)。如果 CRS 显示为EPSG:4326WGS84 地理坐标分辨率不会是10, 10而会接近0.000089, 0.000089度这意味着该栅格是地理坐标系下的 10 米等效分辨率后续算面积前必须重投影具体方案在第 4 章讲。还要确认 nodata有的产品把 0 当作背景有的用 255有的没有设置后续统计面积时漏掉 nodata 会凭空多出成百万个像元这个细节最容易被忽略。2.3 建立分类映射表从 DN 值到地类名称10 米分辨率土地覆盖产品在栅格里存的是整数编码DN 值每个数字对应一个地类。常见产品如 Esri 全球土地覆盖、FROM-GLC 等各自有独立的编码体系拿到数据后先读取像元值的分布再确认是否与包内文档一致。import numpy as np with rasterio.open(data/zhejiang_landcover_2020_10m.tif) as src: data src.read(1) # 统计各类别像元数量 unique, counts np.unique(data, return_countsTrue) for value, count in zip(unique, counts): print(fDN{value}: {count} 像元)统计结果里如果出现文档中不存在的 DN 值说明数据经过了后处理或与文档版本不一致以实际直方图为准。常见的 10 米覆盖分类通常包含水体、树木、草地、作物、不透水面、裸地等一级类有的体系会多出“被洪水淹没的植被”“冰雪”“灌丛”等类别。建议先建一个 Python 字典做 DN 值到中文名称的映射后续所有统计与出图都基于这个字典不要反复手写数字判断减少低级错误。3. 分类体系对齐10m 土地覆盖如何落到土地利用分类3.1 土地覆盖与土地利用的语义差异这一步是整篇文章的核心也是最容易被资深工程师跳过的地方。土地覆盖分类是“所见即所得”一张 2020 年 10 米影像里杭嘉湖平原的水稻田和旱地如果光谱特征相似模型会统一判为“作物”但土地利用分类里它们是“水田”和“旱地”用途和价值差别巨大。反过来一个被林木覆盖的坟地或公园在覆盖分类里是“林地”但在土地利用分类里可能是“特殊用地”或“公园绿地”。10 米产品只靠哨兵-2 的光谱与时间序列信号能稳定区分的是常绿/落叶林、草本、水体、不透水面、裸地、作物。凡是依赖地块形状、权属、社会经济属性的类别——水田、果园、设施农用地、农村宅基地——都无法直接从像素上解译出来。所以当你把这张图用于业务现状分析时必须接受一个现实图能给出的答案是“哪里有树、哪里有水、哪里被硬化”而不是“这块地是水田还是水浇地”。3.2 覆盖分类到国标一级类的映射与风险如果业务要求输出符合《土地利用现状分类》GB/T 21010-2017一级类的结果常见的做法是按照对应用途做一次“粗粒度映射”。这里的关键原则是能 1:1 映射的类直接映射无法 1:1 的先归入最接近的一级类并在报告中明确标注不确定性。覆盖产品 DN 值覆盖类别映射到土地利用一级类映射风险1水体水域及水利设施用地低注意坑塘与河流的区别2树木林地中园地果园、茶园会被误并3草地草地中城市公园绿地会混入4被洪水淹没的植被湿地/水域按当地规程高浙江沿海滩涂归属需人工判断5作物耕地高无法拆分水田/旱地/水浇地6不透水面/建成区建设用地低但农村居民点与城镇工矿无法区分7裸地其他土地低浙江的典型问题是茶园和果园10 米光谱上茶园往往呈现为灌木/树木混合信号分类器极易把它并入“树木”类映射后归入林地但国土调查中茶园属于园地。另一个是“被洪水淹没的植被”类钱塘江口和象山港的滩涂、互花米草在模型里经常落在这个类映射到湿地还是水域缺乏统一标准。这种映射的不确定性不是 bug而是数据本身语义上限决定的。3.3 用面积构成表快速校验数据是否合理映射关系建好后先做一次全省面积构成统计用常识判断数据是否合理。浙江的地类结构大致特征林地占比最高约 60% 上下、耕地次之约 15%-20%、建设用地和水域各占一定比例、草地占比很小。如果统计结果里草地占比超过 10%通常说明分类器把低矮灌木或稀疏林地归入了草地需要调整映射或做后处理。import pandas as pd pixel_area 10 * 10 # 假设已重投影为10米分辨率 area_ha {cls: count * pixel_area / 10000 for cls, count in zip(unique, counts)} df pd.DataFrame(list(area_ha.items()), columns[DN, 面积(公顷)]) df[地类] df[DN].map(DN2CN) print(df.sort_values(面积(公顷), ascendingFalse)).rar数据包里的单位通常没有明说面积统计前必须确认像素分辨率。40 米分辨率算出来是 1600 平方米/像元10 米是 100 平方米/像元差一个数量级。浙江全省陆地面积约 10.55 万平方公里含海岛统计结果与这个量级相差超过 10% 时优先排查投影与 nodata 问题再怀疑数据本身。4. 省级栅格数据的面积统计与变化检测4.1 重投影到等积投影算面积前的必要前提在EPSG:4326地理坐标系下直接统计面积是错误做法因为经度方向长度随纬度变化浙江跨越约 2.5 个纬度纬度越高面积偏差越大。省级尺度的面积统计需要把数据投影到等积投影Albers 或等面积圆锥投影。全国常用的 Albers 参数是中央经线 105°E双标准纬线 25°N 和 47°N。gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs \ -tr 10 10 -r near -overwrite \ data/zhejiang_landcover_2020_10m.tif data/zj_lc2020_albers.tif参数说明-t_srs指定目标坐标系这里用 proj4 字符串定义 Albers 等积投影不依赖 EPSG 编码任何 GDAL 版本都能识别-tr 10 10强制输出像素为 10 米防止重投影过程改变分辨率-r near用最近邻重采样——分类栅格是离散值bilinear或cubic会在类别边界产生不存在的插值这是新手最容易犯的错误-overwrite允许覆盖已存在的输出文件。4.2 分县统计zonal_stats 的操作与关键参数省级分析经常要落到县级比如统计每个县的建设用地面积、林地覆盖率。最常见的做法是用geopandas读入县级行政区划矢量配合rasterstats做分区统计。注意要选用与业务口径一致的县域边界并检查边界坐标系与栅格是否一致。import geopandas as gpd from rasterstats import zonal_stats counties gpd.read_file(data/zhejiang_county.shp) counties counties.to_crs(EPSG:32650) # 或转成与栅格一致 stats zonal_stats( counties, data/zj_lc2020_albers.tif, categoricalTrue, # 输出每个DN值的像元计数 nodata0, # 与栅格 nodata 保持一致 all_touchedFalse, # 边界像元处理方式 geojson_outTrue # 把属性合并到 GeoJSON 中 ) result gpd.GeoDataFrame.from_features(stats) result[建设用地_公顷] result[6] * 100 / 10000categoricalTrue时各列名是 DN 值的字符串如1、2列值是该县内该类别的像元计数all_touchedFalse表示只有像元中心落入县域才计入统计更精确但可能遗漏窄长县域边界的零散像元all_touchedTrue则会把所有与边界相交的像元计入面积会偏大对沿海岛屿多、边界破碎的县域两者差异可能达到 2%-5%nodata0必须显式声明否则背景像元会被统计进某个类别。输出结果里的几何对象来自geojson_outTrue直接就是可直接制图的GeoDataFrame。4.3 与历史数据做变化矩阵先对齐再比较做 2015 年到 2020 年的土地利用变化检测时直接对两个栅格做a ! b是错误做法。两个数据源的生产机构不同投影、分辨率、分类体系都可能不一致必须先重投影到同一坐标系、重采样到同一分辨率、统一 nodata 定义。import numpy as np import rasterio with rasterio.open(data/lc_2015_albers.tif) as a_src: a a_src.read(1) with rasterio.open(data/lc_2020_albers.tif) as b_src: b b_src.read(1) # 只统计两期均为有效值的像元 valid (a 0) (b 0) (a ! 255) (b ! 255) trans np.zeros((8, 8), dtypenp.int64) np.add.at(trans, (a[valid] - 1, b[valid] - 1), 1) # 变化面积像元数 x 100 平方米转公顷 change_area np.count_nonzero((a ! b) valid) * 100 / 10000 print(变化总面积(公顷):, change_area) print(变化矩阵(行2015, 列2020):) print(trans)用np.add.at而不是trans[a, b] 1是因为直接索引加法遇到重复坐标时会覆盖丢失计数np.add.at能正确累加。变化矩阵的解读重点是“转换对”而非“变化总量”比如“耕地→建设用地”才是城市扩张的真实信号“林地→草地”可能是采伐后还未恢复也可能是分类噪声。两张 10 米产品的时间序列做变化检测单个像素的变化里混有一定比例的伪变化一般建议先做 3x3 众数滤波再比较能明显压低椒盐噪声造成的假变化。5. 直接能用的三个落地技巧重编码、出图与精度抽样5.1 类别重编码处理浙江沿海滩涂的特殊归属浙江的“被洪水淹没的植被”类集中在慈溪、上虞、象山港等地的潮间带滩涂分布着大片互花米草与海三棱藨草。在做生态评价时这一类的归属常常引发争议并入“水域”会低估湿地面积并入“草地”又会失真。我一般按业务目标处理做湿地保护评估就单独保留一个“滩涂/滨海湿地”类不并入任何国标一级类做国土空间现状分析则参考最新国土调查工作分类中湿地的定义把它与内陆滩涂归入湿地大类。无论选择哪种关键是重编码操作要可逆——保留一份原始 DN 值副本在副本上做重编码不要覆盖原始栅格。# 将 DN4(被洪水淹没的植被) 重编码为 51(自定义滩涂湿地类) reclassified data.copy() reclassified[reclassified 4] 51 np.save(data/zj_lc2020_reclass.npy, reclassified)5.2 离散色带出图土地覆盖图必须用离散色带连续渐变色会让人误读类别边界。匹配浙江实际地物特征的配色一套可直接用于matplotlib输出。import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap, BoundaryNorm colors [#a6cee3, #1f8d49, #c5e0b4, #8dd3c7, #ffff99, #d95f02, #bdbdbd, #ffffff] cmap ListedColormap(colors) norm BoundaryNorm([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5], cmap.N) fig, ax plt.subplots(figsize(10, 12)) im ax.imshow(data, cmapcmap, normnorm) cbar fig.colorbar(im, ticksrange(1, 9)) cbar.ax.set_yticklabels([水体, 树木, 草地, 洪泛植被, 作物, 不透水面, 裸地, 冰雪]) ax.set_axis_off() fig.savefig(zj_lc2020.png, dpi150, bbox_inchestight)5.3 精度抽样与置信区间估算最终业务报告里只放一张面积统计表是不够的。随机抽取 200-300 个样本点叠加同期高分影像如资源三号或天地图影像人工判读可以快速估算分类精度。当抽样点满足随机独立时用正态近似给出总体精度的 95% 置信区间公式是p ± 1.96 * sqrt(p*(1-p)/n)。如果 300 个样本里判对 255 个精度 85%置信区间约为 80.9%-89.1%。这个区间如果跨过了你业务的接受阈值就需要补样本或换生产过程。一套完整的上游到下游的处理链路做到这里数据包里的栅格才真正变成了可以用来写结论的成果。本文还有配套的精品资源点击获取