企业官网建设流程全解析

简介：2019年10m精度云南省土地覆盖土地利用数据包，基于哨兵影像与深度学习制作，分类涵盖耕地、林地、草地、灌木、湿地、水体、建筑用地、裸地及雪/冰等十类。资源已完成坐标系转换与行政边界裁剪，统一为WGS84地理坐标系，并按云南省各地市分别输出，便于研究或应用中直接调用。压缩包共112个文件，包含16个TIF栅格文件及配套的tfw定位文件、xml元数据、dbf/cpg属性表、png预览图与xlsx表格，合计约115.71MB，结构清晰，适合城乡规划、生态环境及GIS学习者作为底图数据或训练样本。目前已有327人学习下载，可满足对省级尺度高精度土地覆盖数据的快速获取需求。

1. 拿到云南10m土地覆盖数据之后，先别急着出图

2019年10m精度云南省土地覆盖土地利用.rar这类压缩包，近年在GIS从业者手里流转得很多。它本质上是全球10米分辨率土地覆盖产品按云南省边界裁剪出来的本地化数据，常见来源是欧空局或ESRI基于Sentinel-2制作的2019年分类结果。很多人下载完第一件事就是拖进ArcGIS或QGIS拉伸显示，结果要么颜色乱成调色盘，要么把水体显示成林地，然后开始怀疑数据坏了。这里想先立一个判断：10m精度指的是像元大小为10米，而不是每个地类都精确到10米。分类结果受训练样本、云覆盖和地形阴影影响，局部误差可能超过一个像元。它适合做宏观概览、面积统计、变化趋势，不适合当作地块级取证依据。

我一般拿到这类数据会先花半小时做四件事：查分类体系、验投影、算有效覆盖、做一张能稳定出图的样式。这套流程能避掉后续大部分坑。下面从一个rar包的落地路径展开，从数据底细讲到面积统计，再落到云南地形带来的各种翻车点，最后给一个用NDVI做自检的土办法。新手可以照着命令走，熟手可以直接跳去第5章看踩坑清单。

2. 数据底细与选型：2019年10m土地覆盖数据到底是什么

2.1 分类体系与栅格文件组织

先看压缩包内部的结构。常见的10米土地覆盖产品是一个GeoTIFF加一个样式文件（.lyr或.qml），有的还带一个CSV格式的分类描述。云南省的文件名里通常有Yunnan、2019、LC等字样，栅格值从0到11，或从1到12，不同产品对应关系完全不同。这里要特别警惕：ESRI产品中1是水体、2是树木、3是草地；但ESA WorldCover中10是树木、20是灌木。拿ESRI的图例去渲染ESA的产品，整张图都会是错的。

打开文件后用gdalinfo看波段数和像素类型。多数10米产品是单波段8bit分类栅格，一个像素只占一个字节，文件大小可控。云南面积约39万平方公里，10米分辨率下理论像元数接近40亿，但压缩后tif通常在1到2GB，因为分类图有大量游程压缩。如果看到的是一个三波段RGB预览图而不是单波段分类图，说明压缩包里带的是可视化版本，真正的分类图层在子目录里。过去我懒得翻子目录，直接用RGB预览图做统计，结果全部作废，这个亏吃过一次就不会再吃。

2.2 10m与30m、250m的真实权衡

选分辨率不能只看越细越好。30米Landsat数据历史长，能回溯到上世纪80年代，适合做长时间序列；250米MODIS数据时间频次高，适合做物候；而10米Sentinel-2数据从2015年后才稳定获取，2019年的产品已经是比较成熟的批次。对云南这种山高谷深、地块破碎的地形来说，30米像元跨过一条河谷时，混入的类别可能超过一半；10米像元虽然也混，但至少能把梯田、林窗、村寨边界大致分开。

但10米也带来噪声问题。山区阴影、单次观测云遮挡、同物异谱都会让分类结果出现大量椒盐噪点。做面积统计时，直接用原始10米栅格统计，会发现灌木面积忽大忽小，原因在于训练样本里灌木与草地本身就难分。实战中我常做一步：3x3众数滤波，把孤立像元去掉后再统计，结果更接近业务口径。这一步属于常见做法，不算严谨科学，但确实有效。顺便说一句，千万别用均值滤波处理分类栅格，否则会出现0.6之类的地类值，后期查都查不出来。

