简介这份资源为2021年度辽宁省土地利用数据集由ESRI提供采用10米高精度栅格描绘农田、森林、建设用地、水域等土地用途坐标系为WGS84地理坐标系并已按中国大陆架整理、依省级行政区裁剪便于直接用于区域分析。它面向政策制定者、规划师、环境与农业研究人员以及需要开展土地利用规划、环境影响评估和城市发展趋势分析的用户属于中高级GIS应用场景。资源包整体约55.3MB以TIF格式栅格文件为主适合在ArcGIS等平台中加载、叠加与统计支持按省份进行空间比较与变化研究。目前已有134人学习下载可帮助读者快速获取权威、规范的土地覆盖基础数据减少自行裁剪与坐标转换的工作量为后续制图、建模和专题分析提供可靠底图。1. 拿到“2021年土地利用ESRI辽宁省10m精度”这份数据先搞清楚它能干什么如果你手头正好有一份“2021年土地利用ESRI辽宁省10m精度”的数据或者正打算去找这样一份数据来做辽宁省的国土空间分析、耕地变化监测、城市扩张研究那这篇文章就是写给你的。这份数据的核心信息其实就藏在标题里时间锁定2021年空间范围是辽宁省分辨率10米分类体系来自ESRI。它解决的是一个很实际的问题——在省级尺度上你需要一份足够细、足够新、又能直接和全球其他区域对比的土地利用底图而不是拿30米精度的全球产品去硬撑一个县域级别的分析。10米分辨率意味着什么意味着一个像素代表地面上10米×10米的方块辽宁省约14.8万平方公里算下来大概是14.8亿个像素。这个量级的数据普通笔记本打开会卡但放到GIS软件里做分区统计、做变化检测又刚好在可操作范围内。ESRI这套分类体系通常包含水体、林地、草地、耕地、建设用地、裸地等大类2021年的时效性对于“十四五”开局年的相关研究来说时间节点也对得上。适合谁用做辽宁省内市县国土空间规划前期分析的、做黑土地保护相关课题的、做辽河流域生态评估的以及需要一份基准年土地利用底图来叠加其他图层的人。不适合谁如果你要做的是地块级别的权属调查10米还是太粗如果你要的是逐月变化单一年份也不够。先想清楚这个边界再往下看怎么把它跑起来。2. 10米分辨率下的ESRI分类体系从像素值到地类名称的映射逻辑2.1 ESRI土地利用分类的编码规则与辽宁省的适配性ESRI这套土地利用数据底层是一个整型栅格每个像素值对应一个地类编码。常见的编码方案里1通常代表水体2代表林地4代表耕地5代表建设用地7代表草地8代表裸地11代表湿地等。不同年份或不同区域的产品在编码细节上可能有微调但整体框架是稳定的。辽宁省的特殊性在于它有漫长的海岸线、有大面积的辽东湾湿地、有辽西的低山丘陵、还有中部城市群带来的建设用地斑块。10米分辨率下这些地类的边界会比30米产品清晰很多尤其是建设用地和耕地的交错带不会像30米那样糊成一片。但这里有一个容易被忽略的点ESRI的分类体系是面向全球的它的“耕地”定义可能和国内二调、三调的“耕地”口径不完全一致。比如一些间作、套种的地块或者临时休耕的地块在ESRI体系里可能被归到草地或裸地。如果你直接拿这份数据去和统计年鉴的耕地面积对数字对不上是正常的。正确的做法是把它当作一份“空间分布参考”而不是“面积统计权威”。需要面积的时候用像素计数乘以100平方米10米×10米再根据你的研究目的决定要不要做地类归并。2.2 用Python读取栅格并统计各地类像素数拿到数据后第一步不是急着出图而是先确认像素值和地类的对应关系然后统计一遍各地类的像素数看看有没有明显的异常值。下面这段代码用rasterio读取栅格用numpy做统计不依赖ArcGIS在普通Python环境里就能跑。import rasterio import numpy as np from collections import Counter # 打开2021年辽宁省10米土地利用栅格 # 假设文件名为 liaoning_2021_landuse_10m.tif with rasterio.open(liaoning_2021_landuse_10m.tif) as src: print(栅格尺寸:, src.width, x, src.height) print(波段数:, src.count) print(坐标系:, src.crs) print(像素分辨率:, src.res) # 应该是 (10.0, 10.0) 左右 print(无效值:, src.nodata) # 读取第一波段 band src.read(1) # 统计每个像素值出现的次数 unique, counts np.unique(band, return_countsTrue) # 过滤掉无效值如果有nodata nodata src.nodata for val, cnt in zip(unique, counts): if nodata is not None and val nodata: continue # 面积 像素数 × 100 平方米再换算成公顷除以10000 area_ha cnt * 100 / 10000 print(f像素值 {val}: {cnt} 个像素, 约 {area_ha:.2f} 公顷)这段代码的逻辑很直接先确认栅格的基本元信息尤其是分辨率是不是真的10米、坐标系是不是投影坐标系。如果坐标系是地理坐标系比如WGS84那像素的“10米”只是名义上的实际面积会随纬度变化必须重投影到投影坐标系如UTM 51N或Albers等面积投影后再做面积统计。参数上src.res返回的是像素在地面单位的尺寸如果是投影坐标系单位通常是米如果是地理坐标系单位是度这时候直接算面积就是错的。src.nodata是无效值统计时要排除否则会把背景值也算成某个地类。跑完这段你会得到一张各地类的像素数和面积表。如果发现某个地类面积大得离谱比如“裸地”占了全省一半那大概率是编码映射搞错了或者数据里混入了云掩膜、阴影之类的值。这时候别急着往下做先回去核对编码表。3. 把10米栅格用起来从裁剪到重分类的完整操作链3.1 按行政区裁剪用GeoPandas和Rasterio做精确掩膜辽宁省的行政边界可以从公开的矢量数据里拿到比如从全国省级行政区划里筛出“辽宁省”。拿到边界后下一步是把全省的栅格裁剪到行政边界内去掉边界外的像素。这一步用rasterio的mask功能最方便配合geopandas读取矢量。import geopandas as gpd import rasterio from rasterio.mask import mask import numpy as np # 读取辽宁省行政边界矢量 # 假设文件为 liaoning_boundary.shp里面只有辽宁省一个要素 gdf gpd.read_file(liaoning_boundary.shp) # 确认坐标系一致如果不一致需要先投影转换 with rasterio.open(liaoning_2021_landuse_10m.tif) as src: if gdf.crs ! src.crs: gdf gdf.to_crs(src.crs) # 提取几何体 geoms gdf.geometry.values # 执行掩膜裁剪 out_image, out_transform mask(src, geoms, cropTrue, nodata0) out_meta src.meta.copy() # 更新元数据 out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: 0 }) # 写出裁剪后的栅格 with rasterio.open(liaoning_2021_landuse_10m_clipped.tif, w, **out_meta) as dest: dest.write(out_image)这里的关键参数是cropTrue它会把输出栅格的范围收紧到矢量边界的外接矩形减少数据量。nodata0指定裁剪后边界外的像素为0后续统计时要排除。注意如果矢量边界和栅格坐标系不一致必须先做to_crs转换否则掩膜结果会偏移甚至为空。另一个坑是有些行政边界矢量在海岸线附近和栅格的水体边界不完全重合裁剪后会出现一些细碎的无效像素这是正常的不影响整体分析。3.2 重分类把ESRI编码映射到你需要的地类体系裁剪完之后你拿到的还是ESRI的原始编码。如果你的研究只需要“耕地、林地、建设用地、水体”四大类就需要做一次重分类。用numpy的向量化操作比逐像素循环快得多。import rasterio import numpy as np # 定义ESRI编码到目标类别的映射 # 假设原始编码1水体, 2林地, 4耕地, 5建设用地, 7草地, 8裸地, 11湿地 # 目标类别1水体, 2林地, 3耕地, 4建设用地, 5其他 mapping { 1: 1, # 水体 - 水体 2: 2, # 林地 - 林地 4: 3, # 耕地 - 耕地 5: 4, # 建设用地 - 建设用地 7: 5, # 草地 - 其他 8: 5, # 裸地 - 其他 11: 5, # 湿地 - 其他 } with rasterio.open(liaoning_2021_landuse_10m_clipped.tif) as src: band src.read(1) meta src.meta.copy() # 创建输出数组默认0无效值 reclass np.zeros_like(band, dtypenp.uint8) # 逐类映射 for old_val, new_val in mapping.items(): reclass[band old_val] new_val # 保留原始nodata为0 reclass[band 0] 0 meta.update(dtyperasterio.uint8, nodata0) with rasterio.open(liaoning_2021_landuse_reclass.tif, w, **meta) as dest: dest.write(reclass, 1)这段代码的核心是reclass[band old_val] new_val它利用numpy的布尔索引一次性完成所有同类像素的赋值比Python循环快几个数量级。参数上dtype改成uint8是因为目标类别只有5类用8位无符号整数足够能显著减小文件体积。nodata0保持不变确保后续统计时能正确排除无效区域。如果你需要保留更多类别比如把湿地单独列出来只需要在mapping里给11映射一个新值同时把dtype改成uint16。重分类之后建议再跑一次第2章里的统计代码确认各地类面积比例合理。比如辽宁省的耕地占比通常在30%左右林地占比也在30%上下建设用地占比在10%以内。如果重分类后建设用地突然占了30%那肯定是映射搞反了。4. 避坑与排查10米土地利用数据在辽宁省场景下的5个血泪教训4.1 现象面积统计结果和统计年鉴对不上差出好几倍原因最常见的是坐标系问题。如果栅格是地理坐标系WGS84像素的“10米”是名义值实际面积随纬度变化在辽宁省北纬38°到43°误差能达到20%以上。另一个原因是nodata值没排除干净把背景值当成了某个地类。解决先用src.crs确认坐标系如果是地理坐标系用gdalwarp或rasterio的calculate_default_transform重投影到UTM 51N辽宁省大部分区域适用或Albers等面积投影。重投影后再统计面积。nodata值在统计前用band[band nodata] 0统一处理。4.2 现象裁剪后边界外出现大量0值像素统计时被算成“裸地”原因裁剪时nodata0但后续重分类或统计时忘了排除0。有些统计脚本直接把0当成某个地类编码导致边界外区域被计入。解决在任何统计或重分类之前先执行band[band 0] np.nan或者显式排除0。如果用的是rasterio的mask功能输出栅格的nodata已经设为0读取后直接band[band src.nodata] 0即可。4.3 现象建设用地和耕地边界模糊10米分辨率下仍然糊成一片原因ESRI的全球分类模型在城乡交错带的表现有限10米分辨率虽然比30米好但混合像素问题依然存在。一个10米像素里可能一半是宅基地一半是菜地模型只能给一个标签。解决如果研究重点是城乡交错带不要依赖单一像素的分类结果。常见做法是结合更高分辨率影像如Sentinel-2的10米多光谱做局部校验或者用3×3窗口的众数滤波平滑边界。但滤波会损失细节用之前想清楚你的分析尺度。4.4 现象用ArcGIS打开栅格显示正常用Python读取后像素值全变了原因ArcGIS可能自动应用了拉伸或色彩映射你看到的颜色对应的像素值和你以为的不一样。另外如果栅格有多个波段Python默认读第一波段而ArcGIS可能显示的是RGB合成。解决在Python里先用src.read()读所有波段看src.count是几。如果是单波段src.read(1)就是分类编码。如果是多波段确认哪个波段是分类结果。另外用src.colormap(1)可以查看色彩映射表确认像素值和颜色的对应关系。4.5 现象重分类后文件体积暴涨从几百MB变成几个GB原因原始栅格可能是压缩的整型重分类时如果用了float32或int32文件体积会翻好几倍。另外如果输出时没有指定压缩选项默认不压缩。解决重分类时dtype选uint8或uint16写出时加上compresslzw或compressdeflate。在rasterio里meta.update(compresslzw)即可。如果数据量实在太大可以考虑分块处理用window参数逐块读写。5. 进阶技巧用10米数据做辽宁省耕地变化的快速验证5.1 单年份数据的验证思路和公开产品做交叉比对你手里只有2021年一年没法直接做变化检测。但可以用它来验证其他年份的数据质量。比如找一份2020年的30米全球土地利用产品重采样到10米然后和你的2021年数据做空间交叉表。如果两个产品在耕地分布上的一致性超过85%说明你的2021年数据在耕地这一类的空间格局上是可信的。如果一致性低于70%要么是分类体系差异太大要么是某一方的数据有问题。交叉表的做法很简单把两个栅格对齐到同一网格然后统计每个像素对的出现次数。用pandas的crosstab就能出表。import rasterio import numpy as np import pandas as pd # 读取两个已对齐的栅格 with rasterio.open(liaoning_2021_landuse_reclass.tif) as src1: arr1 src1.read(1).flatten() with rasterio.open(liaoning_2020_landuse_30m_resampled.tif) as src2: arr2 src2.read(1).flatten() # 排除无效值 mask (arr1 0) (arr2 0) arr1_valid arr1[mask] arr2_valid arr2[mask] # 构建交叉表 df pd.DataFrame({2021_10m: arr1_valid, 2020_30m: arr2_valid}) ct pd.crosstab(df[2021_10m], df[2020_30m]) print(ct) # 计算耕地假设编码3的一致性 if 3 in ct.index and 3 in ct.columns: 耕地一致 ct.loc[3, 3] 耕地总数_2021 ct.loc[3, :].sum() 一致性 耕地一致 / 耕地总数_2021 print(f耕地空间一致性: {一致性:.2%})这段代码的关键是先把两个栅格对齐到同一网格否则flatten后的数组长度不一样没法直接比较。对齐可以用rasterio.warp.reproject或者先用ArcGIS的“重采样”工具。交叉表出来之后重点看对角线上的数字对角线占比越高说明两个产品越一致。如果某个地类的一致性特别低比如草地对裸地那说明两个产品的分类边界在这类地物上分歧很大用的时候要小心。5.2 我自己的习惯先做小区域试跑再推全省辽宁省东西跨度大辽东和辽西的地貌差异明显。我一般不会一上来就跑全省而是先选两个典型县——比如一个辽中平原的农业县一个辽西山区的林业县——分别裁剪、统计、出图。如果这两个县的统计结果合理地类边界和影像对得上再推全省。这样做的好处是万一编码映射错了或者坐标系有问题在小区域上就能发现不用等全省跑完再返工。另外10米数据在县级尺度上做制图出图比例尺建议控制在1:50000到1:100000之间。再放大像素的锯齿感就出来了再缩小10米的优势又体现不出来。出图时用matplotlib的imshow配合interpolationnearest保持像素边界清晰不要用双线性插值否则地类边界会糊掉。希望帮到你。本文还有配套的精品资源点击获取