简介面向GISer与地理科研工作者的广东省地形地貌分类数据集采用30米精度按海拔、起伏程度、陆地地貌类型及中国地貌类型等标准划分可直接用于区域地貌制图、生态环境分析及国土空间规划。压缩包共30个文件主要以tif栅格数据为核心配套dbf属性表、tfw坐标参考、xml元数据等辅助文件还有2个rar及示意图可辅助理解整体仅9.51MB轻量易下载。已有200人学习、下载实用性获得一定认可。数据内细分低/中/高海拔、丘陵/小起伏/大起伏等地貌形态同时涵盖冲积、洪积、海积等成因类型配合WGS84与Albers双坐标系便于直接接入ArcGIS、QGIS等平台完成专题图输出。1. 一份30m精度地形栅格为什么还要拆开看拿到“广东省地形地貌最新30m精度.rar”第一反应通常是解压找DEM然后直接算坡度坡向。但解压之后你会看到文件列表里根本没有传统的DEM高程而是五套TIF综合地貌、海拔分类、陆地地貌类型、中国地貌类型、起伏程度分类。也就是说这包资源已经把“地貌形态”做成了分类结果不是让你从高程开始推演的原始数据。对做GIS开发、工程选址、环境影响评价的工程师来说30m精度意味着一个像元对应地面30米乘30米能够支撑县域以上级别的资源环境分析。但如果你直接打开landform_广东省.tif会发现问题很多成因类型、形态类型、海拔类型混在一个图例里直接做专题图并不好用。真正的价值在于那三张独立的分类栅格——海拔、起伏度、地貌成因它们才是能直接统计、重分类、接入建模的干净数据。这篇博文就围绕这包RAR资源从文件结构、坐标系统、Python读取、分类应用到排错把它变成你工程里可复用的GeoData。2. 解包后的文件清单TIF、TFW与属性表的含义2.1 五类核心栅格文件的分工这包RAR里的文件名非常直白但第一次接触的人很容易被重复的文件名后缀搞晕。比如landform_广东省.tif是综合地貌分类的合成结果而海拔分类、陆地地貌类型、中国地貌类型、起伏程度分类则分别对应不同分类体系。下表列出核心文件与用途文件类型实际内容landform_广东省.tif综合分类栅格融合海拔、起伏、成因的最终分类海拔分类_广东省.tif分类栅格低海拔、中海拔、中高海拔、高海拔、极高海拔起伏程度分类_广东省.tif分类栅格丘陵、小起伏、中起伏、大起伏、极大起伏陆地地貌类型_广东省.tif分类栅格山脉、丘陵、平原、沟壑等形态类型中国地貌类型_广东省.tif分类栅格按成因分为海积、湖积、冲积、洪积、风积、冰碛、剥蚀等各.tfw世界文件栅格到地理坐标的仿射变换参数各.aux.xmlGDAL元数据坐标系、NoData、统计信息各.vat.dbf属性表栅格值到中文类别名的映射值得注意陆地地貌类型和中国地貌类型是两个维度。前者回答“地形长什么样”比如山脉、丘陵后者回答“是怎么形成的”比如冲积、海积、冰碛。广东以珠江三角洲冲积平原和丘陵山地为主但字段保留的是全国统一分类标准所以海积、湖积、冰碛这些类别都会出现。以前我把陆地地貌类型当成因直接查表结果统计出来的“冰碛”面积在粤北居然不是零后来才发现查错了文件。另一个常见误区是这组TIF没有高程波段只有分类编码。你要算坡度、坡向必须另外找ASTER GDEM或SRTM30m分辨率对应的多为30米目。这包资源能做的是直接告诉你这个格子是低海拔丘陵那个格子是冲积平原。2.2 TFW与aux.xml栅格地理定位的钥匙TFW是ESRI的Tiff World File六行纯文本。很多人以为把它删了也能凑合用但实际上GIS无外扩文件时只有像素坐标没有空间坐标。六行数值分别对应第1行像元在x方向的实际宽度第2行x方向旋转参数正北对齐时为0第3行y方向旋转参数第4行像元在y方向的实际高度通常为负数第5行左上角像元中心的x坐标第6行左上角像元中心的y坐标比如Albers投影下常见的30和-30加上左上角坐标就能恢复完整的仿射变换。但TFW不存坐标系坐标系信息在aux.xml里。用QGIS打开某张tif如果显示“未知投影”大概率是aux.xml被单独剔除了。所以解压时一定要把所有附属文件放在同一目录不要只拷贝tif否则后续坐标转换会吃大亏。2.3 vat.dbf把像素值翻译成中文类别名分类栅格的值是整数编码必须对应vat.dbf才能知道1到底代表低海拔还是中海拔。读取DBF最简单的方式是dbfreadfrom dbfread import DBF table DBF(海拔分类_广东省.tif.vat.dbf) for record in table: print(record[Value], record[Name])这段代码遍历属性表输出像素值与类别名的对应关系。参数说明Value是栅格内的像素值Name是分类名称。在制作专题图前务必先执行这一段确认类别顺序。因为不同生产单位对低海拔、中高海拔的编码顺序可能不同直接拿数值去做重分类错一位后面所有统计就全错。3. 用PythonGDAL解析30m精度栅格元数据、数值与可视化3.1 先跑一遍gdalinfo坐标系和像元尺寸一目了然拿到TIF第一件事不是直接读数组而是跑一遍gdalinfo看元数据gdalinfo -stats 海拔分类_广东省.tif输出里重点看四块Size行列数、Origin左上角坐标、Pixel Size像元尺寸、Coordinate System坐标系。这包资源同时给出WGS84和Albers_Conic_Equal_Area两套坐标所以Pixel Size有两种可能一个接近0.00027经纬度单位一个接近30米。如果你发现同一个文件在ArcGIS和GDAL里看到的分辨率不同多半是投影单位没注意。在Python里读取同样信息from osgeo import gdal ds gdal.Open(海拔分类_广东省.tif) print(尺寸:, ds.RasterXSize, ds.RasterYSize) print(仿射变换:, ds.GetGeoTransform()) print(投影:, ds.GetProjection()) band ds.GetRasterBand(1) print(数据类型:, gdal.GetDataTypeName(band.DataType)) print(NoData:, band.GetNoDataValue()) print(统计:, band.ComputeStatistics(True))这里GetGeoTransform()返回六个浮点数顺序与TFW六行对应GetProjection()返回WKT格式的投影定义ComputeStatistics(True)会强制重新计算统计值。注意如果GetNoDataValue()返回None说明TIF没有显式设置NoData但对海洋区域可能用0填充需要在后续处理中手动掩膜。3.2 让分类栅格真正“看得见”离散色带渲染分类栅格是离散整数值不能用连续拉伸显示否则低海拔和中海拔的颜色过渡会让人误以为存在渐变。常见做法是自定义离散色带用matplotlib把数组画出来import numpy as np from osgeo import gdal import matplotlib.pyplot as plt arr gdal.Open(海拔分类_广东省.tif).ReadAsArray() # 假设0无数据1低海拔2中海拔3中高海拔4高海拔5极高海拔 colors [#ffffff, #a1d99b, #74c476, #41ab5d, #238b45, #005a32] cmap plt.matplotlib.colors.ListedColormap(colors, elevation) bounds [-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5] norm plt.matplotlib.colors.BoundaryNorm(bounds, cmap.N) plt.figure(figsize(10, 12)) plt.imshow(arr, cmapcmap, normnorm) plt.colorbar(ticks[0, 1, 2, 3, 4, 5], label海拔分类编码) plt.title(广东省海拔分类 (30m)) plt.axis(off) plt.savefig(海拔分类_广东.png, dpi200)这里的BoundaryNorm按类别边界映射颜色ticks显示的是像素值对应的图例。如果某个类别没有分布图例上仍然会显示这是正常现象。广东最高峰石坑崆约1902米所以高海拔和极高海拔两个类别基本不会出现但保留它们可以让整套数据和全国其他省份拼接时图例保持一致。3.3 转Albers_Conic_Equal_Area后再做面积统计WGS84坐标下的栅格像元面积不等大纬度越高单个像元代表的实际面积越小。要统计各类地貌面积必须转换到等积投影。Albers_Conic_Equal_Area是中国中纬度区域的常用选择。命令行转换gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 \ -tr 30 30 -r near \ 海拔分类_广东省.tif 海拔分类_广东省_aea.tif参数说明-t_srs设置Albers双标准纬线和中央经线-tr 30 30强制像元重采样为30米-r near使用最邻近法插值。最邻近法保证重采样后的像素值仍然是原始分类编码如果用bilinear会出现1.2这种无效值。这一点对分类栅格是致命的但对高程栅格反而需要bilinear两者不能混用。Python中对应调用的是gdal.Warp参数名与命令行几乎一致。转换完成后再统计面积时直接使用像元数乘以30*30即可精度足够工程分析。4. 海拔与起伏度分类的实战重分类、叠加与区域统计4.1 广东地形特征与分类阈值参考广东地势北高南低最高峰在粤北南岭地区海拔接近1902米因此全省以低海拔和中海拔为主。数据中的分类字段沿用了全国统一标准但阈值在不同生产单位之间会有差异。常见参考如下海拔分类参考范围米起伏度分类参考范围米低海拔 500丘陵 200中海拔500 ~ 1000小起伏200 ~ 500中高海拔1000 ~ 2000中起伏500 ~ 1000高海拔2000 ~ 4000大起伏1000 ~ 2500极高海拔 4000极大起伏 2500这份阈值表可以用于快速理解数据但不能作为调整分类的依据。真正要改阈值应该打开资源自带的vat.dbf看原始范围说明。之前接过一个项目对方把“中高海拔”误读为海拔3000米以上导致粤北地区大面积被划入“低海拔”整个生态红线评估全废。4.2 用NumPy叠加海拔与起伏度生成“适建区”掩模工程上很少直接使用五级分类更常见的是把若干类别合并成二值掩膜。比如适宜建设区域通常要求低海拔、起伏度小。用NumPy叠加两张栅格import numpy as np from osgeo import gdal # 读取海拔分类 ds_hgt gdal.Open(海拔分类_广东省.tif) hgt ds_hgt.GetRasterBand(1).ReadAsArray().astype(np.float32) # 读取起伏度分类 ds_und gdal.Open(起伏程度分类_广东省.tif) und ds_und.GetRasterBand(1).ReadAsArray().astype(np.float32) # 假设海拔编码1低海拔起伏度编码1丘陵 # 需要先打印vat确认这里以常见编码为例 suit np.where((hgt 1) (und 1), 1, 0) # 写为GeoTIFF driver gdal.GetDriverByName(GTiff) out_ds driver.Create(适宜建设区_test.tif, ds_hgt.RasterXSize, ds_hgt.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds_hgt.GetGeoTransform()) out_ds.SetProjection(ds_hgt.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(suit) out_band.SetNoDataValue(0) out_ds.FlushCache()代码逻辑分三步先读两张分类图再用np.where把满足条件的像素赋值为1最后写为GeoTIFF。关键参数是hgt 1和und 1中的编码值必须根据vat.dbf确认。如果两张栅格的行列数不一致比如一张被重采样过在叠加前需要对齐边界。最简单的对齐方法是用gdal.Warp的-cutline或-te参数统一范围避免数组维度错位。4.3 用zonal_stats做行政区尺度的分类统计掩膜之后就要统计面积统计到地级市或县域。比较顺手的工具是rasterstats配合GeoPandas可以直接得到每个行政区的分类像元数import geopandas as gpd from rasterstats import zonal_stats gdf gpd.read_file(广东省地级市.shp) stats zonal_stats( gdf, 海拔分类_广东省_aea.tif, band_num1, categoricalTrue, geojson_outTrue ) result gpd.GeoDataFrame.from_features(stats) result[低海拔面积_km2] result[低海拔] * 30 * 30 / 1e6参数说明categoricalTrue会返回一个字典键是像素值值是像元数量后面把像元数乘以30米乘以30米再除以1e6转成平方千米。这里有一个前提输入的栅格必须已经是Albers等积投影否则直接按面积计算会偏差很大。geojson_outTrue让返回结果保留原矢量属性方便后续join回边界数据。如果内存吃紧可以用rasterio搭配numpy.bincount做更轻量的统计但rasterstats在矢量多边形数量不多时足够用。广东21个地级市、几百个乡镇几分钟就能跑完。5. 30m栅格排错NoData、坐标系混淆与批量切片5.1 先检查proj4字符串别被显示坐标骗了WGS84和Albers_Conic_Equal_Area混在同一包资源里时最典型的错误是把经纬度当米去算面积。打开GDAL的GetProjection()返回的WKT找到PROJCS和PROJCRS关键字。如果PROJCS是Albers单位是米如果CORRDERED_ATTR是[2]且没有PROJCS那是经纬度。一旦确认源坐标再决定是否做投影转换。不要只看地图软件底部的坐标值因为QGIS可能会动态重投影显示的是经纬度但并不代表原始数据是WGS84。5.2 NoData与0值统计前必须掩膜分类栅格在沿海区域可能出现0值或NoData。如果统计时不过滤会被算作未知类别污染面积。用GDAL命令行快速生成掩膜gdal_calc.py -A 海拔分类_广东省.tif --outfilemask.tif --calcA0 --NoDataValue0这段命令把栅格中大于0的像元设为1其余设为0并显式指定NoDataValue0。后续做面积统计时先乘上mask再统计非零值就能排除海洋与无效区域。注意--NoDataValue0与calc里的0需要保持一致否则会出现“有分类但显示透明”的诡异结果。5.3 批量输出TMS瓦片供前端地图加载30m全省数据压缩后仍可能有几百MB直接扔到前端地图不现实。常见做法是先重投影到Web Mercator再用gdal2tiles切瓦片gdal2tiles.py -z 8-14 -s EPSG:4326 -w none 海拔分类_广东省_aea.tif tiles/-z 8-14表示只切8级到14级瓦片-w none去掉网页空白说明文件输出目录是tiles/。这里源数据用Albers也可以gdal2tiles会自动重采样到3857但建议先用gdalwarp转一次避免切图到一半内存爆掉。切完后用Nginx直接指向tiles目录前端Leaflet通过L.tileLayer(tiles/{z}/{x}/{y}.png)加载即可。30m精度在14级时每个瓦片看起来依然清晰适合做省级或县域范围内的快速浏览。本文还有配套的精品资源点击获取