2.3 云覆盖与有效像元：云南数据的隐形参数

云南地处低纬高原，干季晴天多，雨季云量常年在60%以上。2019年的合成产品理论上用到全年多时相合成，但部分月份数据源缺口不小。判断数据能不能用，首先要看有效像元比例，看文件大小没用。有效像元指分类值落在合法范围内的像元，如果某个类别像元占比低于0.1%，且恰好集中在高海拔山区，那多半是云阴影或雪被误分。这类区域在成图时建议标注为“数据不确定区”，而不是强行修图。

我见过不少项目把云南西北部的雪冰类像元直接归零，理由是当地冬季确实有雪，但分类产品中的雪冰面积会随年份浮动。如果你做土地利用变化，而基准年恰好选了多雪年份，森林减少面积会被严重高估。处理办法是统计前把坡度大于35度且海拔高于4000米的像元单独拎出来看，这些地方的分类误差自带Buff，不要直接并进森林或裸地。云南这种地形，任何遥感分类产品都需要结合DEM做二次判读，这是项目开始前就应该写进技术方案里的。

3. 从rar到可用TIF：解压、校验与转坐标系实操

3.1 解压与完整性校验：用命令而不是双击

拿到一个rar包，第一步不是双击解压到桌面，而是用命令行做一次完整校验。rar文件在网盘里传几轮后，文件头损坏概率不低。Windows上可以用WinRAR自带的rar t，Linux下用unrar t。下面以Linux环境为例，因为这后面处理GDAL大多在Linux或WSL下跑。

# 校验压缩包完整性，不释放文件 unrar t "2019年10m精度云南省土地覆盖土地利用.rar" # 查看包内文件列表，避免解压出奇怪的路径 unrar lb "2019年10m精度云南省土地覆盖土地利用.rar" # 解压到data目录，保留原路径结构 unrar x "2019年10m精度云南省土地覆盖土地利用.rar" ./data/

t参数是test，只校验CRC不写盘；lb列出文件路径；x是完整解压并保留目录结构。逻辑很简单：先测后解，能省掉后期缺文件的麻烦。解压后立刻用ls -lh查看每个tif的大小，若某个tif只有几KB，那多半是空文件或损坏文件，趁早重新下载。另外注意RAR包里的文件名如果包含中文或空格，后续脚本处理前统一重命名成y2019_yunnan_lc.tif这类格式，能避免大量编码问题。这不是玄学，是Python读取中文路径时在Windows控制台下经常乱码。

3.2 检查投影与波段：gdalinfo是照妖镜

解压出的tif很可能落在Web墨卡托坐标系下（EPSG:3857），因为很多在线发布源直接吐瓦片坐标。云南东西跨度大，Web墨卡托虽然低纬变形不大，但面积统计需要等积投影。先用gdalinfo看元数据，再决定是否投影。

# 查看栅格基本结构 gdalinfo ./data/yunnan_lc2019.tif # 如果栅格太大，只看关键信息 gdalinfo -stats ./data/yunnan_lc2019.tif | head -60

看输出的关键三行：Size is后面的行列数；Coordinate System is后面的EPSG代码；NoData Value=后面的值可能是255，也可能是0。最容易翻车的就是NoData。很多产品把海洋或境外区域设为0，而0又恰好是合法分类值（水体可能为0或1），导致裁剪边界出现一圈奇怪的黑边或白边。如果发现NoData和合法值重叠，必须在下游处理里重新定义NoData。处理办法是gdal_translate -a_nodata 255，把255设为新的NoData，因为255在大多数8bit分类体系里不是合法地类。

3.3 重投影与裁剪：两个坑一次填平

