ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

ESRI 2021年吉林省10米土地利用数据处理与面积统计实战

ESRI 2021年吉林省10米土地利用数据处理与面积统计实战 简介2021年土地利用ESRI吉林省10m精度数据集面向GIS从业者、城乡规划、农业及生态研究人员提供吉林省以地市为单位的高精度土地覆盖栅格成果。压缩包共63个文件、约44MB含9个地市TIF主影像及配套tfw地理配准文件、aux.xml元数据、dbf/cpg属性表、xlsx地类统计和png快速预览图便于直接叠加分析和出图。数据采用WGS84坐标系统10米分辨率可支撑土地覆盖变化监测、生态环境评估、城乡规划、农业管理等场景已有258人学习下载。借助属性表还能快速提取耕地、林地、水域、建设用地等类型面积适合作为省级尺度土地利用研究的入门底图与教学案例。1. ESRI 2021年吉林省10米土地利用数据一个不需要自己训练就能落地的现成底图2021年土地利用ESRI吉林省10m精度指的是Esri公司基于欧空局Sentinel-2影像制作的全球10米土地覆盖数据在2021年版本里提供了吉林省范围的栅格成果。我第一次拿它做黑土区耕地结构快速识别时只花了半天就把全省耕地分布和建设用地边界拉出了底数比自己跑标注快了一个量级。它不是高精度的地块图斑更替代不了外业调查但作为免费、现成、覆盖全年的10米级分类底图在省级前期研判、农业结构监测、生态评价和自然资源审计里非常能打。适合谁手上没有高质量标注样本、又想在短时间里拿到一个可解释地表覆盖格局的从业者。2. 看懂ESRI 10米分类底图层、类别编码与吉林省三大地貌区的真实表现2.1 Sentinel-2底图ESRI的10米分辨率是怎么来的ESRI这套全球土地覆盖产品的原始输入是欧空局Sentinel-2的L2A级地表反射率数据。Sentinel-2的波段里B2蓝、B3绿、B4红、B8近红外是原生10米分辨率B5、B6、B7、B8A这些红边波段和B11、B12短波红外波段是20米重采样。ESRI在制作2021年版时把全年可用的影像按季节合成、去云筛选再用深度学习模型逐像元分类。这里最关键的是“全年合成”四个字。因为使用的是全年多期影像冬季覆盖好的地区会引入大量雪盖信息夏季云雨多的地区则可能选不到完全无云的一帧。实际用下来吉林东部林区经常混入冰雪类别中部平原的耕地又可能在春播窗口被识别成裸地。这就是10米产品看上去漂亮但要付出的代价类别一共就那么十几个地物边界被模型平滑了季节性的错分却依然明显。还要注意一点ESRI的模型不是用传统决策树做的而是深度卷积网络训练标签本身是黑匣子官方只公布了分类结果和总体精度没公开每类的样本量。所以拿到tif以后别把分类码当成标准答案它更像一组“带先验概率的地表猜想”。2.2 类别编码表拿到tif先看唯一值再谈统计ESRI这套数据的类别编码在全球是统一的拿到吉林省切片后第一件事不是直接做面积统计而是把唯一值表拉出来看一眼。常见类别如下编码类别名称吉林省地面覆盖对应的典型对象1Water松花湖、查干湖、大型水库、松花江干流2Trees长白山阔叶红松林、落叶松人工林、村屯片林3Grass西部草原、路旁荒草、休耕草地4Flooded vegetation稻田泡田期、芦苇湿地、河漫滩植被5Crops玉米、大豆、水稻非灌水期6Built Area城镇建成区、农村宅基地、工矿厂房7Bare Ground春播前裸土、沙地、取土场、采矿迹地8Snow/Ice冬季残雪、冰面9Clouds合成时段残留云及云影需要特别说明的是不同年份版本的类别有微调。有些2021年再加工切片里会出现标号为10的类别对应灌丛或稀树草原某些渠道下载的省域tif为了压缩体积把无数据区域写成989或999还有分块之间NoData值不统一的情况。这些非常规值必须在统计前全部处理掉否则汇总面积时会出现一个比吉林省总面积还大的异常值。和武汉大学CLCD 30米数据、GlobeLand30对比来看ESRI这版最大的优势是10米分辨率对农村宅基地、田间道路、小片林带的刻画更接近真实地表劣势是类别数少林地细分、草地退化分级基本没有。如果你只做土地覆盖大类ESRI是免费方案里性价比最高的选择。2.3 吉林省三大地貌区与这套数据的适配度吉林省大致可以分成三块来看。东部延边、通化、白山一带是长白山林区天然林连片2类树木在连续林区里识别稳定但会被道路和小块农田打断形成大量长条状碎图斑。中部长春、四平、吉林市周边是典型黑土农业区田块规整5类农作物在大田生长期很稳定可每年四月底到五月初的播种窗口地表完全裸露7类裸地的比重会异常上升。西部白城、松原一带是农牧交错区分布着向海、莫莫格等湿地3类草地和4类洪水植被随季节水位变化互相串扰。这套数据在中部大农田区的表现最让人放心在西部湿地和东部过渡地带则需要额外验证。理解这三块的地貌差异之后后面下载裁剪时才知道该在哪个局部多做检查。ESRI数据是“快速底图”不是“地面真值”这一点越早想明白后期踩坑越少。3. 下载与裁剪吉林省范围从Living Atlas到本地tif的完整路径3.1 获取2021年吉林省切片的两种途径最常规的做法是打开Esri官方的Land Cover Explorer应用选年份2021在地图上框住吉林省范围等待后台生成tif下载。这个应用输出的是Web墨卡托分块的栅格一个省往往拆成好几块到手之后要拼接。不想用浏览器的话也可以在ArcGIS Pro的Catalog里从Living Atlas搜索“Sentinel-2 10-Meter Land Use/Land Cover”把对应图层拖进内容列表右键导出吉林范围。这里有个血泪经验在地图上手动拖矩形框选吉林省时由于吉林省东西跨度大下载结果经常被切成带状的若干块坐标系锁死为WGS84 Web墨卡托。块与块拼接时容易出现半个像元的边缘错位后续和国土调查数据叠加就会“差一点点”。我现在的习惯是直接用吉林省面状矢量作为范围参数提交下载而不是手动拖框省得后期拼接返工。3.2 换投影再裁剪吉林省为什么不能用Web墨卡托做统计ESRI默认输出的是WGS84 Web墨卡托投影。这个投影在吉林省所在的北纬41度到46度区间会把面积放大很多当屏幕底图没问题拿来做面积统计就会翻车。正确顺序是先转等积投影再做掩膜裁剪顺序反过来也能用但转投影时栅格重采样会把边界弄毛。我的标准流程是先把拼接好的原始栅格用“投影栅格”工具转到Albers等积投影吉林省范围我习惯设中央经线129°E双标准纬线43°N和46°N。分类栅格重采样必须用NEAREST最邻近法不要用BILINEAR或CUBIC否则原本编码是2的树木像元会被平均成小数后续统计直接报错。import arcpy src_raster rD:\lc\jl_esri2021_raw.tif albers_raster rD:\lc\jl_esri2021_albers.tif sr arcpy.SpatialReference() sr.loadFromString( PROJCS[China_Albers_129E, GEOGCS[GCS_China_Geodetic_Coordinate_System_2000, DATUM[D_China_2000,SPHEROID[CGCS2000,6378137.0,298.257222101]], PRIMEM[Greenwich,0],UNIT[Degree,0.0174532925199433]], PROJECTION[Albers_Conic_Equal_Area], PARAMETER[False_Easting,0], PARAMETER[False_Northing,0], PARAMETER[Central_Meridian,129.0], PARAMETER[Standard_Parallel_1,43.0], PARAMETER[Standard_Parallel_2,46.0], PARAMETER[Latitude_Of_Origin,40.0], UNIT[Meter,1.0]]) arcpy.ProjectRaster_management( src_raster, albers_raster, sr, resampling_typeNEAREST, cell_size10 10 ) clip_raster rD:\lc\jl_esri2021_clip.tif mask rD:\lc\jl_boundary.shp arcpy.Clip_management( albers_raster, #, clip_raster, mask, nodata_value0, clipping_geometryNODATA, maintain_clipping_extent# ) print(done)这段代码的逻辑是两步走第一步做投影转换把Web墨卡托全球分块变成适合吉林范围的等积投影同时把像元重采样为10米乘10米第二步用省界矢量做掩膜裁剪把省界外的像元去掉NoData统一写成0。需要留意的参数有三个一是resampling_type必须NEAREST二是cell_size要显式给“10 10”否则ArcGIS会沿用原始栅格尺寸导致输出像元不是正方形三是nodata_value写0只是把掩膜外的值置0栅格内部原本存在的989、999无数据值还要单独处理。3.3 内部无数据值989、999与0的三重筛查有同行问我为什么裁剪完以后唯一值统计表里出现989、999。原因是ESRI全球tif的NoData在不同分块里定义不一致。吉林省域tif如果跨越多个分块拼接时NoData值可能是989、999、0三种混在一起而ArcGIS一个栅格通常只认一个NoData值。处理办法是在裁剪前先做一次“复制栅格”把NoData设为统一值。如果已经裁完了就用“重分类”把989、999、0之外的真实值保留下来再把这几个值设为NoData。这里有个常见误用直接在栅格计算器里用Con(IsNull(raster), 0, raster)会把NoData变成0但0在ESRI分类体系里并不是有效类别统计时0就会混进去。正确做法是在重分类界面里把989、999、0这三类都勾选为NoData输出以后再用TabulateArea做统计并且确保工具选项里勾选“忽略NoData”这样内部云洞和外部边界不会参与像元计数。提示处理前先打开栅格属性面板查看“NoData值”那一栏记录原始tif的NoData定义再决定用哪个值做统一。别一上来就重采样先摸清底细再动手。4. 面积统计与专题图把10米栅格变成省域能用的结论4.1 不换投影面积到底会差多少吉林省全境面积大约18.7万平方公里。用Web墨卡托对全省做简单面积统计结果会比真实面积高出不少具体差异随纬度变化白城、松原一带纬度更高偏大更明显而Albers等积投影在43°N、46°N两条标准纬线之间面积变形能控制在很小的范围内对省级宏观评估完全够用。如果工作区域只限长春或吉林市这类中等范围也可以用CGCS2000高斯-克吕格3度带投影局部变形更小但一旦把范围扩大到全省跨带问题就来了。所以我对吉林省省域统计只用Albers对单一地市的精细图斑分析才换高斯投影。这一点决定面积统计结果的可靠度别偷懒。4.2 用ArcPy做分县面积统计从栅格像元到Excel表格做省域统计时最常干的事是把10米栅格叠加到县级行政区界上算出每个县每一类地类的面积。这里用TabulateArea面积制表工具它会按面要素对栅格像元累计输出一个dbf表。import arcpy lc_raster rD:\lc\jl_esri2021_clip.tif county_fc rD:\lc\jl_county.shp out_table rD:\lc\jl_county_area.dbf arcpy.env.overwriteOutput True # zone_field用行政区划代码保证唯一 arcpy.gp.TabulateArea( county_fc, PAC, lc_raster, out_table, 10 ) # 把统计表的单位从“像元数”换算成“平方公里” # 10米×10米像元 100平方米除以1e6得到km² arcpy.AddField_management(out_table, AREA_KM2, DOUBLE) fields arcpy.ListFields(out_table) count_fields [f.name for f in fields if f.name.startswith(C_)] with arcpy.da.UpdateCursor(out_table, [AREA_KM2] count_fields) as cursor: for row in cursor: total_km2 0.0 for i, v in enumerate(row[1:]): if v is not None: total_km2 v / 1e6 row[0] total_km2 cursor.updateRow(row)逻辑说明TabulateArea输出表的字段命名是“C_类别编码”值代表该县内这个类别的像元个数。dbf里的数值单位是像元不是面积。上面代码加了双精度字段把每个像元按100平方米换算成平方公里再汇总方便直接导进Excel做占比矩阵。参数说明第一个参数是区划面要素第二个是区划字段最好用稳定的行政区划代码PACcell_size传10是为了让统计工具明确栅格像元大小防止NoData混乱时像元数被重复计算。一个容易忽略的坑是TabulateArea对面要素覆盖范围的界定。当区划边界和栅格边界不完全对齐时边界像元按“中心点落入”原则归属因此统计结果和统计年鉴可能有百分之几的出入。这不代表方法错只能说明栅格化误差存在。做趋势比较可以做绝对值年报不建议直接引用。4.3 专题出图唯一值渲染与碎图斑处理统计完之后制图是硬需求。ESRI栅格直接加载进ArcGIS Pro时系统会默认按连续色带渲染必须改成“唯一值”。把1到9分别赋色水体用普蓝、树木用深绿、草地用中绿、洪水植被用青绿、农作物用橙黄、建设用地用暗红、裸地用土灰、冰雪用浅紫、云用白色。这套配色接近土地利用现状分类的习惯色出图时看图的人不需要翻图例就能认。但10米栅格直接出图会有大量椒盐状碎图斑全省缩到1比50万还凑合放大到1比5万就很花。我一般会做一个多数滤波Majority Filter再配合区域分组Region Group把面积小于1公顷的碎斑合并到相邻大图斑这步必须在面积统计之后做否则统计结果会偏离真实。制图归制图统计归统计两者用不同的数据副本避免互相污染。5. 吉林省10米数据的五个避坑记录从条带云到水稻田误分5.1 大片耕地在春播窗口被判成裸地现象4月中下旬吉林省中西部平原在ESRI分类里出现大面积7类裸地尤其梨树、农安一带外业看其实是等待播种的玉米田。原因ESRI全年合成影像会优先挑选云量少、反射清晰的时相而东北黑土区在整地播种期地表完全裸露模型很容易把无植被的农田归入裸地。解决不要单看2021年一帧。把2020、2022年分类结果叠加同一位置连续两年都是5类农作物就把它视为耕地。如果项目只做单年叠加4到9月的NDVI最大值合成峰值高的区域直接改判为农作物。5.2 东部林区的雪被误认为冰雪类现象长白山周边的和龙、安图、抚松一带有成片8类冰雪和2类树木镶嵌一眼看去像雪崩。原因高海拔林区冬季影像云雪难分Sentinel-2合成时保留了部分积雪帧模型把雪盖判成了冰雪类别。解决冰雪类在吉林省低海拔区几乎不存在所以直接用高程掩膜加类别替换把8类改判为相邻类别或按该区域占优类别重新投票。这里要保留ArcGIS的Con函数脚本做记录保证处理过程可复现。5.3 水稻田到底是洪水植被还是农作物现象吉林省水田多分布在东部低山盆地和西部灌区以延吉、珲春、前郭灌区为代表。同一块水稻田泡田插秧阶段被判为4类洪水植被成熟期又变成5类农作物两张专题图对不上。原因模型训练样本里对季节性淹没农田的标注不一致导致时间序列上同一地块编码波动。解决省级统计时把4类洪水植被中与5类农作物相邻的像元做一次时空一致性修正若某像元在两年内至少一年为5类则优先取5。也可以在水稻识别窗口用NDVI阈值覆盖分类结果把泡田期和成熟期的特征揉进去。5.4 989/999混进面积统计占比多出十几个点现象用ZonalHistogram统计吉林省各类面积各类面积之和比全省总面积大了10%以上。原因ESRI全球tif拼接后存在989、999这类无数据值它们没有统一被ArcGIS识别为NoData被当成真实类别参与了求和。解决统计前先做重分类把989、999、0都设为NoData再用TabulateArea并确认“忽略NoData”。怎么快速发现统计前打印唯一值表看最大值是否超过9一看到989就该知道问题在哪。5.5 碎图斑让地图没法看现象制图时图面密密麻麻的碎斑单像元100平方米的独立类别有几千处。原因10米分辨率对线性地物敏感田埂、沟渠、林带都被拆成单像元碎片分类模型对边缘像元本来就有振荡。解决制图副本上执行Majority Filter参数设八邻域多数权重再配合Region Group合并小图斑。注意这步会改变面积精度所以只用于制图副本不用于统计结果。5.6 通用检查顺序拿到新数据后先做这四步第一唯一值表确认类别编码有没有超过9的异常值。第二检查NoData值定义不一致就先统一。第三转等积投影再裁剪别拿Web墨卡托直接统计。第四做一次快速随机样点抽检每个地貌区抽20个点看高分影像确认大类对得上。这四步花不了半小时但能省掉后面所有返工的麻烦。6. 精度自检用随机样点混淆矩阵给ESRI数据算一次细账6.1 布点与判读精度自检不必全吉林铺开。我通常会按东、中、西三个地貌区分别抽取2到3个乡镇做验证区每区随机布100到200个点。判读底图用0.5米分辨率的高分影像吉林省林区参考二类调查小班边界耕地参考永久基本农田图斑。布点要避开边界像元尽量落在类别内部稳定区域。6.2 混淆矩阵与Kappa把ESRI分类结果和人工判读结果做成交叉表就是混淆矩阵。总体精度等于对角线之和除以总样本数Kappa系数等于总体精度减去期望精度再除以1减去期望精度。当总体精度低于75%时这张图只能做空间分布趋势参考不能用来报面积绝对值。要特别盯住“农作物-草地”“草地-洪水植被”两对混淆东北地区这两对是最容易打架的。6.3 用验证结果反哺后续流程我现在的习惯是每接一个省域底图任务先跑一遍自检把总体精度、Kappa和混淆严重的类别写进项目说明后面所有统计结论都注明“基于2021年ESRI 10米分类自检结果”。这个习惯救过我很多次。后来发现个别县水稻面积统计偏小不是计算流程错了而是底图把水田分到了洪水植被类拿自检记录一查就定位到了根因。坦诚地把底图分类误差写进报告比硬着头皮报数字更让甲方信任。希望帮到你。本文还有配套的精品资源点击获取
返回列表