简介这份资源是中国1km土地利用遥感监测数据集覆盖1980至2015年间共45个年份以1km×1km格网记录耕地、林地、草地、建设用地、水域等土地类型面向地理学、环境科学、城乡规划与经济学等领域的研究者及政策分析人员用于分析近四十年土地利用变迁、城市扩张、生态保护与农业动态。压缩包共126个文件约21.83MB以ArcInfo栅格格式为主包含adf、nit、dat等数据文件以及log日志、doc说明文档和xml元数据便于在GIS软件中直接读取与二次处理。目前已有2896人学习下载。数据时间序列连续、覆盖面广可用于对比不同年份识别耕地转建设用地、林地恢复等变化评估城市化速度与生态政策成效也可与社会经济数据融合构建土地利用变化模型为气候变化适应、灾害评估与可持续发展研究提供基础支撑。1. 拿到“中国1km土地利用数据1980-2015.zip”之后先想清楚三件事很多人第一次接触这套数据是被“1km”“1980-2015”“全国覆盖”这几个词吸引的——分辨率够用、时间跨度长、空间范围大看起来是做土地利用变化、生态评估、城市扩张分析的理想底料。但真正解压之后面对一堆按年份或按图幅命名的栅格文件第一反应往往是我该从哪一年开始看投影是什么分类体系怎么读值域代表什么如果这三个问题没搞清楚后面所有分析都是空中楼阁。这套数据本质上是一套长时间序列的全国土地利用栅格产品空间分辨率约1km时间覆盖1980年代到2015年前后通常按若干期发布。它适合做区域尺度的土地利用结构变化、耕地与建设用地转移、生态用地演变等分析不适合做地块级精细测绘。使用者主要是地理信息、生态、城市规划、农业资源方向的从业者和研究生。下面按“先立住理论、再动手复现”的顺序把从拿到压缩包到跑出第一张变化图的完整路径拆开讲。2. 解压之后先别急着打开搞清数据组织与分类体系2.1 常见目录结构与年份命名规律不同来源的打包方式略有差异但常见做法是按年份分文件夹每个文件夹里放若干分幅栅格或者直接放一张全国拼接好的大栅格。先列目录别急着用软件打开。# 查看压缩包内文件列表不直接解压 unzip -l 中国1km土地利用数据1980-2015.zip | head -50 # 解压到指定目录 mkdir -p ./landuse_data unzip 中国1km土地利用数据1980-2015.zip -d ./landuse_data # 查看解压后的目录树只看两层 find ./landuse_data -maxdepth 2 -type d | sort逻辑说明unzip -l先看清单避免解压出一堆不需要的中间文件。find只看目录层级快速判断是按年份还是按图幅组织。参数上-maxdepth 2控制递归深度防止输出过长。如果发现是按图幅分幅存放比如N50E110这类命名说明需要先拼接再分析。如果直接是landuse_1990.tif、landuse_2000.tif这种说明已经拼接好可以直接进入投影和分类检查。2.2 分类体系值域与含义必须逐一对上这套数据最常见的分类体系是6个一级类耕地、林地、草地、水域、建设用地、未利用地。但不同期次可能用不同编码比如1代表耕地、2代表林地也可能用10、20、30这种。不要凭经验猜一定要找到随数据附带的说明文件或者用栅格属性表反查。import rasterio import numpy as np # 打开某一年的栅格 with rasterio.open(./landuse_data/landuse_1990.tif) as src: print(投影:, src.crs) print(尺寸:, src.width, x, src.height) print(波段数:, src.count) print(数据类型:, src.dtypes) # 读取第一波段 arr src.read(1) # 统计唯一值和出现频次 unique, counts np.unique(arr, return_countsTrue) for u, c in zip(unique, counts): print(f值 {u}: {c} 像元)逻辑说明rasterio打开栅格后先看crs确认投影再看dtypes确认是整型还是浮点。np.unique统计值域分布如果出现0、1、2、3、4、5、6基本就是6类编码如果出现0、10、20、30就是另一套。参数上src.read(1)读第一波段土地利用数据通常只有单波段。注意如果唯一值里出现255或-9999通常是填充值或无效值后续分析前必须掩膜掉否则会污染统计结果。2.3 投影与坐标系别等面积算错了才回头1km分辨率的数据常见投影是Albers等面积投影或经纬度投影。如果是经纬度投影直接算面积会随纬度变化北方和南方的1km像元实际面积不一样。做面积统计前必须确认投影类型。with rasterio.open(./landuse_data/landuse_1990.tif) as src: print(src.crs) # 如果是地理坐标系像元大小单位是度 print(像元大小:, src.res) # 如果是投影坐标系像元大小单位是米如果src.crs显示的是EPSG:4326这类地理坐标系而src.res是0.01度左右说明需要先投影转换再算面积。常见做法是转到Albers等面积投影中央经线取105°E双标准纬线取25°N和47°N。这一步不做后面所有面积数字都只能当相对值看。3. 从单年读取到多期变化用Python跑通最小分析链路3.1 读取单年数据并做无效值掩膜先跑通一年确认分类和投影没问题再批量处理。下面这段代码读取1990年数据掩膜无效值统计各类面积占比。import rasterio import numpy as np import pandas as pd def read_landuse(path, nodata_values[0, 255, -9999]): with rasterio.open(path) as src: arr src.read(1).astype(np.float32) # 掩膜无效值 mask np.isin(arr, nodata_values) arr[mask] np.nan return arr, src.crs, src.transform arr, crs, transform read_landuse(./landuse_data/landuse_1990.tif) # 统计各类像元数 unique, counts np.unique(arr[~np.isnan(arr)], return_countsTrue) total counts.sum() for u, c in zip(unique, counts): print(f类别 {int(u)}: {c} 像元, 占比 {c/total*100:.2f}%)逻辑说明nodata_values根据实际情况调整常见的是0、255、-9999。np.isin做批量掩膜比循环快。统计时排除NaN避免无效值拉低占比。参数上astype(np.float32)是为了后续能存NaN整型不能直接存NaN。3.2 两期数据做变化检测转移矩阵怎么算变化检测最常用的输出是转移矩阵从1990年到2015年有多少耕地变成了建设用地多少林地变成了草地。用pandas.crosstab可以快速算。def change_matrix(path_t1, path_t2, nodata_values[0, 255, -9999]): arr1, _, _ read_landuse(path_t1, nodata_values) arr2, _, _ read_landuse(path_t2, nodata_values) # 只保留两期都有效的像元 valid (~np.isnan(arr1)) (~np.isnan(arr2)) a1 arr1[valid].astype(int) a2 arr2[valid].astype(int) # 计算转移矩阵 df pd.crosstab(a1, a2, rownames[1990], colnames[2015]) return df matrix change_matrix(./landuse_data/landuse_1990.tif, ./landuse_data/landuse_2015.tif) print(matrix)逻辑说明valid确保两期都有效才参与统计避免一期无效导致误判。crosstab默认按出现值建表如果某类在两期中都没出现表里不会显示需要手动补全类别。参数上astype(int)是为了让crosstab按整数分类而不是浮点。注意如果两期分类编码不一致比如一期用1-6另一期用10-60必须先做编码映射否则转移矩阵全是错位。3.3 批量处理多年数据循环与内存控制1980-2015通常有5到7期批量处理时不要一次性把所有栅格读进内存。逐期读取、逐期统计最后合并结果。import glob import os files sorted(glob.glob(./landuse_data/landuse_*.tif)) results [] for f in files: year os.path.basename(f).replace(landuse_, ).replace(.tif, ) arr, _, _ read_landuse(f) unique, counts np.unique(arr[~np.isnan(arr)], return_countsTrue) total counts.sum() for u, c in zip(unique, counts): results.append({year: year, class: int(u), count: c, pct: c/total*100}) df pd.DataFrame(results) print(df.pivot(indexyear, columnsclass, valuespct))逻辑说明glob按文件名排序确保年份顺序。每期读完只保留统计结果不保留原始数组内存占用可控。pivot输出年份×类别的占比表方便看趋势。参数上如果文件命名不规则需要手动指定年份列表。4. 避坑与排查这套数据最容易翻车的5个地方4.1 现象面积统计出来是经纬度度数不是平方公里原因数据是地理坐标系像元大小单位是度直接乘像元数得到的是“平方度”不是面积。解决先投影到等面积投影再统计像元数乘以像元面积比如1000m×1000m1km²。用rasterio.warp或gdalwarp做批量投影转换。4.2 现象两期数据叠加时像元对不齐转移矩阵出现大量错位原因不同期次可能用了不同投影或不同分辨率或者虽然投影相同但原点有偏移。解决用rasterio的reproject或gdalwarp统一到同一网格重采样方法用最近邻因为土地利用是分类数据不能用双线性。4.3 现象某一年数据打开全是0或255原因可能是分幅数据没拼接只打开了一个空白图幅或者无效值编码不是常见的0/255而是其他值。解决先看元数据里的nodata字段再用np.unique看实际值域。如果是分幅用gdal_merge.py或rasterio.merge拼接。4.4 现象分类编码在不同年份之间不一致原因不同期次可能用了不同的分类标准比如早期6类后期8类。解决建立映射表把后期类别合并到早期体系或者反过来。映射表要写进脚本注释别靠记忆。4.5 现象统计结果里出现小数类别比如1.5原因重采样时用了双线性或立方卷积导致分类值被插值成小数。解决分类数据重采样必须用最近邻。如果已经错了重新用最近邻跑一遍。5. 进阶技巧用变化热点图快速定位重点区域跑完转移矩阵之后下一步通常是找“哪里变化最剧烈”。一个实用技巧是计算变化强度图对每个像元统计1980-2015年间类别发生变化的次数变化次数越多颜色越深。def change_intensity(files, nodata_values[0, 255, -9999]): # 读取第一期作为基准 base, crs, transform read_landuse(files[0], nodata_values) intensity np.zeros_like(base, dtypenp.float32) prev base for f in files[1:]: curr, _, _ read_landuse(f, nodata_values) # 两期都有效且类别不同强度加1 changed (~np.isnan(prev)) (~np.isnan(curr)) (prev ! curr) intensity[changed] 1 prev curr return intensity, crs, transform files sorted(glob.glob(./landuse_data/landuse_*.tif)) intensity, crs, transform change_intensity(files) # 保存为GeoTIFF with rasterio.open(./change_intensity.tif, w, driverGTiff, heightintensity.shape[0], widthintensity.shape[1], count1, dtypenp.float32, crscrs, transformtransform) as dst: dst.write(intensity, 1)逻辑说明intensity初始为0每发现一次类别变化加1。prev逐期更新保证比较的是相邻两期。保存时保持原投影和变换方便在GIS里叠加。参数上dtypenp.float32足够存变化次数不需要双精度。拿到这张图之后可以叠加行政区划或生态保护区边界快速定位变化热点。我一般会再算一个“变化类型多样性”每个像元上出现过几种不同类别多样性高的区域往往是农牧交错带或城市边缘值得进一步下钻。最后说一个血泪经验这套数据的时间跨度大不同期次的精度和分类标准可能有细微差异做长时间序列分析时不要假设所有年份完全一致。先做一致性检查再下结论。希望帮到你。本文还有配套的精品资源点击获取