云南常用的投影是Albers等积圆锥投影或UTM 47N（EPSG:32647）。如果只是做省级统计，我一般直接用UTM 47N，它覆盖云南大部分区域，像元面积接近常数，统计方便。下面的命令先把数据从3857重投影到32647，再用省界矢量裁剪。注意这里用-r near而不是cubic，因为分类栅格不允许插值，三次卷积会把地类值变成小数。

# 重投影到UTM 47N，最近邻重采样 gdalwarp -overwrite \ -t_srs EPSG:32647 \ -r near \ -dstnodata 255 \ ./data/yunnan_lc2019.tif ./data/yunnan_lc2019_utm47.tif # 用云南省界矢量裁剪，同时保持投影不变 gdalwarp -overwrite \ -cutline ./data/yunnan_boundary.shp \ -crop_to_cutline \ -dstnodata 255 \ -of GTiff \ ./data/yunnan_lc2019_utm47.tif ./data/yunnan_lc2019_clip.tif

参数说明：-r near使用最近邻重采样，分类值不会被插值，这是处理分类栅格的基本纪律。-dstnodata 255把掩膜区设为255，统计时直接剔除。-cutline指定边界矢量，-crop_to_cutline让输出范围与边界完全一致，不会留下矩形白底。如果数据本身已经是按云南裁剪好的，重投影时不要再次裁剪，否则会在边界处产生第二次重采样，导致边界像元值变化。判断是否已经裁剪，看gdalinfo里的角点坐标是否落在云南边界附近即可。

重投影后影像的像元大小可能变成10.2米或9.8米，这是地图投影变形造成的，不是错了。如果后面要做变化检测，所有年份的数据都必须重投影到同一坐标系，否则像元错位会让你后悔没有做好这一步。

3.4 为输出建金字塔与颜色表

处理完成后，建立金字塔能让你在QGIS里缩放不卡。分类栅格用平均采样会得到奇怪颜色，所以要用最邻近采样。同时写一个简单的颜色表，把地类颜色固定下来，后续所有图都用这一套颜色，省得每次调样式。

# 建立金字塔，最近邻采样 gdaladdo -r nearest ./data/yunnan_lc2019_clip.tif 2 4 8 16 # 用文本颜色表直接生成带颜色的渲染tif gdaldem color-relief -of GTiff \ -nearest_color_entry \ ./data/yunnan_lc2019_clip.tif \ ./data/lc_color.txt ./data/yunnan_lc2019_color.tif

颜色表lc_color.txt的格式是每行一个值加对应R G B，例如1 31 79 255表示水体蓝色。-nearest_color_entry保证像素值落在某区间时取最近的颜色，不会插出中间色。这个带颜色的tif可以拖进软件当底图，但项目交付时一定要保留一份未做渲染的原始分类tif。我见过有人只存了color tif，结果业务方一问“你这里的类4是什么”，他只能对着颜色猜，非常被动。

4. 分类图转矢量与面积统计：让土地覆盖数据变成业务口径

4.1 像元面积统计：不要直接数像元数

统计面积时很多人直接用像元数乘以100平方米，这在局部小范围内成立，但全省范围会累积投影变形误差。更稳的做法是读取栅格后统计每类像元数，再乘以像元实际面积。UTM 47N下像元面积随纬度和位置变化很小，但严谨起见，我们可以用栅格的地理变换参数计算面积。

from osgeo import gdal import numpy as np ds = gdal.Open('./data/yunnan_lc2019_clip.tif') band = ds.GetRasterBand(1) # 读取分类栅格，把NoData之外的像元作为有效区域 clc = band.ReadAsArray().astype(np.uint8) nodata = band.GetNoDataValue() valid_mask = clc != nodata # 统计各类像元数 classes, counts = np.unique(clc[valid_mask], return_counts=True) # 获取像元尺寸（UTM下单位为米） gt = ds.GetGeoTransform() pixel_area = abs(gt[1] * gt[5]) # 输出每个类别的面积（平方公里） for cls, cnt in zip(classes, counts): area_km2 = cnt * pixel_area / 1e6 print(f'类 {cls}: {cnt} 像元, {area_km2:.2f} km²')

