简介这份资源为新疆维吾尔自治区30m精度地形地貌栅格数据面向地理信息、遥感、水文、生态及国土空间规划等方向的研究者与工程人员用于区域地貌分类、地形起伏分析与建模制图。数据依据海拔划分为低海拔、中海拔、中高海拔、高海拔、极高海拔依据起伏程度区分丘陵、小起伏、中起伏、大起伏、极大起伏并涵盖丘陵、山脉、平原、沟壑等地形类型以及海积、湖积、冲积、洪积、风积、冰碛、剥蚀等成因类型坐标采用WGS84与Albers_Conic_Equal_Area已按省整理为tif格式。压缩包约37.79MB文件总数暂未提供主要文件类型为tif栅格数据可直接在ArcGIS、QGIS等平台加载使用。目前已有175人学习下载适合需要快速获取新疆地貌基础数据、开展空间分析与专题制图的读者参考使用。1. 新疆30m地形地貌数据从拿到压缩包到跑通分析中间隔着什么如果你手头正好有一个「新疆维吾尔自治区地形地貌最新30m精度.rar」大概率不是单纯想收藏一份数据而是想把它用起来——做地形起伏度统计、坡度坡向提取、地貌类型分区或者给水文、生态、选址类模型当底图。新疆地域辽阔从阿尔泰山到昆仑山从塔里木盆地到吐鲁番洼地高差跨度极大30m 分辨率意味着单省数据量就可能达到几十 GB 级别直接双击打开往往不是正确姿势。这份数据通常以栅格 DEM 或地貌分类栅格的形式存在可能附带投影信息、NoData 定义和分层打包结构。真正要解决的问题是怎么解压不炸盘、怎么确认坐标系、怎么裁出研究区、怎么让 GIS 软件和 Python 都能稳定读取。适合做区域地理分析、遥感建模、工程选址和教学复现的从业者也适合刚接触大栅格数据、想少走弯路的新手。下面按我实际处理这类数据的顺序把每一步拆开讲。2. 先搞清楚压缩包里装的是什么30m 地形地貌数据的结构判断2.1 从文件名和目录树反推数据组织方式拿到一个 .rar 包第一件事不是解压到桌面而是先看目录树。Windows 下可以用 7-Zip 或 WinRAR 的「查看」功能Linux/macOS 下用unrar l列出内容。常见结构有三种第一种是单幅大栅格比如xinjiang_dem_30m.tif整省一张图第二种是按图幅分块比如N40E080.tif这种百万图幅命名第三种是按地州或流域分文件夹每个文件夹里再放若干 tif。判断依据很简单看文件数量和单个文件大小。如果只有一个几 GB 的 tif说明是整省拼接如果有几百个几十 MB 的 tif说明是分幅存储后续必须做镶嵌或按需读取。# Linux/macOS 下列出 rar 内容不实际解压 unrar l 新疆维吾尔自治区地形地貌最新30m精度.rar | head -50 # 只看文件数量和总大小 unrar l 新疆维吾尔自治区地形地貌最新30m精度.rar | tail -5逻辑说明unrar l只列目录不写盘适合先判断数据规模。参数上head -50看前 50 行目录tail -5看汇总行里的文件数和总大小。如果总大小超过你剩余磁盘空间的 2 倍先别解压考虑只解压研究区涉及的图幅。2.2 用 gdalinfo 确认坐标系、分辨率和 NoData解压出 tif 后不要急着拖进 ArcGIS 或 QGIS。先用gdalinfo看一眼元数据这是避免后续投影错乱最便宜的一步。重点看四个字段Coordinate System、Pixel Size、NoData Value、Type。新疆常用投影有 WGS84 地理坐标系EPSG:4326和 CGCS2000 高斯克吕格投影如 EPSG:4544 等 3 度带。30m 分辨率在地理坐标系下表现为0.000277...度在投影坐标系下才是真正的 30 米。如果 Pixel Size 是 0.000277 左右说明是经纬度存储做面积和距离计算前必须投影。gdalinfo xinjiang_dem_30m.tif | grep -E Coordinate System|Pixel Size|NoData|Type逻辑说明grep过滤关键行避免被直方图信息刷屏。参数上如果 NoData 是-9999或-32768后续统计时要显式屏蔽如果 Type 是Int16说明是整数高程做坡度计算前建议转Float32否则梯度运算会丢精度。2.3 判断是 DEM 还是地貌分类栅格「地形地貌」四个字容易混。DEM 是连续高程像素值代表海拔地貌分类栅格是离散类别像素值代表地貌类型编码比如 1 代表山地、2 代表丘陵、3 代表平原。判断方法看像素值范围。DEM 的范围通常从负几百米到七八千米分类栅格的范围往往只有几十个整数值。也可以用gdalinfo -hist看直方图连续分布的是 DEM几个尖峰的是分类。这个判断直接决定后续分析路径DEM 做坡度坡向、起伏度分类栅格做面积统计和转移矩阵。提示如果压缩包里同时有 DEM 和地貌分类先确认两者分辨率和范围是否一致不一致时以 DEM 为基准做重采样对齐。3. 把数据读进 Python大栅格分块读取与内存控制3.1 用 rasterio 做窗口读取而不是整幅加载新疆 30m DEM 整幅加载到内存按 Float32 算面积约 166 万平方公里像素数约 18 亿单个波段就是 7 GB 以上加上坡度计算中间变量16 GB 内存的机器直接爆。正确做法是用rasterio的窗口读取每次只读一个 block 或一个研究区窗口。下面是一个按窗口计算坡度的最小示例。import rasterio from rasterio.windows import Window import numpy as np path xinjiang_dem_30m.tif with rasterio.open(path) as src: # 只读左上角 2000x2000 的窗口做测试 window Window(col_off0, row_off0, width2000, height2000) dem src.read(1, windowwindow).astype(float32) transform src.window_transform(window) nodata src.nodata # 把 NoData 替换为 nan避免参与梯度计算 if nodata is not None: dem[dem nodata] np.nan # 计算 x 和 y 方向梯度 dy, dx np.gradient(dem) slope np.degrees(np.arctan(np.sqrt(dx**2 dy**2))) print(窗口内坡度范围, np.nanmin(slope), np.nanmax(slope))逻辑说明Window指定列偏移、行偏移、宽高只把这一块读进内存。window_transform拿到该窗口的仿射变换后续写结果时要用。np.gradient默认假设像素间距为 1如果投影坐标系下像素是 30 米需要把dx和dy分别除以 30否则坡度值会偏大。参数上窗口大小根据内存调整8 GB 内存建议不超过 3000x3000。3.2 分块写入结果并保留地理元数据计算完一个窗口不能只打印要写回磁盘。用rasterio.open的w模式创建输出文件复制源文件的 CRS、transform 和 dtype然后按窗口写入。关键点是输出文件的 transform 要用窗口的 transform而不是整幅的 transform否则结果会错位。with rasterio.open(path) as src: profile src.profile.copy() profile.update(dtypefloat32, nodatanp.nan, count1) with rasterio.open(slope_test.tif, w, **profile) as dst: window Window(col_off0, row_off0, width2000, height2000) dem src.read(1, windowwindow).astype(float32) if src.nodata is not None: dem[dem src.nodata] np.nan dy, dx np.gradient(dem) slope np.degrees(np.arctan(np.sqrt(dx**2 dy**2))) dst.write(slope, 1, windowwindow)逻辑说明profile继承源文件的驱动、CRS、transform 等只改 dtype 和 nodata。dst.write的window参数确保写入位置正确。参数上如果源文件是 Int16输出坡度必须用 Float32因为坡度有小数。NoData 设为 nan 后后续用np.nanmean统计不会把无效值算进去。3.3 用 dask 或 rioxarray 做惰性计算如果不想手写窗口循环可以用rioxarray打开文件它基于 dask 做惰性分块。适合做全省统计而不是逐窗口调试。下面是一个用rioxarray计算平均高程的示例。import rioxarray as rxr dem rxr.open_rasterio(xinjiang_dem_30m.tif, chunks{x: 2000, y: 2000}) dem dem.where(dem ! dem.rio.nodata) mean_elev dem.mean().compute() print(新疆平均高程, float(mean_elev))逻辑说明chunks指定分块大小dask 会自动按块读取和计算。where屏蔽 NoData。compute()触发实际计算。参数上chunks 太小会增加调度开销太大吃内存2000x2000 是常见折中。注意rioxarray对某些压缩格式支持不如rasterio直接遇到读取报错时回退到rasterio窗口方式。4. 投影、裁剪与重采样让数据适配你的研究区4.1 地理坐标系转投影坐标系以 CGCS2000 3 度带为例如果gdalinfo显示是 EPSG:4326做坡度、面积、距离之前必须投影。新疆跨多个 3 度带常见做法是按研究区中央经线选带号。比如乌鲁木齐附近约 87°E对应 3 度带带号 29EPSG 可能为 4544 或类似。用gdalwarp一行完成。gdalwarp -t_srs EPSG:4544 -r bilinear -of GTiff \ xinjiang_dem_30m.tif xinjiang_dem_30m_proj.tif逻辑说明-t_srs指定目标投影-r bilinear是重采样方法连续高程用 bilinear分类栅格用 near。参数上如果研究区跨带不要强行用一个带号考虑分带处理或使用 Albers 等面积投影。-of GTiff指定输出格式大文件建议加-co COMPRESSLZW压缩。4.2 用研究区矢量裁剪gdalwarp 的 cutline 用法拿到投影后的数据下一步是按行政区或流域裁剪。gdalwarp支持-cutline参数直接传 shp 文件。注意 shp 的坐标系要和目标数据一致否则裁出来是空的。gdalwarp -cutline study_area.shp -crop_to_cutline \ -dstnodata -9999 -co COMPRESSLZW \ xinjiang_dem_30m_proj.tif study_area_dem.tif逻辑说明-cutline指定裁剪边界-crop_to_cutline让输出范围紧贴边界-dstnodata设置裁剪后空白区的值。参数上如果 shp 是地理坐标系而栅格是投影坐标系先ogr2ogr转 shp 或加-s_srs/-t_srs。裁剪后建议用gdalinfo -stats检查有效像素比例比例过低说明坐标系不匹配。4.3 重采样到统一分辨率30m 到 90m 的取舍有些模型要求 90m 或 1km 分辨率直接重采样会丢细节。常见做法是先用 30m 计算坡度、起伏度等派生指标再对派生指标重采样而不是对 DEM 重采样后再算。因为 DEM 重采样会平滑地形导致坡度偏小。如果必须重采样 DEM用average而不是bilinear后者在降采样时容易产生锯齿。gdalwarp -tr 90 90 -r average \ xinjiang_dem_30m_proj.tif xinjiang_dem_90m.tif逻辑说明-tr 90 90指定输出分辨率-r average做均值聚合。参数上-tr的单位跟随目标投影投影坐标系下是米地理坐标系下是度。重采样后要重新检查 NoData 是否被平均成有效值必要时在重采样前把 NoData 设为 nan 并用-srcnodata声明。注意对分类栅格重采样必须用near用average会产生 1.5 这种无意义类别值。5. 避坑与排查处理新疆30m数据时最容易翻车的 5 个点5.1 解压后文件名乱码导致脚本读不到现象在 Windows 下解压后文件名显示正常但传到 Linux 服务器或用 Python 读取时报No such file or directory。原因rar 包内文件名可能是 GBK 编码Linux 默认 UTF-8解压工具没有正确转码。解决用unrar x -p指定密码如果有或在 Windows 下用 7-Zip 解压后手动重命名更稳妥的是在 Linux 下用unrar的-sc参数尝试转换或者直接用 Python 的rarfile库读取并指定编码。5.2 坡度计算结果全是 0 或异常大现象用np.gradient算坡度结果要么全是 0要么出现 89 度这种极端值。原因一是像素间距没有代入np.gradient默认间距为 1投影坐标系下 30m 像素会导致梯度被放大 30 倍二是 NoData 没有屏蔽-9999 参与梯度计算产生巨大差值。解决把dx和dy分别除以像素宽度和高度并在计算前把 NoData 替换为 nan。如果投影是地理坐标系还要把度转米纬度方向约 111km/度经度方向乘以 cos(纬度)。5.3 裁剪后数据范围对不上裁出来是空白现象gdalwarp -cutline执行成功但输出文件全是 NoData。原因shp 和栅格坐标系不一致或者 shp 的几何类型不是多边形比如是线或点。解决先用ogrinfo -al -so study_area.shp看 shp 的 CRS 和几何类型确保是多边形且 CRS 与栅格一致。不一致时用ogr2ogr -t_srs转换 shp。另外检查 shp 是否有无效几何用ogr2ogr -makevalid修复。5.4 内存溢出整幅读取导致进程被 kill现象Python 脚本运行到src.read()时被系统 kill或者报MemoryError。原因整幅 30m 新疆 DEM 超过 7 GB加上中间变量超过物理内存。解决改用窗口读取或 dask 分块设置chunks大小不超过内存的 1/4。如果必须整幅处理考虑用gdal_translate先切分图幅或者用rasterio的block_windows逐块处理。监控内存可以用htop或 Python 的psutil。5.5 投影带号选错面积统计偏差巨大现象投影后做面积统计结果和官方公布的新疆面积差很多。原因新疆跨 73°E 到 96°E跨多个 3 度带用一个带号投影会导致边缘区域变形大。解决做全省面积统计时用 Albers 等面积投影如 EPSG:4526 或自定义中央经线 85°E而不是高斯克吕格。做局部研究时按研究区中央经线选带号。验证方法投影后算一个已知县市的面积和公开数据对比偏差超过 1% 就要检查投影。6. 进阶技巧用地貌分类栅格做地形起伏度分区统计6.1 从 DEM 计算起伏度并重分类地形起伏度是窗口内最大高程减最小高程常用 3x3 或 5x5 窗口。用scipy.ndimage的maximum_filter和minimum_filter可以快速实现。下面是一个窗口起伏度计算示例。import rasterio import numpy as np from scipy.ndimage import maximum_filter, minimum_filter with rasterio.open(study_area_dem.tif) as src: dem src.read(1).astype(float32) nodata src.nodata profile src.profile.copy() if nodata is not None: dem[dem nodata] np.nan # 用 5x5 窗口计算起伏度 max_dem maximum_filter(dem, size5) min_dem minimum_filter(dem, size5) relief max_dem - min_dem profile.update(dtypefloat32, nodatanp.nan) with rasterio.open(relief_5x5.tif, w, **profile) as dst: dst.write(relief, 1)逻辑说明maximum_filter和minimum_filter分别取窗口内最大和最小相减得到起伏度。参数上size5对应 5x5 窗口在 30m 分辨率下约 150m 范围。窗口越大起伏度越平滑但会掩盖小地形。NoData 区域在滤波后可能被邻域有效值填充建议在滤波前用np.nan屏蔽滤波后再把原 NoData 位置恢复为 nan。6.2 按地貌类型分区统计起伏度如果压缩包里附带地貌分类栅格可以用rasterio读取分类用numpy的bincount或pandas做分组统计。下面是一个按地貌类型统计平均起伏度的示例。import rasterio import numpy as np with rasterio.open(relief_5x5.tif) as src: relief src.read(1) with rasterio.open(geomorph_type.tif) as src: geom src.read(1) mask np.isfinite(relief) (geom 0) types np.unique(geom[mask]) for t in types: vals relief[mask (geom t)] print(f地貌类型 {t}: 平均起伏度 {np.mean(vals):.2f} 米, 像元数 {len(vals)})逻辑说明mask同时屏蔽起伏度无效值和地貌类型无效值。np.unique拿到所有类型循环统计。参数上如果地貌分类编码有对照表把t替换成类型名称。统计前确认两个栅格的 shape 和 transform 一致不一致先用gdalwarp对齐。6.3 用 zonal statistics 做更规范的分区统计如果研究区是行政区矢量用rasterstats做 zonal statistics 更规范。下面是一个按县统计平均高程的示例。from rasterstats import zonal_stats stats zonal_stats( counties.shp, xinjiang_dem_30m_proj.tif, stats[mean, min, max, std], nodata-9999, geojson_outFalse ) for i, s in enumerate(stats[:5]): print(f县 {i}: 平均高程 {s[mean]:.1f} 米)逻辑说明zonal_stats自动处理矢量与栅格的对齐stats指定要计算的指标nodata声明无效值。参数上矢量坐标系必须和栅格一致否则结果为空。大矢量建议先按研究区裁剪栅格减少计算量。rasterstats底层用rasterio和shapely对大数据量可能较慢可以分县循环。6.4 我踩过的坑和现在的习惯最早处理这类数据时我习惯先解压到桌面再拖进 QGIS结果 30m 新疆 DEM 让 QGIS 卡死好几次后来改成先用gdalinfo看元数据再用gdalwarp裁出研究区最后才进 QGIS 出图。另一个血泪经验是投影带号曾经用 3 度带 28 带做全省统计面积比实际少了 8%后来改用 Albers 才对齐。现在我的固定流程是unrar l看结构 →gdalinfo看元数据 →gdalwarp投影和裁剪 → Python 窗口读取做分析 → 结果用rasterio写回并检查 NoData。这套流程不新鲜但能避开九成以上的翻车。希望帮到你。本文还有配套的精品资源点击获取