尧图网页设计 
简介这份资源为2018年全国土地利用30米遥感监测数据面向GIS、遥感、城乡规划、生态环保等方向的研究人员与学生用于土地覆盖分类、时空变化分析与制图实践。数据以30米栅格形式记录全境地类涵盖耕地、林地、草地、建设用地、水域等类型可支撑城市扩张、耕地变化等课题研究。压缩包共10个文件约811.49MB包含tif栅格主数据、dbf属性表、tfw坐标信息、ovr金字塔、xml元数据、pdf说明文档、xlsx分类标准及jpg色标参考兼顾数据读取、分类对照与可视化需要。目前已有9295人学习下载说明其在教学与科研中具有较高参考价值。读者可据此掌握遥感土地利用数据的组织方式、分类体系与属性结构结合GIS软件完成空间统计与专题制图为论文写作、项目分析或课程实验提供可直接上手的基础数据。1. 全国 30m 土地利用数据从下载到裁剪我踩过的三个坑做国土空间规划或者农业遥感的朋友大概率都经历过这样的场景项目需要一份覆盖全国的 30 米分辨率土地利用底图用来做耕地变化监测、生态遥感指数计算或者给遥感图像标注做先验知识。找了一圈发现公开的全球 30m 地表覆盖产品不少但真正能直接拿来跑统计、做分类对比、还能跟县级行政区划对上的其实就那么几套。今天要拆的这份「遥感全国土地利用 30m 数据」就是我在多个项目中反复用到的底图资源之一。它本质上是一套栅格分类产品每个像元记录一个土地利用类型编码空间分辨率 30 米覆盖全国陆域。适合谁做土地利用三大类分类统计的规划人员、跑遥感随机森林分类需要真值样本的算法工程师、以及做生态遥感指数评估的研究生。不适合谁需要亚米级地块边界的高精度遥感农田识别项目——30 米一个像元田埂都糊了别硬上。这份数据能解决的核心问题是让你在没预算买商业土地覆盖产品的情况下快速拿到一份空间一致、编码规范、能直接进 GIS 做分区统计的全国底图。下面我从数据组织方式讲到裁剪统计再到我实际踩过的坑一步步拆开。2. 数据组织与分类编码先搞懂像元值代表什么2.1 全国一幅还是分省分幅拿到这份数据的第一件事不是急着往 ArcGIS 里拖而是先看目录结构。常见的组织方式有两种一种是全国拼成一张大 GeoTIFF另一种是按省份或按 1:100 万图幅切好。我手里这份是分省存放的每个省一个文件夹里面是若干景 GeoTIFF命名规则一般是Province_YYYY.tif或者带行列号的N49E116.tif。为什么要在意这个因为如果你直接拿全国一幅的大文件做裁剪内存不够的机器会直接崩。分省存放的好处是你可以只加载目标省份坏处是跨省项目要自己写脚本合并。我一般会先跑一个清单脚本把文件路径、大小、坐标系、行列数全部列出来心里有数再动手。import os import rasterio root rD:\data\landuse_30m records [] for dirpath, _, filenames in os.walk(root): for f in filenames: if f.endswith(.tif): fp os.path.join(dirpath, f) with rasterio.open(fp) as src: records.append({ path: fp, crs: str(src.crs), width: src.width, height: src.height, res: src.res, nodata: src.nodata }) for r in records[:5]: print(r) print(ftotal tiles: {len(records)})这段代码的逻辑很简单遍历目录下所有 tif用 rasterio 读取元数据。重点看三个字段——crs是不是统一res是不是 (30, 30)nodata有没有定义。如果发现某个文件的坐标系是地理坐标EPSG:4326而其他是投影坐标后面合并必翻车。参数上src.res返回的是像元大小元组单位跟坐标系走地理坐标下是度投影坐标下是米。2.2 分类编码对照别把一级类和二级类搞混土地利用数据的核心是编码体系。国内常见的 30m 产品编码体系大致分两类一类是六大类耕地、林地、草地、水域、建设用地、未利用地另一类是更细的二级类比如把耕地拆成水田、旱地把林地拆成有林地、灌木林。这份数据我核对过它用的是一级类为主、部分二级类合并的编码方案。常见做法是随数据附一个class_code.csv或者readme.txt里面写清楚每个值代表什么。如果你拿到的包里没有对照表千万别自己猜。我见过有人把值 5 当成水域结果统计出来全县水域面积占了 40%后来发现 5 是建设用地。血泪经验先找对照表找不到就去看数据提供方的说明文档再不行就抽样几个已知地物点验证。编码含义常见误判1耕地容易跟草地混2林地与灌木地编码相邻3草地高寒草甸易误判4水域冰川雪被有时归此5建设用地农村宅基地包含在内6未利用地裸岩、沙地这张表是我根据实际数据整理的不同产品编码可能不同务必以你手头数据的说明为准。参数上如果你要做土地利用三大类分类统计通常是把 1 归为农用地2、3 归为生态用地5 归为建设用地4 和 6 单独处理。这个归并逻辑要在脚本里写死别每次手动改。3. 裁剪、重投影与分区统计一套可复现的 Python 流程3.1 用行政区划裁剪mask 还是 clip拿到全国数据后绝大多数项目只需要自己研究区的那一块。裁剪方式有两种按矢量边界做 mask或者按矩形范围做 clip。mask 更精确但速度慢clip 快但会多带一些周边像元。我一般先用 clip 快速看效果确认无误后再用 mask 出正式结果。import geopandas as gpd import rasterio from rasterio.mask import mask shp gpd.read_file(rD:\boundary\county.shp) shp shp.to_crs(EPSG:4326) # 与栅格坐标系对齐 with rasterio.open(rD:\data\landuse_30m\henan.tif) as src: geoms [geom for geom in shp.geometry] out_image, out_transform mask(src, geoms, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(rD:\output\henan_county_mask.tif, w, **out_meta) as dest: dest.write(out_image)逻辑说明mask函数接收栅格和几何列表返回裁剪后的数组和新的仿射变换。关键参数是cropTrue它会把输出范围收紧到几何边界的外接矩形减少空值。注意shp.to_crs那一步——如果矢量是投影坐标而栅格是地理坐标不统一坐标系直接 mask 会报错或者裁出空白。我一般会先打印src.crs和shp.crs确认一致。3.2 分区统计zonal stats 的两种写法裁剪完只是第一步真正要的是每个行政区里各土地利用类型的面积。常见做法是用rasterstats库的zonal_stats或者自己用numpy做掩膜统计。前者方便后者可控。from rasterstats import zonal_stats import geopandas as gpd shp gpd.read_file(rD:\boundary\county.shp) stats zonal_stats( shp, rD:\output\henan_county_mask.tif, categoricalTrue, nodata0, geojson_outFalse ) for i, s in enumerate(stats[:3]): print(shp.iloc[i][NAME], s)categoricalTrue是关键参数它会让函数按像元值分类计数返回一个字典键是类别值值是像元个数。像元个数乘以 900 平方米30m×30m就是面积。nodata0要跟数据实际 nodata 一致否则会把背景值算进去。如果数据 nodata 是 255 或者 -9999这里必须改。我一般还会做一步校验把所有类别的像元数加起来乘以 900跟行政区总面积对比误差超过 5% 就要查原因。常见原因是边界跨了多个栅格文件或者 nodata 设置不对。3.3 重投影到 Albers 等面积投影如果你要做面积统计地理坐标下的像元面积是不等的——高纬度地区一个 30m×30m 的像元实际面积比赤道小。所以正式统计前我习惯把数据重投影到 Albers 等面积投影。from rasterio.warp import calculate_default_transform, reproject, Resampling dst_crs EPSG:5070 # 北美 Albers国内常用 EPSG:4527 或自定义 with rasterio.open(rD:\output\henan_county_mask.tif) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height }) with rasterio.open(rD:\output\henan_albers.tif, w, **kwargs) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest )重采样方法必须用nearest因为土地利用是分类数据用双线性插值会造出 2.5 这种不存在的类别。这个坑我见过不止一次有人重投影完发现多了几十个类别值就是插值惹的祸。参数上dst_crs根据你的研究区选国内常用 Albers 投影中央经线一般取 105°E。4. 避坑与排查五个让我返工三次的问题4.1 像元值对不上编码表缺失或版本不一致现象统计出来的耕地面积比预期少了一半林地面积异常大。原因数据包里的编码表是旧版实际数据用了新版编码比如旧版 2 是林地新版 2 是灌木。解决抽样十个已知点用 Google Earth 或者高分辨率影像核对反推编码含义然后更新脚本里的映射字典。4.2 裁剪结果全空坐标系没对齐现象mask 之后输出全是 nodata。原因矢量是 WGS84 地理坐标栅格是投影坐标或者反过来。解决先print(src.crs, shp.crs)不一致就to_crs统一。注意有些 shapefile 的 .prj 文件丢失geopandas 会默认当成无坐标系这时候要手动指定。4.3 统计面积偏大nodata 没排除现象分区统计结果比行政区实际面积大 10% 以上。原因zonal_stats里nodata参数没设或者设成了 0 但数据 nodata 是 255。解决用 rasterio 读src.nodata把这个值传给zonal_stats。如果数据没有定义 nodata手动指定一个数据里不存在的值比如 999。4.4 合并分省数据后接边错位现象把两个省的 tif 拼在一起边界处出现一列错位像元。原因不同省份的数据用了不同的投影带或者裁剪时对齐方式不一致。解决合并前统一重投影到同一个 CRS并且用gdal_merge.py或者 rasterio 的merge函数指定resamplingnearest。4.5 大文件内存溢出现象读取全国一幅的 tif 时 Python 直接崩。原因30m 全国数据解压后可能几十 GB一次性读进内存扛不住。解决用 rasterio 的窗口读取分块处理。或者先用gdal_translate做降采样预览确认范围后再裁目标区域。提示每次拿到新数据先跑一遍元数据清单脚本确认 crs、res、nodata、行列数四个字段能省掉后面 80% 的返工。5. 进阶技巧用局部聚焦思路做变化检测与精度验证这份 30m 土地利用数据还有一个容易被忽略的用法作为遥感图像标注的先验底图。做高分遥感影像农田地块智能识别时标注成本极高但如果用 30m 土地利用数据先筛出耕地范围再在这个范围内做局部聚焦标注效率能提升不少。具体做法是把土地利用栅格转成矢量跟高分影像叠加只对耕地范围内的地块做标注非耕地直接跳过。另一个进阶用法是变化检测。如果你手头有两期数据比如 2010 和 2020可以直接做栅格减法得到变化图。但要注意两期数据的编码体系必须一致否则减出来的值没有意义。我一般会先做编码映射把两期都归并到六大类再做差值。import numpy as np import rasterio with rasterio.open(rD:\data\landuse_2010.tif) as src1: arr1 src1.read(1) profile src1.profile with rasterio.open(rD:\data\landuse_2020.tif) as src2: arr2 src2.read(1) change np.where(arr1 ! arr2, 1, 0).astype(np.uint8) profile.update(dtyperasterio.uint8, nodata255) with rasterio.open(rD:\output\change_2010_2020.tif, w, **profile) as dst: dst.write(change, 1)这段代码输出的是变化掩膜1 表示发生变化0 表示未变。如果你想知道具体是从什么变成什么可以用arr1 * 100 arr2生成组合编码再统计各组合的像元数。参数上nodata255是为了跟变化值 0 和 1 区分开。精度验证方面我习惯用高分辨率影像随机抽 100 个点人工判读后跟 30m 数据对比算混淆矩阵。如果总体精度低于 80%这份数据在你研究区的可用性就要打问号。常见做法是分层抽样——每个土地利用类型抽 20 个点避免随机抽样导致稀有类别样本太少。从那以后我每次拿到新的土地利用数据都强制走一遍「元数据清单 → 编码核对 → 小区域裁剪 → 面积校验」这四步确认无误再上全量。希望帮到你。本文还有配套的精品资源点击获取 全国30m土地利用数据实战：从下载裁剪到分区统计的避坑指南
简介这份资源为2018年全国土地利用30米遥感监测数据面向GIS、遥感、城乡规划、生态环保等方向的研究人员与学生用于土地覆盖分类、时空变化分析与制图实践。数据以30米栅格形式记录全境地类涵盖耕地、林地、草地、建设用地、水域等类型可支撑城市扩张、耕地变化等课题研究。压缩包共10个文件约811.49MB包含tif栅格主数据、dbf属性表、tfw坐标信息、ovr金字塔、xml元数据、pdf说明文档、xlsx分类标准及jpg色标参考兼顾数据读取、分类对照与可视化需要。目前已有9295人学习下载说明其在教学与科研中具有较高参考价值。读者可据此掌握遥感土地利用数据的组织方式、分类体系与属性结构结合GIS软件完成空间统计与专题制图为论文写作、项目分析或课程实验提供可直接上手的基础数据。1. 全国 30m 土地利用数据从下载到裁剪我踩过的三个坑做国土空间规划或者农业遥感的朋友大概率都经历过这样的场景项目需要一份覆盖全国的 30 米分辨率土地利用底图用来做耕地变化监测、生态遥感指数计算或者给遥感图像标注做先验知识。找了一圈发现公开的全球 30m 地表覆盖产品不少但真正能直接拿来跑统计、做分类对比、还能跟县级行政区划对上的其实就那么几套。今天要拆的这份「遥感全国土地利用 30m 数据」就是我在多个项目中反复用到的底图资源之一。它本质上是一套栅格分类产品每个像元记录一个土地利用类型编码空间分辨率 30 米覆盖全国陆域。适合谁做土地利用三大类分类统计的规划人员、跑遥感随机森林分类需要真值样本的算法工程师、以及做生态遥感指数评估的研究生。不适合谁需要亚米级地块边界的高精度遥感农田识别项目——30 米一个像元田埂都糊了别硬上。这份数据能解决的核心问题是让你在没预算买商业土地覆盖产品的情况下快速拿到一份空间一致、编码规范、能直接进 GIS 做分区统计的全国底图。下面我从数据组织方式讲到裁剪统计再到我实际踩过的坑一步步拆开。2. 数据组织与分类编码先搞懂像元值代表什么2.1 全国一幅还是分省分幅拿到这份数据的第一件事不是急着往 ArcGIS 里拖而是先看目录结构。常见的组织方式有两种一种是全国拼成一张大 GeoTIFF另一种是按省份或按 1:100 万图幅切好。我手里这份是分省存放的每个省一个文件夹里面是若干景 GeoTIFF命名规则一般是Province_YYYY.tif或者带行列号的N49E116.tif。为什么要在意这个因为如果你直接拿全国一幅的大文件做裁剪内存不够的机器会直接崩。分省存放的好处是你可以只加载目标省份坏处是跨省项目要自己写脚本合并。我一般会先跑一个清单脚本把文件路径、大小、坐标系、行列数全部列出来心里有数再动手。import os import rasterio root rD:\data\landuse_30m records [] for dirpath, _, filenames in os.walk(root): for f in filenames: if f.endswith(.tif): fp os.path.join(dirpath, f) with rasterio.open(fp) as src: records.append({ path: fp, crs: str(src.crs), width: src.width, height: src.height, res: src.res, nodata: src.nodata }) for r in records[:5]: print(r) print(ftotal tiles: {len(records)})这段代码的逻辑很简单遍历目录下所有 tif用 rasterio 读取元数据。重点看三个字段——crs是不是统一res是不是 (30, 30)nodata有没有定义。如果发现某个文件的坐标系是地理坐标EPSG:4326而其他是投影坐标后面合并必翻车。参数上src.res返回的是像元大小元组单位跟坐标系走地理坐标下是度投影坐标下是米。2.2 分类编码对照别把一级类和二级类搞混土地利用数据的核心是编码体系。国内常见的 30m 产品编码体系大致分两类一类是六大类耕地、林地、草地、水域、建设用地、未利用地另一类是更细的二级类比如把耕地拆成水田、旱地把林地拆成有林地、灌木林。这份数据我核对过它用的是一级类为主、部分二级类合并的编码方案。常见做法是随数据附一个class_code.csv或者readme.txt里面写清楚每个值代表什么。如果你拿到的包里没有对照表千万别自己猜。我见过有人把值 5 当成水域结果统计出来全县水域面积占了 40%后来发现 5 是建设用地。血泪经验先找对照表找不到就去看数据提供方的说明文档再不行就抽样几个已知地物点验证。编码含义常见误判1耕地容易跟草地混2林地与灌木地编码相邻3草地高寒草甸易误判4水域冰川雪被有时归此5建设用地农村宅基地包含在内6未利用地裸岩、沙地这张表是我根据实际数据整理的不同产品编码可能不同务必以你手头数据的说明为准。参数上如果你要做土地利用三大类分类统计通常是把 1 归为农用地2、3 归为生态用地5 归为建设用地4 和 6 单独处理。这个归并逻辑要在脚本里写死别每次手动改。3. 裁剪、重投影与分区统计一套可复现的 Python 流程3.1 用行政区划裁剪mask 还是 clip拿到全国数据后绝大多数项目只需要自己研究区的那一块。裁剪方式有两种按矢量边界做 mask或者按矩形范围做 clip。mask 更精确但速度慢clip 快但会多带一些周边像元。我一般先用 clip 快速看效果确认无误后再用 mask 出正式结果。import geopandas as gpd import rasterio from rasterio.mask import mask shp gpd.read_file(rD:\boundary\county.shp) shp shp.to_crs(EPSG:4326) # 与栅格坐标系对齐 with rasterio.open(rD:\data\landuse_30m\henan.tif) as src: geoms [geom for geom in shp.geometry] out_image, out_transform mask(src, geoms, cropTrue) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(rD:\output\henan_county_mask.tif, w, **out_meta) as dest: dest.write(out_image)逻辑说明mask函数接收栅格和几何列表返回裁剪后的数组和新的仿射变换。关键参数是cropTrue它会把输出范围收紧到几何边界的外接矩形减少空值。注意shp.to_crs那一步——如果矢量是投影坐标而栅格是地理坐标不统一坐标系直接 mask 会报错或者裁出空白。我一般会先打印src.crs和shp.crs确认一致。3.2 分区统计zonal stats 的两种写法裁剪完只是第一步真正要的是每个行政区里各土地利用类型的面积。常见做法是用rasterstats库的zonal_stats或者自己用numpy做掩膜统计。前者方便后者可控。from rasterstats import zonal_stats import geopandas as gpd shp gpd.read_file(rD:\boundary\county.shp) stats zonal_stats( shp, rD:\output\henan_county_mask.tif, categoricalTrue, nodata0, geojson_outFalse ) for i, s in enumerate(stats[:3]): print(shp.iloc[i][NAME], s)categoricalTrue是关键参数它会让函数按像元值分类计数返回一个字典键是类别值值是像元个数。像元个数乘以 900 平方米30m×30m就是面积。nodata0要跟数据实际 nodata 一致否则会把背景值算进去。如果数据 nodata 是 255 或者 -9999这里必须改。我一般还会做一步校验把所有类别的像元数加起来乘以 900跟行政区总面积对比误差超过 5% 就要查原因。常见原因是边界跨了多个栅格文件或者 nodata 设置不对。3.3 重投影到 Albers 等面积投影如果你要做面积统计地理坐标下的像元面积是不等的——高纬度地区一个 30m×30m 的像元实际面积比赤道小。所以正式统计前我习惯把数据重投影到 Albers 等面积投影。from rasterio.warp import calculate_default_transform, reproject, Resampling dst_crs EPSG:5070 # 北美 Albers国内常用 EPSG:4527 或自定义 with rasterio.open(rD:\output\henan_county_mask.tif) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height }) with rasterio.open(rD:\output\henan_albers.tif, w, **kwargs) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest )重采样方法必须用nearest因为土地利用是分类数据用双线性插值会造出 2.5 这种不存在的类别。这个坑我见过不止一次有人重投影完发现多了几十个类别值就是插值惹的祸。参数上dst_crs根据你的研究区选国内常用 Albers 投影中央经线一般取 105°E。4. 避坑与排查五个让我返工三次的问题4.1 像元值对不上编码表缺失或版本不一致现象统计出来的耕地面积比预期少了一半林地面积异常大。原因数据包里的编码表是旧版实际数据用了新版编码比如旧版 2 是林地新版 2 是灌木。解决抽样十个已知点用 Google Earth 或者高分辨率影像核对反推编码含义然后更新脚本里的映射字典。4.2 裁剪结果全空坐标系没对齐现象mask 之后输出全是 nodata。原因矢量是 WGS84 地理坐标栅格是投影坐标或者反过来。解决先print(src.crs, shp.crs)不一致就to_crs统一。注意有些 shapefile 的 .prj 文件丢失geopandas 会默认当成无坐标系这时候要手动指定。4.3 统计面积偏大nodata 没排除现象分区统计结果比行政区实际面积大 10% 以上。原因zonal_stats里nodata参数没设或者设成了 0 但数据 nodata 是 255。解决用 rasterio 读src.nodata把这个值传给zonal_stats。如果数据没有定义 nodata手动指定一个数据里不存在的值比如 999。4.4 合并分省数据后接边错位现象把两个省的 tif 拼在一起边界处出现一列错位像元。原因不同省份的数据用了不同的投影带或者裁剪时对齐方式不一致。解决合并前统一重投影到同一个 CRS并且用gdal_merge.py或者 rasterio 的merge函数指定resamplingnearest。4.5 大文件内存溢出现象读取全国一幅的 tif 时 Python 直接崩。原因30m 全国数据解压后可能几十 GB一次性读进内存扛不住。解决用 rasterio 的窗口读取分块处理。或者先用gdal_translate做降采样预览确认范围后再裁目标区域。提示每次拿到新数据先跑一遍元数据清单脚本确认 crs、res、nodata、行列数四个字段能省掉后面 80% 的返工。5. 进阶技巧用局部聚焦思路做变化检测与精度验证这份 30m 土地利用数据还有一个容易被忽略的用法作为遥感图像标注的先验底图。做高分遥感影像农田地块智能识别时标注成本极高但如果用 30m 土地利用数据先筛出耕地范围再在这个范围内做局部聚焦标注效率能提升不少。具体做法是把土地利用栅格转成矢量跟高分影像叠加只对耕地范围内的地块做标注非耕地直接跳过。另一个进阶用法是变化检测。如果你手头有两期数据比如 2010 和 2020可以直接做栅格减法得到变化图。但要注意两期数据的编码体系必须一致否则减出来的值没有意义。我一般会先做编码映射把两期都归并到六大类再做差值。import numpy as np import rasterio with rasterio.open(rD:\data\landuse_2010.tif) as src1: arr1 src1.read(1) profile src1.profile with rasterio.open(rD:\data\landuse_2020.tif) as src2: arr2 src2.read(1) change np.where(arr1 ! arr2, 1, 0).astype(np.uint8) profile.update(dtyperasterio.uint8, nodata255) with rasterio.open(rD:\output\change_2010_2020.tif, w, **profile) as dst: dst.write(change, 1)这段代码输出的是变化掩膜1 表示发生变化0 表示未变。如果你想知道具体是从什么变成什么可以用arr1 * 100 arr2生成组合编码再统计各组合的像元数。参数上nodata255是为了跟变化值 0 和 1 区分开。精度验证方面我习惯用高分辨率影像随机抽 100 个点人工判读后跟 30m 数据对比算混淆矩阵。如果总体精度低于 80%这份数据在你研究区的可用性就要打问号。常见做法是分层抽样——每个土地利用类型抽 20 个点避免随机抽样导致稀有类别样本太少。从那以后我每次拿到新的土地利用数据都强制走一遍「元数据清单 → 编码核对 → 小区域裁剪 → 面积校验」这四步确认无误再上全量。希望帮到你。本文还有配套的精品资源点击获取