参数说明：gt[1]是东西方向像元宽，gt[5]是南北方向像元高（投影坐标系下通常为负），乘积绝对值就是单像元面积。这个代码的隐含假设是投影后像元是规则矩形，UTM下近似成立。统计前先剔除NoData，否则边界外的黑色区域会被计入类0或类255。如果你发现某类面积大得离谱，先回头检查NoData值，而不是怀疑算法。

4.2 栅格转矢量：设置聚合参数避免碎面

业务方经常要shp，不要tif。但直接把分类栅格转矢量，会产生密密麻麻的碎面，云南这种陡峭地带更严重，一个10米像元的变化就是一个多边形。常见做法是先做众数滤波，再转矢量，最后按面积过滤碎面。

# 先做3x3众数滤波，去掉孤立像元 gdal_fillnodata.py -md 3 -si 1 ./data/yunnan_lc2019_clip.tif ./data/yunnan_lc2019_fill.tif # 使用gdal_polygonize.py转矢量 gdal_polygonize.py \ ./data/yunnan_lc2019_fill.tif \ -f "ESRI Shapefile" \ ./data/yunnan_lc2019_poly.shp \ ./data/yunnan_lc2019_poly layer1

这里有个隐藏坑：gdal_polygonize.py输出的shp字段只有DN，随后用ogr2ogr按面积过滤时，如果shp没有定义投影坐标系，面积字段可能是经纬度平方，不是平方米。所以更稳的做法是转出后，先定义投影再计算面积并过滤。

# 定义投影并过滤碎面，保留面积大于1公顷的图斑 ogr2ogr -overwrite \ -t_srs EPSG:32647 \ -dialect sqlite \ -sql "SELECT *, ST_Area(geometry) AS area_m2 FROM yunnan_lc2019_poly WHERE ST_Area(geometry) > 10000" \ ./data/yunnan_lc2019_poly_clean.shp \ ./data/yunnan_lc2019_poly.shp

ST_Area在投影坐标系下返回平方米；-t_srs确保输出shp自带投影信息。过滤阈值10000平方米换来的是图面整洁和更快的渲染速度。代价是丢掉了小于1公顷的独立地类图斑，如果你的项目关注小规模零散地块，就不要过滤这么狠，改到1000平方米或保留全部。这个权衡没有标准答案，取决于业务口径。

4.3 按行政区汇总：统计表长什么样

业务上最常要的统计是“云南省各州市土地利用面积”或“某流域地类构成”。用rasterstats库可以直接对矢量分区做分类统计，输出一个CSV表。

import geopandas as gpd import pandas as pd from rasterstats import zonal_stats # 读取州市界线矢量，并统一投影到UTM 47N zones = gpd.read_file('./data/yunnan_cities.shp').to_crs('EPSG:32647') # 对分类栅格做分区统计，categorical=True会统计每个类别的像元数 stats = zonal_stats( zones, './data/yunnan_lc2019_clip.tif', categorical=True, nodata=255, geojson_out=True, ) rows = [] for st in stats: props = st['properties'] row = {'市州': props['name']} for k, v in props.items(): # 过滤掉非统计字段和NoData类 if k.isdigit() and int(k) != 255: row[f'类{k}_km2'] = round(v * 100 / 1e6, 2) rows.append(row) df = pd.DataFrame(rows) df.to_csv('./data/yunnan_lc2019_by_city.csv', index=False)

categorical=True会让zonal_stats直接返回每个分类值的像元数。乘以100平方米再除以1e6得到平方公里。这段代码输出的表会非常整齐，每行一个市州，每列一个地类面积。报数据的时候，最好同时输出一个面积占比列，例如“类2占该州市总面积的比例”，这样领导才不用自己拿计算器去按。占比可以直接计算：row[f'类{k}_km2'] / total_area * 100。

另外提醒一句：行政边界矢量版本不同，统计结果会有几个百分点的差异。如果项目跨年度对比，所有年份必须使用同一版本的边界文件，否则变化量会掺入边界修订的噪声。

5. 避坑：云南地形导致的五个常见问题与排查

5.1 全黑或全白：先看NoData再看直方图

