尧图网站建设 
简介资源为ESRI发布的2021年云南省10米分辨率土地利用栅格数据集面向地理信息研究者、城乡规划人员及生态环保从业者可用于土地覆盖分类、环境变化监测、区域对比分析与GIS基础数据准备。压缩包共112个文件以tif栅格影像为核心其中xlsx与dbf存放属性记录tfw提供坐标配准信息xml保存元数据cpg标注编码格式png可作为快速预览图整体包体115.96MB。数据集按云南省各州市拆分组织西双版纳、文山、普洱等地均有独立文件便于按行政区局部调用同时基于WGS84坐标系并依大陆架范围裁剪可直接与全球其他地区数据叠加。已有247人学习。高分辨率栅格能够清晰反映10米尺度的土地利用差异为土地资源调查、环境影响评估、城市规划等应用提供可靠的基础数据。1. 一张 10 米土地利用图为什么值得单独拆开讲做云南项目的同行普遍有个抱怨10 m 分辨率的土地利用数据覆盖全省动不动就是几十 GB下载慢、处理卡、出图乱。但 ESRI 在 2021 年发布的这版云南省 10 m 精度土地利用产品反而值得你专门腾出时间拆一下。它比 30 m 的 GlobeLand30 在梯田、村庄、林缘灌丛这类小图斑上清楚得多而且类别体系统一省掉你自己标注和训练样本的环节。适合做国土空间监测、生态评估、农业估产和论文底图的人。问题是很多人只知道它精度高却不知道它的坐标系、类别编码和下载方式里藏着多少坑。2. ESRI 2021 云南 10m 土地利用先把数据源、类别和投影读懂2.1 数据源与产品定位不是影像是别人替你分好的结果我最初拿到这份云南省数据时第一反应是去翻原始影像结果浪费了两天。后来才明白ESRI 2021 年的土地利用产品并不是把 Sentinel-2 原始影像打包给你而是用深度学习模型对 2021 年度 Sentinel-2 L2A 时序影像做完推理之后输出的分类结果。一个像元的值代表地类编码比如 1 是水体、2 是树木。它解决的是“没有样本怎么做全云南地类底图”的问题帮你省掉自己标注和训练环节。为什么是 Sentinel-2 而不是 LandsatLandsat 只有 30 m 分辨率云南这种山高谷深、地块破碎的地方30 m 很容易把一条乡村公路和两边的灌丛混成一个像元梯田的层级关系也不明显。Sentinel-2 的 10 m 波段至少能把小地块边缘勾出来对识别村庄、林窗、小农田有很大帮助。ESRI 官方给的是 2021 年度合成分类结果不是单一天影像所以比抽查单景更能代表全年地表状态。数据由 Impact Observatory 生产、ESRI 分发在 ArcGIS Living Atlas 上以影像服务形式发布。你在 ArcGIS Pro 里搜 “ESRI Land Cover 2021” 能直接找到对个人研究和教学来说直接用它自带的分类结果比从头跑模型要快得多。但要注意它还没细到能代替国土调查那种人工审计图斑本质是模型推理结果2021 年也只代表该年度多数时相的“主要地类”。这句话决定后面所有操作方式你不能像处理遥感影像一样对它做平滑、插值和均值滤波因为分类结果是名义变量。像元值 1 和 2 加起来没有意义1.5 也是错误值。凡是改变像元值的重采样都要用最近邻法nearest凡是统计面积必须先明确 NoData 和投影接下来两节就讲这两个关键点。2.2 类别编码与字段先看懂 0~9 再动手ESRI 这个产品把全球土地覆盖归成 10 类官方字段名一般是“Class”像素类型多为 Int8。云南大部分地方能见到 1、2、3、5、6、7、8 这几种4 在湿地和坝塘附近9 在高海拔常年积雪和冰川区域。下表是我一直在用的编码对照也是后面做样式表、面积统计时的基准。编码英文类名中文示意云南常见分布0No Data空值/无数据边界外、云区、无效观测1Water水体江河、湖泊、水库、坝塘2Trees乔木/森林滇西北、滇南、哀牢山等林区3Grass草地高山草甸、林间草地4Flooded Vegetation洪水淹没植被/湿地湖滨、湿地、河边5Crops耕地坝区农田、梯田、台地6Scrub/Shrub灌丛干热河谷、退化林地7Built Area建设用地城镇、村庄、道路周边8Bare Ground裸地高山流石滩、裸露坡面、矿山9Snow/Ice冰雪/冰川高海拔常年积雪区这套类别体系的问题是太粗。你没法分出常绿阔叶林和针叶林也没法分出城市和农村居民点。做省级分析可以直接用做到县级或项目级就要再用辅助数据重编码。比如我做云南的农业项目时会把 5 保留把 3 和 6 在坡度小于 8° 且靠近 5 的区域重新判作休耕或撂荒农田这个规则只能结合地形自己做。还有一点容易踩坑0 在分类结果里是“无数据”类不是传统影像里的背景值。它可能代表原始影像被云覆盖也可能代表剪切后位于边界外。任何统计命令里0 都不能被当成一个真正的土地覆盖类别去算面积但它在显示时又要保留否则边界和云洞会变成空洞。2.3 坐标系与瓦片结构服务是 Web Mercator统计面积别在里面做Living Atlas 上这套影像服务发布坐标系是 WGS 1984 Web Mercator也就是 EPSG:3857。瓦片金字塔是正方形瓦片但 Web Mercator 有个特性纬度越高地面面积被放大得越厉害。在云南北纬 21°~29° 的区域如果直接在 3857 下计算像元面积面积会比真实值偏大 15%~22%这已经不是可以忽略的误差。所以下载之后的第一件事是把它重投影到等积坐标系。全省尺度的我一般用 WGS 84 / Albers Conic Equal Area标准纬线取 25°N 和 47°N中央经线取 105°E。这个选择和云南省地理跨度是配合的两条标准纬线把变形均匀分配到全省105°E 大致压住云南的中部边界东西两端的面积误差在可接受范围内。如果你是做某个州市或流域也可以改用 UTM 47N 或 48N但全省统计我就认 Albers。瓦片结构还限制了你下载的方式。影像服务默认按 256×256 瓦片组织请求整省范围容易撞上服务端最大像元数限制。常见的做法是先按行政边界切出范围再让 GDAL 分块读取和重投影而不是一次把全图塞进内存。这也是下一章那条命令为什么把重投影、裁剪和压缩放在一起的原因。3. 下载与预处理把服务变成云南本地的 Albers 等积 TIFF3.1 先确认你拿到的是服务地址还是本地栅格如果你是第一次接触这份数据先别急着跑重投影。我一般按两条路线走如果只有影像服务地址就在 ArcGIS Pro 里加进来或者直接用 QGIS 的 ArcGIS Image Server 连接把它当图层挂载如果已经拿到按图幅切好的 GeoTIFF就直接用 GDAL 做本地处理。从 Living Atlas 服务导出时不要简单点 Export Raster 就完事。服务端会默认沿用 Web Mercator 和全球范围整幅导出的结果可能大到难以处理而且分辨率会被金字塔抽稀。我更建议先用云南省边界确定导出范围在导出界面里设置 Processing Extent 为 yunnan_boundaryCell Size 填 10Resampling Method 选 Nearest NeighborOutput NoData 为 0。这个步骤只是把数据从服务端搬下来不承担投影转换。如果你已经有本地栅格第一步改成跑gdalinfo重点看四样东西投影、像元大小、像素类型、NoData 值。这三个信息决定后面所有参数。命令行贴在这里方便你快速试gdalinfo esri_yn_2021_raw.tif | grep -E Size|Pixel Size|Data type|NoData|Coordinate这段命令会输出栅格宽高、像元分辨率、数据类型、NoData 和坐标系。看到EPSG:3857不需要慌它只是说明你还没到统计阶段。看到像素类型是Byte或UInt8说明类别码可以直接用整数索引不需要转 float。3.2 一条 GDAL 命令完成重投影、裁剪和压缩确认原始文件没问题后我通常用下面这条gdalwarp把三件事一次做完。它把 Web Mercator 重投影到 Albers 等积投影同时按云南省边界裁剪最后输出压缩过的 GeoTIFF。命令如下gdalwarp -overwrite \ -s_srs EPSG:3857 \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 10 10 -r near \ -cutline yunnan_boundary.shp -crop_to_cutline \ -dstnodata 0 \ -co COMPRESSDEFLATE -co ZLEVEL9 -co TILEDYES \ esri_yn_2021_raw.tif esri_yn_2021_albers.tif-t_srs后面的 proj 字符串是 Albers 等积投影的核心参数。lat_125和lat_247是两条标准纬线lon_0105是中央经线x_00和y_00是假东和假北。这套参数覆盖云南全省时面积变形控制在较小范围内。如果你要严格对齐某套已有坐标可以直接把这里的datumWGS84换成datumCGCS2000或者直接填 EPSG 代码但多数项目用 WGS84 够用。-tr 10 10强制像元大小为 10m×10m。-r near指定最近邻重采样。类别栅格是离散变量不能用 bilinear 或 cubic否则会在森林和草地之间造出“2.5 类”这是很多人拿到结果后发现图斑边缘发花的直接原因。-cutline指定云南省界-crop_to_cutline让输出只保留边界内像元。-dstnodata 0把边界外和无效区域写成 0和官方类别 0 保持一致。-co COMPRESSDEFLATE -co ZLEVEL9 -co TILEDYES是输出参数。DEFLATE 压缩比高ZLEVEL 设为 9 会让文件更小但压缩稍慢TILEDYES 能让后续裁剪、统计和显示更快。如果你的原始文件范围特别大还可以加-wm 512限制缓存避免内存占满。3.3 用 Python 校核类别和面积重投影裁剪完之后我会立刻跑一段 Python 做检查。目的是确认类别码没有出现异常值也确认 0 值没有把统计口径污染掉。代码用 rasterio 读栅格对有效像元做计数并换算成公顷import numpy as np import rasterio with rasterio.open(esri_yn_2021_albers.tif) as src: data src.read(1) valid data 0 counts np.bincount(data[valid], minlength10) cell_m2 src.res[0] * src.res[1] for c in range(1, 10): if counts[c] 0: print(fClass {c}: {counts[c]} px, {counts[c] * cell_m2 / 1e4:.1f} ha)data 0是关键条件。它把类别 0 排除在统计外避免把边界外像元也算成地类面积。np.bincount对 0~9 的整数做计数minlength10 保证缺少某个类别时也能输出 0。cell_m2由 Albers 投影下的 10m 像元直接算出 100²公顷就除以 10000。如果输出里某个类别的像素数比预想高很多或者出现超过 9 的像素值回头查 gdalwarp 的-r near参数有没有写错。要是看到 255 或负值说明输入文件被浮点转化过需要重新从原始服务导出 Int8 数据。4. QGIS/ArcGIS 出图与面积统计别让图例和面积对不上4.1 唯一值渲染动手前先把颜色固定QGIS 打开 GeoTIFF 默认是按灰度拉伸显示的这时候裸地的 8 可能被显示成刺眼的亮白色建设用地的 7 反而偏暗如果你直接拿去出图就会翻车。正确做法是把图层渲染方式改成唯一值按类别码 0~9 分别给颜色。这张色表是我常用的版本主要照顾色盲友好和打印对比度编码颜色RGB说明0黑色0,0,0NoData1蓝色65,157,211水体2深绿58,110,74森林3浅绿136,170,119草地4青绿154,208,163湿地5黄褐色228,192,112耕地6土黄165,143,111灌丛7红色196,75,75建设用地8灰白214,204,192裸地9白色230,237,243冰雪在 QGIS 里双击图层Symbology 选择 Categorized字段选 Class然后新建 10 个类别逐个输入上面的色值。在 ArcGIS Pro 里对应的是 Unique Values 渲染同样按 VALUE 字段分 10 类。这个动作会直接决定你后面导出的图例名称是否正确千万不要用默认 ESRI 分类色带因为 ESRI 官方的符号库在不同版本里对同一编码的显示颜色并不一致。色表固定的另一个好处是当你把 2020 和 2021 两年的图放在同一版面时同一地类颜色稳定读者不会因为图例不一致而产生误读。尤其在做时间对比专题图时颜色语义一致比“好看”重要得多。4.2 面积统计Albers 投影下按类别汇总如果你只是想快速得到全省各类面积建议用上一章的 Python 方法。如果你需要输出带图斑边界的统计表可以先转成矢量再汇总也可以直接在栅格上做分区统计。我这里给一个更适用于机构内部项目的做法用 rasterio 的 mask 配合县界 shp逐县统计各类面积。import geopandas as gpd import numpy as np import rasterio from rasterio.mask import mask county gpd.read_file(yunnan_counties.shp) with rasterio.open(esri_yn_2021_albers.tif) as src: for _, row in county.iterrows(): out, transform mask(src, [row.geometry], cropTrue, nodata0, all_touchedFalse) data out[0] valid data 0 counts np.bincount(data[valid], minlength10) cell_m2 transform[0] * transform[4] print(row[NAME], counts, counts * cell_m2 / 1e4)mask函数里的cropTrue会按当前县界范围裁剪栅格nodata0把裁剪后边界外填成 0all_touchedFalse表示只有中心点落在县界内的像元才被保留。如果某个县边界特别细长像元中心刚好不落在边界内的概率会变大此时可以把all_touchedTrue代价是边界附近像元会被重复计入相邻县统计总和会略偏大。做省级汇总时我一般全图统计做县级台账时则明确记录用了哪种参数。这里要再次强调以上统计必须发生在 Albers 等积投影上。如果你拿 Web Mercator 做同样操作面积偏大 15%~22%而且纬度越高偏得越多最后写报告时很难跟同事解释。4.3 分区统计的顺序先掩膜再计数别先裁剪再拼接有一个常见误操作是先按县界把全省裁剪成 129 个县发到群里让别人分别统计最后再把结果拼起来。这样很容易出现接边处类别码错位和面积重复计算。我一般会先保留一份全省完整 tif再用rasterio.mask逐个县拉取而不是把同一个原始数据切碎。如果项目必须分幅处理记住两件事每个分幅文件要保留至少 50m 重叠区最后镶嵌时用gdal_merge而不是直接栅格覆盖分幅处理后的 NoData 必须统一写成 0否则镶嵌后边界会出现一条一条的接缝。这条血泪经验我吃过亏后面会专门写进避坑清单。5. 避坑与排查云南山地数据最容易翻车的五个地方5.1 导出影像出现大面积 0 值空洞现象从影像服务导出后森林或农田区域中间出现一大片规则或不规则的黑色空洞统计面积时该类面积明显偏小。原因服务端是按瓦片返回数据的高分辨率请求如果超过服务限制会静默丢弃部分瓦片另外 Sentinel-2 在云南多云季节有大量云掩膜模型把云区直接标成 0这也会造成真实空洞。解决先看空洞是否沿瓦片边界分布。沿瓦片边界的是下载中断重新用 gdalwarp 并加上-co TILEDYES分段重试即可。云掩膜造成的 0 则需要保留因为它代表“该位置没有可用观测”如果你非要补洞只能拿相邻年份同月产品或辅助影像做变化填充不要用中值滤波否则会把真实地类边界抹掉。5.2 在 Web Mercator 里统计面积结果整体偏大现象我用同一份云南省边界统计森林面积Albers 投影下是 14.2 万平方公里Web Mercator 投影下算出来接近 17 万平方公里整整大了 20%。原因Web Mercator 在赤道附近变形小纬度越高面积放大越明显。云南大部分地区在北纬 21°~29°面积变形达到 15%~22%直接统计必然偏大。解决所有面积统计前先检查src.crs。只要看到 EPSG:3857就执行一次gdalwarp到等积投影。统计完在报告里注明“面积按 WGS 84 / Albers Conic Equal Area 计算”这样别人复算时不会怀疑你数据有问题。5.3 农田和灌丛混淆坝区耕地被分成碎块现象在元江、版纳、怒江河谷等干热河谷区域大片 5耕地和 6灌丛交错出现梯田边界像锯齿一样单独看某一块图斑根本分不清是田还是灌丛。原因ESRI 分类模型依赖 Sentinel-2 时序光谱雨季和旱季的植被状态变化会造成同物异谱云南山地坡度大阴影和裸岩也会让模型把作物误判为灌丛。解决别单独用 5 或 6 做精细农业分析。可以结合坡度图和 NDVI 时序做后处理把坡度小于 8°、且被 5 包围的 6 重新归类为农田把高海拔、坡度大于 25° 的 5 重新归类为灌丛或草地。这一步需要你自己有地形数据ESRI 官方分类结果解决不了。5.4 图例颜色和类别码错位现象同样的 7 编码在某个配色文件里是建设用地换到另一个 lyr 样式后变成了裸地打印出来的图例上 8 写的是水体但实际是裸地。原因唯一值渲染按整数匹配如果样式文件来自旧版产品或某个自定义 LUT类别 7、8 的颜色和名称会被带偏。尤其是从 GeoTIFF 直接拖进 QGIS 时QGIS 默认的色带完全不知道 1 是水、7 是城市。解决每次换机器、换项目前先加载上文那张 0~9 编码表确认一遍。我习惯把颜色表存成独立 CSV 或 qml 文件和 tif 放在同一个数据目录里不依赖 ArcGIS 的默认样式库。这样即使同事打开数据也不会出现“红色是裸地”这种低级错误。5.5 裁剪后黑边被当成有效面积现象按县界或省界裁剪后输出的 tif 边缘有一圈黑色多边形面积统计时总类面积偏大或者矢量转面后会沿边界生成一圈细长图斑。原因-crop_to_cutline本身会把边界外写成 NoData但如果你的-dstnodata没写好边界外的值可能是 0而 0 同时又是官方类别 0。统计时如果不做data 0过滤这些边界像元就被当成一个“类别”参与计算。解决统计前像 3.3 那样用valid data 0做掩膜。如果转矢量先执行一步gdal_translate -mask_mode或者直接删除值为 0 的图斑。还要检查行政边界本身是否覆盖了邻省飞地云南有几个县界在实地勘界后会动态调整直接用最新省界比用旧版稳定得多。6. 进阶把分类结果转成矢量地块并做两年变化检测如果你的分析不满足于“全云南有几个像素是森林”而是想拿到具体地块边界下一步就是把分类栅格矢量化成多边形。GDAL 自带gdal_polygonize.py对 10 m 栅格来说生成的面数量会非常大直接转出来的图斑通常碎到没法用。我一般先做一次 Majority Filter把小于 4~8 个像元的碎斑合并到邻域主类再转面最后用面积阈值删掉小于 0.01 平方公里的零碎图斑。# 先做众数滤波消除孤立碎斑 gdal_fillnodata.py -md 200 esri_yn_2021_albers.tif filled.tif # 或者用 OpenCV 的 mode filter但 GDAL 内置更方便 gdal_polygonize.py esri_yn_2021_albers.tif -mask 0 -8 \ -b 1 classes.gpkg -f GPKG-mask 0意思是把像素值 0 当 NoData不参与转面-8字段名是整理后的类名你也可以改成class。转面后通常要跑一次 Dissolve按DN字段合并相邻同类别多边形再在 QGIS 里按面积字段筛选。这个流程做完你手里就是一套带边界的土地利用矢量图层可以直接出专题图或做叠置分析。如果手头还有 2020 年的同类数据变化检测就变得很直观把两年栅格重投影到同一个网格上逐像元对比类别码变化区域取new_code ! old_code。要注意重投影网格的锚点必须对齐我一般用-te和-tr强制指定相同的范围和像元大小否则两年的边界会因为半个像元的偏移产生大量“伪变化”。算完变化后把变化区域做成掩膜按行政区统计各类转入转出面积这份结果比直接看遥感影像更接近土地利用转移矩阵。从那以后我每次做云南土地利用底图都会先花十分钟把投影、0 值、类别对应这三件事写进项目说明再开始出图省掉了后来无数返工。希望帮到你。本文还有配套的精品资源点击获取 云南省10米土地利用数据下载与预处理：从投影转换到面积统计全攻略
简介资源为ESRI发布的2021年云南省10米分辨率土地利用栅格数据集面向地理信息研究者、城乡规划人员及生态环保从业者可用于土地覆盖分类、环境变化监测、区域对比分析与GIS基础数据准备。压缩包共112个文件以tif栅格影像为核心其中xlsx与dbf存放属性记录tfw提供坐标配准信息xml保存元数据cpg标注编码格式png可作为快速预览图整体包体115.96MB。数据集按云南省各州市拆分组织西双版纳、文山、普洱等地均有独立文件便于按行政区局部调用同时基于WGS84坐标系并依大陆架范围裁剪可直接与全球其他地区数据叠加。已有247人学习。高分辨率栅格能够清晰反映10米尺度的土地利用差异为土地资源调查、环境影响评估、城市规划等应用提供可靠的基础数据。1. 一张 10 米土地利用图为什么值得单独拆开讲做云南项目的同行普遍有个抱怨10 m 分辨率的土地利用数据覆盖全省动不动就是几十 GB下载慢、处理卡、出图乱。但 ESRI 在 2021 年发布的这版云南省 10 m 精度土地利用产品反而值得你专门腾出时间拆一下。它比 30 m 的 GlobeLand30 在梯田、村庄、林缘灌丛这类小图斑上清楚得多而且类别体系统一省掉你自己标注和训练样本的环节。适合做国土空间监测、生态评估、农业估产和论文底图的人。问题是很多人只知道它精度高却不知道它的坐标系、类别编码和下载方式里藏着多少坑。2. ESRI 2021 云南 10m 土地利用先把数据源、类别和投影读懂2.1 数据源与产品定位不是影像是别人替你分好的结果我最初拿到这份云南省数据时第一反应是去翻原始影像结果浪费了两天。后来才明白ESRI 2021 年的土地利用产品并不是把 Sentinel-2 原始影像打包给你而是用深度学习模型对 2021 年度 Sentinel-2 L2A 时序影像做完推理之后输出的分类结果。一个像元的值代表地类编码比如 1 是水体、2 是树木。它解决的是“没有样本怎么做全云南地类底图”的问题帮你省掉自己标注和训练环节。为什么是 Sentinel-2 而不是 LandsatLandsat 只有 30 m 分辨率云南这种山高谷深、地块破碎的地方30 m 很容易把一条乡村公路和两边的灌丛混成一个像元梯田的层级关系也不明显。Sentinel-2 的 10 m 波段至少能把小地块边缘勾出来对识别村庄、林窗、小农田有很大帮助。ESRI 官方给的是 2021 年度合成分类结果不是单一天影像所以比抽查单景更能代表全年地表状态。数据由 Impact Observatory 生产、ESRI 分发在 ArcGIS Living Atlas 上以影像服务形式发布。你在 ArcGIS Pro 里搜 “ESRI Land Cover 2021” 能直接找到对个人研究和教学来说直接用它自带的分类结果比从头跑模型要快得多。但要注意它还没细到能代替国土调查那种人工审计图斑本质是模型推理结果2021 年也只代表该年度多数时相的“主要地类”。这句话决定后面所有操作方式你不能像处理遥感影像一样对它做平滑、插值和均值滤波因为分类结果是名义变量。像元值 1 和 2 加起来没有意义1.5 也是错误值。凡是改变像元值的重采样都要用最近邻法nearest凡是统计面积必须先明确 NoData 和投影接下来两节就讲这两个关键点。2.2 类别编码与字段先看懂 0~9 再动手ESRI 这个产品把全球土地覆盖归成 10 类官方字段名一般是“Class”像素类型多为 Int8。云南大部分地方能见到 1、2、3、5、6、7、8 这几种4 在湿地和坝塘附近9 在高海拔常年积雪和冰川区域。下表是我一直在用的编码对照也是后面做样式表、面积统计时的基准。编码英文类名中文示意云南常见分布0No Data空值/无数据边界外、云区、无效观测1Water水体江河、湖泊、水库、坝塘2Trees乔木/森林滇西北、滇南、哀牢山等林区3Grass草地高山草甸、林间草地4Flooded Vegetation洪水淹没植被/湿地湖滨、湿地、河边5Crops耕地坝区农田、梯田、台地6Scrub/Shrub灌丛干热河谷、退化林地7Built Area建设用地城镇、村庄、道路周边8Bare Ground裸地高山流石滩、裸露坡面、矿山9Snow/Ice冰雪/冰川高海拔常年积雪区这套类别体系的问题是太粗。你没法分出常绿阔叶林和针叶林也没法分出城市和农村居民点。做省级分析可以直接用做到县级或项目级就要再用辅助数据重编码。比如我做云南的农业项目时会把 5 保留把 3 和 6 在坡度小于 8° 且靠近 5 的区域重新判作休耕或撂荒农田这个规则只能结合地形自己做。还有一点容易踩坑0 在分类结果里是“无数据”类不是传统影像里的背景值。它可能代表原始影像被云覆盖也可能代表剪切后位于边界外。任何统计命令里0 都不能被当成一个真正的土地覆盖类别去算面积但它在显示时又要保留否则边界和云洞会变成空洞。2.3 坐标系与瓦片结构服务是 Web Mercator统计面积别在里面做Living Atlas 上这套影像服务发布坐标系是 WGS 1984 Web Mercator也就是 EPSG:3857。瓦片金字塔是正方形瓦片但 Web Mercator 有个特性纬度越高地面面积被放大得越厉害。在云南北纬 21°~29° 的区域如果直接在 3857 下计算像元面积面积会比真实值偏大 15%~22%这已经不是可以忽略的误差。所以下载之后的第一件事是把它重投影到等积坐标系。全省尺度的我一般用 WGS 84 / Albers Conic Equal Area标准纬线取 25°N 和 47°N中央经线取 105°E。这个选择和云南省地理跨度是配合的两条标准纬线把变形均匀分配到全省105°E 大致压住云南的中部边界东西两端的面积误差在可接受范围内。如果你是做某个州市或流域也可以改用 UTM 47N 或 48N但全省统计我就认 Albers。瓦片结构还限制了你下载的方式。影像服务默认按 256×256 瓦片组织请求整省范围容易撞上服务端最大像元数限制。常见的做法是先按行政边界切出范围再让 GDAL 分块读取和重投影而不是一次把全图塞进内存。这也是下一章那条命令为什么把重投影、裁剪和压缩放在一起的原因。3. 下载与预处理把服务变成云南本地的 Albers 等积 TIFF3.1 先确认你拿到的是服务地址还是本地栅格如果你是第一次接触这份数据先别急着跑重投影。我一般按两条路线走如果只有影像服务地址就在 ArcGIS Pro 里加进来或者直接用 QGIS 的 ArcGIS Image Server 连接把它当图层挂载如果已经拿到按图幅切好的 GeoTIFF就直接用 GDAL 做本地处理。从 Living Atlas 服务导出时不要简单点 Export Raster 就完事。服务端会默认沿用 Web Mercator 和全球范围整幅导出的结果可能大到难以处理而且分辨率会被金字塔抽稀。我更建议先用云南省边界确定导出范围在导出界面里设置 Processing Extent 为 yunnan_boundaryCell Size 填 10Resampling Method 选 Nearest NeighborOutput NoData 为 0。这个步骤只是把数据从服务端搬下来不承担投影转换。如果你已经有本地栅格第一步改成跑gdalinfo重点看四样东西投影、像元大小、像素类型、NoData 值。这三个信息决定后面所有参数。命令行贴在这里方便你快速试gdalinfo esri_yn_2021_raw.tif | grep -E Size|Pixel Size|Data type|NoData|Coordinate这段命令会输出栅格宽高、像元分辨率、数据类型、NoData 和坐标系。看到EPSG:3857不需要慌它只是说明你还没到统计阶段。看到像素类型是Byte或UInt8说明类别码可以直接用整数索引不需要转 float。3.2 一条 GDAL 命令完成重投影、裁剪和压缩确认原始文件没问题后我通常用下面这条gdalwarp把三件事一次做完。它把 Web Mercator 重投影到 Albers 等积投影同时按云南省边界裁剪最后输出压缩过的 GeoTIFF。命令如下gdalwarp -overwrite \ -s_srs EPSG:3857 \ -t_srs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs \ -tr 10 10 -r near \ -cutline yunnan_boundary.shp -crop_to_cutline \ -dstnodata 0 \ -co COMPRESSDEFLATE -co ZLEVEL9 -co TILEDYES \ esri_yn_2021_raw.tif esri_yn_2021_albers.tif-t_srs后面的 proj 字符串是 Albers 等积投影的核心参数。lat_125和lat_247是两条标准纬线lon_0105是中央经线x_00和y_00是假东和假北。这套参数覆盖云南全省时面积变形控制在较小范围内。如果你要严格对齐某套已有坐标可以直接把这里的datumWGS84换成datumCGCS2000或者直接填 EPSG 代码但多数项目用 WGS84 够用。-tr 10 10强制像元大小为 10m×10m。-r near指定最近邻重采样。类别栅格是离散变量不能用 bilinear 或 cubic否则会在森林和草地之间造出“2.5 类”这是很多人拿到结果后发现图斑边缘发花的直接原因。-cutline指定云南省界-crop_to_cutline让输出只保留边界内像元。-dstnodata 0把边界外和无效区域写成 0和官方类别 0 保持一致。-co COMPRESSDEFLATE -co ZLEVEL9 -co TILEDYES是输出参数。DEFLATE 压缩比高ZLEVEL 设为 9 会让文件更小但压缩稍慢TILEDYES 能让后续裁剪、统计和显示更快。如果你的原始文件范围特别大还可以加-wm 512限制缓存避免内存占满。3.3 用 Python 校核类别和面积重投影裁剪完之后我会立刻跑一段 Python 做检查。目的是确认类别码没有出现异常值也确认 0 值没有把统计口径污染掉。代码用 rasterio 读栅格对有效像元做计数并换算成公顷import numpy as np import rasterio with rasterio.open(esri_yn_2021_albers.tif) as src: data src.read(1) valid data 0 counts np.bincount(data[valid], minlength10) cell_m2 src.res[0] * src.res[1] for c in range(1, 10): if counts[c] 0: print(fClass {c}: {counts[c]} px, {counts[c] * cell_m2 / 1e4:.1f} ha)data 0是关键条件。它把类别 0 排除在统计外避免把边界外像元也算成地类面积。np.bincount对 0~9 的整数做计数minlength10 保证缺少某个类别时也能输出 0。cell_m2由 Albers 投影下的 10m 像元直接算出 100²公顷就除以 10000。如果输出里某个类别的像素数比预想高很多或者出现超过 9 的像素值回头查 gdalwarp 的-r near参数有没有写错。要是看到 255 或负值说明输入文件被浮点转化过需要重新从原始服务导出 Int8 数据。4. QGIS/ArcGIS 出图与面积统计别让图例和面积对不上4.1 唯一值渲染动手前先把颜色固定QGIS 打开 GeoTIFF 默认是按灰度拉伸显示的这时候裸地的 8 可能被显示成刺眼的亮白色建设用地的 7 反而偏暗如果你直接拿去出图就会翻车。正确做法是把图层渲染方式改成唯一值按类别码 0~9 分别给颜色。这张色表是我常用的版本主要照顾色盲友好和打印对比度编码颜色RGB说明0黑色0,0,0NoData1蓝色65,157,211水体2深绿58,110,74森林3浅绿136,170,119草地4青绿154,208,163湿地5黄褐色228,192,112耕地6土黄165,143,111灌丛7红色196,75,75建设用地8灰白214,204,192裸地9白色230,237,243冰雪在 QGIS 里双击图层Symbology 选择 Categorized字段选 Class然后新建 10 个类别逐个输入上面的色值。在 ArcGIS Pro 里对应的是 Unique Values 渲染同样按 VALUE 字段分 10 类。这个动作会直接决定你后面导出的图例名称是否正确千万不要用默认 ESRI 分类色带因为 ESRI 官方的符号库在不同版本里对同一编码的显示颜色并不一致。色表固定的另一个好处是当你把 2020 和 2021 两年的图放在同一版面时同一地类颜色稳定读者不会因为图例不一致而产生误读。尤其在做时间对比专题图时颜色语义一致比“好看”重要得多。4.2 面积统计Albers 投影下按类别汇总如果你只是想快速得到全省各类面积建议用上一章的 Python 方法。如果你需要输出带图斑边界的统计表可以先转成矢量再汇总也可以直接在栅格上做分区统计。我这里给一个更适用于机构内部项目的做法用 rasterio 的 mask 配合县界 shp逐县统计各类面积。import geopandas as gpd import numpy as np import rasterio from rasterio.mask import mask county gpd.read_file(yunnan_counties.shp) with rasterio.open(esri_yn_2021_albers.tif) as src: for _, row in county.iterrows(): out, transform mask(src, [row.geometry], cropTrue, nodata0, all_touchedFalse) data out[0] valid data 0 counts np.bincount(data[valid], minlength10) cell_m2 transform[0] * transform[4] print(row[NAME], counts, counts * cell_m2 / 1e4)mask函数里的cropTrue会按当前县界范围裁剪栅格nodata0把裁剪后边界外填成 0all_touchedFalse表示只有中心点落在县界内的像元才被保留。如果某个县边界特别细长像元中心刚好不落在边界内的概率会变大此时可以把all_touchedTrue代价是边界附近像元会被重复计入相邻县统计总和会略偏大。做省级汇总时我一般全图统计做县级台账时则明确记录用了哪种参数。这里要再次强调以上统计必须发生在 Albers 等积投影上。如果你拿 Web Mercator 做同样操作面积偏大 15%~22%而且纬度越高偏得越多最后写报告时很难跟同事解释。4.3 分区统计的顺序先掩膜再计数别先裁剪再拼接有一个常见误操作是先按县界把全省裁剪成 129 个县发到群里让别人分别统计最后再把结果拼起来。这样很容易出现接边处类别码错位和面积重复计算。我一般会先保留一份全省完整 tif再用rasterio.mask逐个县拉取而不是把同一个原始数据切碎。如果项目必须分幅处理记住两件事每个分幅文件要保留至少 50m 重叠区最后镶嵌时用gdal_merge而不是直接栅格覆盖分幅处理后的 NoData 必须统一写成 0否则镶嵌后边界会出现一条一条的接缝。这条血泪经验我吃过亏后面会专门写进避坑清单。5. 避坑与排查云南山地数据最容易翻车的五个地方5.1 导出影像出现大面积 0 值空洞现象从影像服务导出后森林或农田区域中间出现一大片规则或不规则的黑色空洞统计面积时该类面积明显偏小。原因服务端是按瓦片返回数据的高分辨率请求如果超过服务限制会静默丢弃部分瓦片另外 Sentinel-2 在云南多云季节有大量云掩膜模型把云区直接标成 0这也会造成真实空洞。解决先看空洞是否沿瓦片边界分布。沿瓦片边界的是下载中断重新用 gdalwarp 并加上-co TILEDYES分段重试即可。云掩膜造成的 0 则需要保留因为它代表“该位置没有可用观测”如果你非要补洞只能拿相邻年份同月产品或辅助影像做变化填充不要用中值滤波否则会把真实地类边界抹掉。5.2 在 Web Mercator 里统计面积结果整体偏大现象我用同一份云南省边界统计森林面积Albers 投影下是 14.2 万平方公里Web Mercator 投影下算出来接近 17 万平方公里整整大了 20%。原因Web Mercator 在赤道附近变形小纬度越高面积放大越明显。云南大部分地区在北纬 21°~29°面积变形达到 15%~22%直接统计必然偏大。解决所有面积统计前先检查src.crs。只要看到 EPSG:3857就执行一次gdalwarp到等积投影。统计完在报告里注明“面积按 WGS 84 / Albers Conic Equal Area 计算”这样别人复算时不会怀疑你数据有问题。5.3 农田和灌丛混淆坝区耕地被分成碎块现象在元江、版纳、怒江河谷等干热河谷区域大片 5耕地和 6灌丛交错出现梯田边界像锯齿一样单独看某一块图斑根本分不清是田还是灌丛。原因ESRI 分类模型依赖 Sentinel-2 时序光谱雨季和旱季的植被状态变化会造成同物异谱云南山地坡度大阴影和裸岩也会让模型把作物误判为灌丛。解决别单独用 5 或 6 做精细农业分析。可以结合坡度图和 NDVI 时序做后处理把坡度小于 8°、且被 5 包围的 6 重新归类为农田把高海拔、坡度大于 25° 的 5 重新归类为灌丛或草地。这一步需要你自己有地形数据ESRI 官方分类结果解决不了。5.4 图例颜色和类别码错位现象同样的 7 编码在某个配色文件里是建设用地换到另一个 lyr 样式后变成了裸地打印出来的图例上 8 写的是水体但实际是裸地。原因唯一值渲染按整数匹配如果样式文件来自旧版产品或某个自定义 LUT类别 7、8 的颜色和名称会被带偏。尤其是从 GeoTIFF 直接拖进 QGIS 时QGIS 默认的色带完全不知道 1 是水、7 是城市。解决每次换机器、换项目前先加载上文那张 0~9 编码表确认一遍。我习惯把颜色表存成独立 CSV 或 qml 文件和 tif 放在同一个数据目录里不依赖 ArcGIS 的默认样式库。这样即使同事打开数据也不会出现“红色是裸地”这种低级错误。5.5 裁剪后黑边被当成有效面积现象按县界或省界裁剪后输出的 tif 边缘有一圈黑色多边形面积统计时总类面积偏大或者矢量转面后会沿边界生成一圈细长图斑。原因-crop_to_cutline本身会把边界外写成 NoData但如果你的-dstnodata没写好边界外的值可能是 0而 0 同时又是官方类别 0。统计时如果不做data 0过滤这些边界像元就被当成一个“类别”参与计算。解决统计前像 3.3 那样用valid data 0做掩膜。如果转矢量先执行一步gdal_translate -mask_mode或者直接删除值为 0 的图斑。还要检查行政边界本身是否覆盖了邻省飞地云南有几个县界在实地勘界后会动态调整直接用最新省界比用旧版稳定得多。6. 进阶把分类结果转成矢量地块并做两年变化检测如果你的分析不满足于“全云南有几个像素是森林”而是想拿到具体地块边界下一步就是把分类栅格矢量化成多边形。GDAL 自带gdal_polygonize.py对 10 m 栅格来说生成的面数量会非常大直接转出来的图斑通常碎到没法用。我一般先做一次 Majority Filter把小于 4~8 个像元的碎斑合并到邻域主类再转面最后用面积阈值删掉小于 0.01 平方公里的零碎图斑。# 先做众数滤波消除孤立碎斑 gdal_fillnodata.py -md 200 esri_yn_2021_albers.tif filled.tif # 或者用 OpenCV 的 mode filter但 GDAL 内置更方便 gdal_polygonize.py esri_yn_2021_albers.tif -mask 0 -8 \ -b 1 classes.gpkg -f GPKG-mask 0意思是把像素值 0 当 NoData不参与转面-8字段名是整理后的类名你也可以改成class。转面后通常要跑一次 Dissolve按DN字段合并相邻同类别多边形再在 QGIS 里按面积字段筛选。这个流程做完你手里就是一套带边界的土地利用矢量图层可以直接出专题图或做叠置分析。如果手头还有 2020 年的同类数据变化检测就变得很直观把两年栅格重投影到同一个网格上逐像元对比类别码变化区域取new_code ! old_code。要注意重投影网格的锚点必须对齐我一般用-te和-tr强制指定相同的范围和像元大小否则两年的边界会因为半个像元的偏移产生大量“伪变化”。算完变化后把变化区域做成掩膜按行政区统计各类转入转出面积这份结果比直接看遥感影像更接近土地利用转移矩阵。从那以后我每次做云南土地利用底图都会先花十分钟把投影、0 值、类别对应这三件事写进项目说明再开始出图省掉了后来无数返工。希望帮到你。本文还有配套的精品资源点击获取