现象：把tif拖进软件，全黑，拉伸无效。原因有两个：一是NoData被设成0，0又是合法地类值，渲染时整图被当成空值；二是颜色表没有随文件加载，软件找不到渲染映射。排查顺序是先看直方图。

gdalinfo -hist ./data/yunnan_lc2019_clip.tif | grep -A 20 "Histogram"

如果直方图显示0像元占比99%，那文件本身可能下载错了。如果其他值正常，只是渲染黑，就用gdal_translate -a_nodata 255把NoData改掉，再重新加载。注意gdal_translate会重写整个文件，执行前先备份，否则改了后悔没药可吃。

5.2 高山积雪被分进水体

现象：德钦、香格里拉一带统计出的水体面积远高于常年水面。原因：2019年产品中，冰川、雪地与水体在部分合成场景中存在同物异谱，阴坡积雪被分成了水体类。处理办法不是直接改分类栅格，而是叠加DEM做掩膜。

# 用DEM计算坡度和海拔，把海拔4000米以上且坡度大于20度的水体像元重分类为'高山不确定' gdaldem slope ./data/dem_utm47.tif ./data/dem_slope.tif \ -p -s 111120 -a 10 # 用gdal_calc.py做条件赋值 gdal_calc.py -A ./data/yunnan_lc2019_clip.tif \ -B ./data/dem_utm47.tif -C ./data/dem_slope.tif \ --outfile=./data/yunnan_lc2019_masked.tif \ --calc="((A==1) * (B>4000) * (C>20)) * 200 + (A!=1) * A"

解释一下：-A是分类图，-p生成百分度坡度，-s指定水平因子。gdal_calc.py中(A==1)表示原分类为水体，(B>4000)和(C>20)表示高海拔陡坡，满足条件的像元重分类为200（你可以定义成“高山不确定”）。其他像元保留原值。加这样一层掩膜之后，水体统计就不会被积雪地带污染了。

5.3 边界出现连续不自然条带

现象：云南省界线外侧有一圈与界线平行的异常区块，颜色介于两个地类之间。原因：gdalwarp裁剪时，边界外像元被赋予默认的srcnodata 0，而0恰好是合法分类值，导致边界带被填充成错误类别。解决方法是重投影时显式指定源NoData，并在裁剪后重新检查。

# 如果已经出现条带，用缓冲矢量'收边' ogr2ogr -dialect sqlite \ -sql "SELECT ST_Buffer(geometry, -10) FROM yunnan_boundary" \ ./data/yunnan_boundary_inner.shp ./data/yunnan_boundary.shp # 用内缩边界重新裁剪 gdalwarp -overwrite -cutline ./data/yunnan_boundary_inner.shp \ -crop_to_cutline -dstnodata 255 \ ./data/yunnan_lc2019_utm47.tif ./data/yunnan_lc2019_fixed.tif

ST_Buffer(geometry, -10)把边界向内收缩10米，也就是一个像元，把受污染的边缘像元裁掉。这样做会让统计面积比实际略小，但换来的干净边界对出图更有价值。

5.4 统计面积和官方年鉴对不上

现象：用栅格统计的耕地面积，与省自然资源厅年鉴数字差距超过15%。原因：分类器把大量撂荒地、梯田、园地归并到了草地或灌木；同时年鉴口径来自国土调查，分类体系完全不同。不要试图在这份数据上去“对齐”国土调查。正确的定位是：10米覆盖产品适合看空间分布、相对占比和变化趋势，不适合做绝对面积法定统计。如果你必须给一个估算值，那就按自己的重分类表调整：例如把“草地”中的25%估计为撂荒耕地，然后单独写一个说明字段。这种做法属于模型估算，不是数据修正，千万别把栅格里的值直接改了。

5.5 图例颜色和别人发布的风格不一致

现象：加载官方样式文件后，水体是红色、森林是黑色。原因：不同发布版本的颜色板顺序不一样，文件名里都是2019，但实际对应ESA或ESRI不同批次。最稳的做法是自己生成颜色表，然后存成QGIS样式。颜色表的格式很简单：每个分类值一行，写RGB。我常用的配色是：水体31-79-255，树木10-107-53，草地168-212-143，耕地255-211-0，湿地170-170-170。保存为lc_style.qml后，整个项目组的出图效果就统一了，不再有人“突然交出一张反色图”。

5.6 一条命令快速体检

写一个几十行的脚本太累，直接用一个Python一行统计方式来做整体体检：

python3 -c " from osgeo import gdal import numpy as np ds=gdal.Open('./data/yunnan_lc2019_clip.tif') b=ds.GetRasterBand(1) a=b.ReadAsArray() vals,cnts=np.unique(a,return_counts=True) total=a.size for v,c in zip(vals,cnts): print(v, round(c/total*100,2), '%') "

如果输出结果中NoData值255占比超过5%，说明裁剪边界留白太多，或者镶嵌时源数据就存在大量空洞。如果某个类别占比接近0%，结合地理位置判断是合法还是异常。这个命令花不了几秒，但能把分类分布、NoData比例、异常值一次看清。我习惯在每次重分类后都跑一遍，形成一条快速自检的肌肉记忆。

6. 进阶：用NDVI做2019年土地覆盖的精度自检

拿到这份10m数据，如果没有实地样点，怎么判断它到底靠不靠谱？一个低成本的办法是用同一年的哨兵2号影像合成NDVI，对分类结果做交叉验证。原理很简单：不同地类的NDVI分布应该有明显区分。水体NDVI接近0或为负，森林通常高于0.6，草地和耕地集中在0.3到0.6，裸地和建设用地低于0.2。如果某个“森林”像元对应位置的NDVI只有0.15，那这个分类十有八九是错的。

操作上我不建议下载整年影像做逐像元验证，那样计算量太大。更实用的方案是分层随机抽样：每个地类别随机抽200个点，提取这些点位上的NDVI中位数和标准差，画出每类的箱线图。如果某类别的NDVI分布与经验值完全偏离，就该怀疑这类被混淆了。

import geopandas as gpd import numpy as np from osgeo import gdal from rasterio.sample import sample_gen import rasterio # 读取分类栅格，生成每个类别的随机采样点 src = rasterio.open('./data/yunnan_lc2019_clip.tif') clc = src.read(1) valid_mask = (clc != 255) # 用numpy随机采样200个位置 rows, cols = np.where(valid_mask) np.random.seed(42) idx = np.random.choice(len(rows), size=2000, replace=False) samples = [src.transform * (cols[i], rows[i]) for i in idx] class_vals = clc[rows[idx], cols[idx]] # 读取对应的NDVI影像（8月晴天的单景或合成） ndvi_ds = rasterio.open('./data/yunnan_aug_ndvi.tif') ndvi_vals = np.array([x[0] for x in sample_gen(ndvi_ds, samples)]) # 按类别统计NDVI中位数 for cls in np.unique(class_vals): mask = class_vals == cls if mask.sum() > 10: print('类', cls, '中位数NDVI', np.median(ndvi_vals[mask]))

这段代码用了rasterio的sample_gen，直接读取采样点处的NDVI值。注意两个栅格必须保持地理坐标一致，如果分类图已经投影到UTM 47N，NDVI影像也要用同样投影，否则采样点会落到错位位置。这也是为什么我在第3章坚持先把投影统一：后面所有验证步骤都依赖坐标系一致。

做完这步，你会对这份数据的“坑位”有直观认识。比如我跑到云南南部，发现常绿阔叶林的NDVI中位数只有0.55，而北部针叶林能到0.75，这其实是物候差异，不一定分类错；但如果某块“水体”的NDVI中位数是0.45，那绝对是分错了。

我个人的习惯是：无论数据来自哪里，永远保留一份未改动的原始tif，然后把所有中间产物都命名为带处理后缀的版本。这样做是吃过亏之后的教训——有次我直接改原图，结果后面想重新对比不同滤波参数时，发现原始版本已经救不回来了，只能重新解压。数据备份是最后的后悔药，而这套NDVI自检流程，则能在你向项目汇报前，提前把明显错分类的区域识别出来。希望帮到你。

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