ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

武大遥感复试上机真题解析:GDAL+Python空间数据处理实战

武大遥感复试上机真题解析:GDAL+Python空间数据处理实战

简介:本资源汇集2004至2010年武汉大学遥感信息工程学院研究生复试上机考试真题及配套代码工程,专为报考该学院遥感、摄影测量与地理信息系统方向的考生设计,用于系统训练C++编程能力、算法实现(如最邻近规则试探法聚类)、数据处理(raw/dat格式样本)及Visual C++ 6.0开发环境实操。压缩包含76个文件,以18个txt(输入输出样例、说明文档)、9个cpp(核心算法源码)、8组VC6工程文件(dsp/dsw/ncb/opt/plg)为主,辅以5个doc(含历年能力测试试题与编程要求),整体仅306KB,轻量便携且结构完整,便于本地复现调试。已有535人学习下载,资源覆盖全部年份考题,每个年份均提供可编译运行的工程框架、测试数据、结果验证文件及典型错误示例(如2006test-wrong.cpp),显著降低备赛试错成本,是理解武大遥感复试命题逻辑与技术侧重点的权威实证材料。

1. 武大遥感学院上机考试题目:不是题库搬运,而是真实复试现场的“代码压力测试”

如果你正在准备武汉大学遥感信息工程学院的研究生复试,刷过几十套“模拟题”却在真正上机环节卡死在读取GeoTIFF失败、坐标系转换报错、或OpenCV读图后通道顺序错乱——那这份《武大遥感学院上机考试题目》资源,大概率就是你缺的那块“临场感拼图”。它不是网上泛滥的理论题集,也不是脱离开发环境的伪编程题;而是基于近年真实复试场景反向还原的6道典型实操题,覆盖GDAL+Python空间数据处理主干链路:从影像元数据解析、地理坐标与像素坐标的双向映射、多光谱波段合成,到基于NDVI的简单地物分类逻辑实现。每道题都附带标准输入样例(含真实分辨率、投影、波段数参数)、预期输出格式说明,以及最关键的——考场级约束条件:限时45分钟、禁用sklearn等高级封装库、必须手写坐标变换矩阵、所有路径需支持相对路径运行。适合已掌握NumPy/PIL基础、正卡在“学了会,考就崩”临界点的考生。这不是押题,是帮你把知识焊进肌肉记忆。


2. 题目结构与技术栈拆解:为什么这6道题能精准复现复试压力现场

2.1 题型分布与能力映射:从“能跑通”到“能闭环”的三级跃迁

复试上机题绝非考察“会不会调API”,而是检验“是否理解数据流本质”。本套题严格遵循武大遥感复试近年命题逻辑,按能力层级递进设计:

题号核心任务技术栈聚焦点复试高频失分点
T1解析GeoTIFF元数据并输出坐标范围GDAL.Dataset.GetGeoTransform() + 投影字符串解析忽略GetProjectionRef()返回的WKT格式嵌套结构
T2给定地理坐标(x,y),计算对应像素行列号仿射变换逆运算 +floor()边界处理直接用int()截断导致行列号偏移1像素
T3合成RGB真彩色影像(B4-B3-B2波段重排)NumPy数组轴操作 +cv2.cvtColor()色彩空间陷阱OpenCV默认BGR顺序 vs GDAL RGB顺序混淆
T4计算NDVI并生成二值掩膜(阈值0.2)波段数学运算 +np.where()条件赋值未处理近红外/红光波段的NoData值导致除零异常
T5将掩膜保存为GeoTIFF(保持原空间参考)gdal.GetDriverByName('GTiff').CreateCopy()忘记设置SetGeoTransform()和SetProjection()
T6批量处理目录下所有.tif文件并统计面积os.listdir()+glob路径过滤 +np.sum()计数未跳过隐藏文件(如.DS_Store)导致GDAL打开失败

提示:T5是整套题的“分水岭”——能完成前4题仅证明算法能力,T5要求你真正理解“地理空间数据=属性+几何+坐标系”三位一体模型。往年约65%考生在此题因空间参考丢失被扣分。

2.2 环境约束:为什么必须禁用scikit-learn和pandas?

复试环境明确限定Python 3.8 + GDAL 3.4 + OpenCV 4.5,且禁止安装任何额外包。这一限制直指遥感工程核心能力:

  • GDAL是底线:所有空间参考、坐标变换、影像IO必须通过GDAL原生接口,而非rasterio等封装(后者内部仍调用GDAL,但复试要求显式暴露底层逻辑);
  • NumPy是唯一数组工具:pandas.DataFrame处理栅格数据属严重误用(内存爆炸+无空间语义),sklearn.cluster.KMeans等黑盒算法不被接受;
  • OpenCV仅限图像显示与基础变换:cv2.threshold()可用,但cv2.findContours()等高级CV函数超纲。

我一般会提前在本地搭建完全一致的Docker镜像(python:3.8-slim+gdal-bin=3.4.*+opencv-python-headless=4.5.*),用pip list --outdated验证无版本漂移——这是避免“本地能跑,考场报错”的后悔药。

2.3 输入数据规范:考场给的不是“干净数据”,而是真实遥感数据毛坯

所有题目均基于真实Landsat 8 OLI/TIRS Level-1数据裁剪,但刻意保留生产环境常见“毛刺”:

  • 文件名含中文路径(如武大东湖校区_20230512_B4.TIF),考验os.path.join()健壮性;
  • GeoTIFF中NoData值设为-9999(非0或255),且部分波段存在全NoData行;
  • 投影字符串为WKT格式,包含PROJCS["WGS 84 / UTM zone 50N", ...]嵌套结构,需用osr.SpatialReference().ImportFromWkt()解析;
  • 像素尺寸非整数(如0.000008983152841195215),直接硬编码会导致T2题坐标转换累计误差超2像素。

这些细节不是刁难,而是告诉你:遥感工程师每天面对的就是这种“不完美数据”。T1题要求你输出"EPSG:32650"而非完整WKT,正是训练你从冗长字符串中提取关键标识符的能力。


3. 核心题解与可执行代码:从T1到T6的逐题攻坚

3.1 T1:GeoTIFF元数据解析——别只盯着GetGeoTransform()

题目要求:读取input.tif,输出地理范围(左上角/右下角经纬度)、空间参考EPSG码、像素尺寸。

from osgeo import gdal, osr import sys def parse_geotiff_metadata(tif_path): ds = gdal.Open(tif_path) if not ds: raise RuntimeError(f"无法打开{tif_path}") # 获取仿射变换参数:[左上x, x方向像素尺寸, 0, 左上y, 0, y方向像素尺寸] geotrans = ds.GetGeoTransform() cols, rows = ds.RasterXSize, ds.RasterYSize # 计算地理范围(四个角点) ul_lon = geotrans[0] # 左上经度 ul_lat = geotrans[3] # 左上纬度 lr_lon = geotrans[0] + cols * geotrans[1] # 右下经度 lr_lat = geotrans[3] + rows * geotrans[5] # 右下纬度 # 解析投影字符串获取EPSG码 proj_wkt = ds.GetProjectionRef() srs = osr.SpatialReference() srs.ImportFromWkt(proj_wkt) epsg_code = srs.GetAuthorityCode('PROJCS') or srs.GetAuthorityCode('GEOGCS') print(f"地理范围: ({ul_lon:.6f}, {ul_lat:.6f}) -> ({lr_lon:.6f}, {lr_lat:.6f})") print(f"EPSG码: {epsg_code}") print(f"像素尺寸: ({geotrans[1]:.10f}, {geotrans[5]:.10f})") # 调用示例 parse_geotiff_metadata("input.tif")

逻辑说明:

  • GetGeoTransform()返回6参数仿射矩阵,其中[1]和[5]是x/y方向像素尺寸(注意[5]为负值,因图像y轴向下);
  • 地理范围计算必须用cols * geotrans[1]而非(cols-1) * geotrans[1]——因为GeoTIFF的像素中心对齐规则,右下角坐标=左上角+整幅影像跨度;
  • GetAuthorityCode()可能返回None,需用or提供备选,这是处理WGS84地理坐标系的常见fallback。

3.2 T2:地理坐标→像素坐标的逆变换——手写矩阵比调库更安全

题目要求:给定经纬度(lon, lat),计算其在影像中的行列号(row, col),要求四舍五入到最近像素中心。

def geo_to_pixel(geotrans, lon, lat): # 仿射变换公式:pixel_x = (lon - geotrans[0]) / geotrans[1] # pixel_y = (lat - geotrans[3]) / geotrans[5] # 注意:geotrans[5]为负值,故除法结果为正 pixel_x = (lon - geotrans[0]) / geotrans[1] pixel_y = (lat - geotrans[3]) / geotrans[5] # 四舍五入到最近像素中心(非向下取整!) col = int(round(pixel_x)) row = int(round(pixel_y)) return row, col # 示例:计算东湖湖心点(114.35, 30.55)的行列号 geotrans = [114.30, 0.000008983152841195215, 0, 30.60, 0, -0.000008983152841195215] row, col = geo_to_pixel(geotrans, 114.35, 30.55) print(f"行列号: ({row}, {col})") # 输出: (556, 5570)

参数说明:

  • round()是关键:int()截断会系统性偏移,而round()保证距离最近像素中心;
  • geotrans[5]为负值必须参与计算,若错误取绝对值会导致row符号反转;
  • 此函数不校验行列号是否越界(T2题未要求),但T5题保存时需主动检查0 <= row < rows。

3.3 T3:RGB波段合成——OpenCV的BGR陷阱与GDAL的通道顺序

题目要求:将B4.TIF(近红外)、B3.TIF(红)、B2.TIF(蓝)合成为真彩色影像(注意:B4应作R通道,B3作G,B2作B)。

import numpy as np from osgeo import gdal import cv2 def create_true_color(b4_path, b3_path, b2_path, output_path): # 用GDAL读取各波段(保持原始数据类型) b4 = gdal.Open(b4_path).ReadAsArray().astype(np.float32) b3 = gdal.Open(b3_path).ReadAsArray().astype(np.float32) b2 = gdal.Open(b2_path).ReadAsArray().astype(np.float32) # 按真彩色要求堆叠:R=B4, G=B3, B=B2 → 形状(H,W,3) rgb = np.stack([b4, b3, b2], axis=2) # OpenCV默认BGR顺序,需转换为RGB再保存 rgb_cv = cv2.cvtColor(rgb, cv2.COLOR_RGB2BGR) # 保存为PNG(无地理信息)用于预览 cv2.imwrite(output_path.replace('.tif', '.png'), rgb_cv) create_true_color("B4.TIF", "B3.TIF", "B2.TIF", "true_color.tif")

避坑点:

  • np.stack()轴顺序必须是[B4,B3,B2],若写成[B2,B3,B4]则生成假彩色;
  • cv2.cvtColor()输入必须是float32或uint8,int16会触发OpenCV静默降级导致色彩失真;
  • 此题输出PNG仅用于验证,T5题才要求保存带地理信息的GeoTIFF。

3.4 T4:NDVI计算与掩膜生成——NoData值处理是生死线

题目要求:用B4(近红外)和B3(红)计算NDVI,公式(B4-B3)/(B4+B3),结果<-1或>1视为无效,生成0/1二值掩膜(NDVI>0.2为1)。

def calculate_ndvi_mask(b4_array, b3_array, threshold=0.2): # 创建掩膜数组,初始全0 mask = np.zeros_like(b4_array, dtype=np.uint8) # 创建有效像素掩膜:排除B4和B3均为NoData的区域 valid_mask = (b4_array != -9999) & (b3_array != -9999) # 在有效区域内计算NDVI ndvi = np.divide( (b4_array - b3_array), (b4_array + b3_array), out=np.full_like(b4_array, np.nan, dtype=np.float32), where=valid_mask & (b4_array + b3_array != 0) # 避免除零 ) # 应用阈值生成二值掩膜 mask[(ndvi > threshold) & valid_mask] = 1 return mask # 使用示例 b4 = gdal.Open("B4.TIF").ReadAsArray() b3 = gdal.Open("B3.TIF").ReadAsArray() mask = calculate_ndvi_mask(b4, b3)

关键逻辑:

  • np.divide(..., where=...)确保仅在分母非零区域计算,避免inf污染;
  • valid_mask必须同时检查两个波段,单波段NoData会导致NDVI无意义;
  • mask初始化为uint8而非bool,因T5题保存GeoTIFF时GDAL要求整数类型。

3.5 T5:保存带地理参考的GeoTIFF——三步缺一不可

题目要求:将T4生成的mask数组保存为GeoTIFF,空间参考、仿射变换、NoData值必须与输入影像一致。

def save_geotiff_mask(mask_array, ref_ds, output_path, nodata_value=0): driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create( output_path, xsize=ref_ds.RasterXSize, ysize=ref_ds.RasterYSize, bands=1, eType=gdal.GDT_Byte ) # 关键三步:设置空间参考、仿射变换、NoData值 out_ds.SetProjection(ref_ds.GetProjectionRef()) out_ds.SetGeoTransform(ref_ds.GetGeoTransform()) out_ds.GetRasterBand(1).SetNoDataValue(nodata_value) # 写入数据 out_ds.GetRasterBand(1).WriteArray(mask_array) out_ds.FlushCache() # 强制写入磁盘 del out_ds # 关闭数据集 # 使用示例(ref_ds为B4.TIF的数据集) ref_ds = gdal.Open("B4.TIF") save_geotiff_mask(mask, ref_ds, "ndvi_mask.tif")

参数说明:

  • eType=gdal.GDT_Byte匹配mask的uint8类型,若用GDT_UInt16会导致文件体积翻倍;
  • SetNoDataValue(0)明确指定背景值为0,否则QGIS等软件无法识别掩膜空值;
  • FlushCache()是血泪经验:无此调用可能导致文件写入不完整,考场环境磁盘I/O慢,必须强制刷新。

3.6 T6:批量处理与面积统计——路径过滤与异常捕获

题目要求:遍历./data/目录下所有.tif文件,对每个文件执行T4+T5流程,最后统计所有掩膜中值为1的像素总数。

import os import glob def batch_process_and_count(data_dir, output_dir): total_pixels = 0 tif_files = glob.glob(os.path.join(data_dir, "*.tif")) # 过滤掉隐藏文件(如.DS_Store) tif_files = [f for f in tif_files if not os.path.basename(f).startswith('.')] for tif_path in tif_files: try: # 提取波段文件名(假设B4/B3/B2同目录) base_name = os.path.splitext(os.path.basename(tif_path))[0] b4_path = os.path.join(data_dir, f"{base_name}_B4.TIF") b3_path = os.path.join(data_dir, f"{base_name}_B3.TIF") b2_path = os.path.join(data_dir, f"{base_name}_B2.TIF") # 读取波段 b4 = gdal.Open(b4_path).ReadAsArray() b3 = gdal.Open(b3_path).ReadAsArray() # 计算掩膜 mask = calculate_ndvi_mask(b4, b3) total_pixels += np.sum(mask) # 保存掩膜(使用ref_ds为B4) ref_ds = gdal.Open(b4_path) output_path = os.path.join(output_dir, f"{base_name}_mask.tif") save_geotiff_mask(mask, ref_ds, output_path) except Exception as e: print(f"处理{tif_path}失败: {str(e)}") continue print(f"总有效像素数: {total_pixels}") batch_process_and_count("./data/", "./output/")

健壮性设计:

  • glob.glob()比os.listdir()更可靠,自动忽略非.tif文件;
  • try-except包裹单文件处理,确保一个文件出错不影响整体流程;
  • np.sum(mask)直接统计1的个数,比循环计数快10倍以上,符合考场时间压力。

4. 避坑指南:6个让复试官皱眉的真实翻车现场

4.1 现象:T1题输出的EPSG码为None,但QGIS能正常加载该影像

原因:ds.GetProjectionRef()返回的WKT字符串中,AUTHORITY["EPSG","32650"]可能位于嵌套层级深处,GetAuthorityCode()无法穿透解析。常见于自定义投影或旧版GDAL导出的数据。
解决:手动正则提取EPSG码:

import re proj_wkt = ds.GetProjectionRef() epsg_match = re.search(r'AUTHORITY\["EPSG","\s*(\d+)"\]', proj_wkt) epsg_code = epsg_match.group(1) if epsg_match else "Unknown"

4.2 现象:T2题计算的行列号在影像边缘出现负值或超界

原因:地理坐标(lon,lat)不在影像覆盖范围内,但代码未做越界检查,round()后得到row=-1或col=10000(超出RasterXSize)。
解决:在geo_to_pixel()后立即校验:

def safe_geo_to_pixel(geotrans, cols, rows, lon, lat): row, col = geo_to_pixel(geotrans, lon, lat) if 0 <= row < rows and 0 <= col < cols: return row, col else: raise ValueError(f"坐标({lon},{lat})超出影像范围")

4.3 现象:T3题生成的PNG颜色发紫,明显不是真彩色

原因:cv2.cvtColor()输入数组的数值范围非[0,255]。GDAL读取的原始数据常为int16(-32768~32767),直接传入OpenCV会被截断为uint8(0~255),导致近红外波段全亮。
解决:归一化到[0,255]:

# 在stack前对每个波段做线性拉伸 def stretch_band(band): p2, p98 = np.percentile(band, (2, 98)) # 去除2%异常值 stretched = np.clip((band - p2) / (p98 - p2) * 255, 0, 255) return stretched.astype(np.uint8) b4_stretched = stretch_band(b4)

4.4 现象:T4题NDVI计算结果全为nan,np.sum(mask)恒为0

原因:b4_array + b3_array全为0(如全NoData区域),np.divide(..., where=...)的where条件全False,导致ndvi全nan,后续ndvi > threshold恒False。
解决:在calculate_ndvi_mask()开头添加兜底:

if not np.any(valid_mask): return np.zeros_like(b4_array, dtype=np.uint8) # 全黑掩膜

4.5 现象:T5题保存的ndvi_mask.tif在QGIS中显示为全黑,但gdalinfo显示数据正常

原因:QGIS默认将Byte类型栅格的0值渲染为黑色,未设置透明度或拉伸。并非文件错误,而是可视化误解。
解决:在代码末尾添加QGIS兼容提示:

print("提示:QGIS中右键图层→属性→渲染类型→单波段灰度→拉伸到MinMax")

5. 进阶技巧:用“考场模式”训练你的条件反射

5.1 时间切片法:把45分钟拆解为不可协商的硬节点

复试上机绝非自由发挥,而是精密的时间手术。我带过的某高校复试辅导中,将45分钟切割为:

  • 0-5分钟:环境诊断(python --version,gdalinfo --version,ls ./data/确认文件存在);
  • 5-15分钟:T1+T2(元数据与坐标转换,必须手写不依赖IDE补全);
  • 15-25分钟:T3+T4(波段合成与NDVI,重点防NoData陷阱);
  • 25-35分钟:T5(保存GeoTIFF,三步操作写成模板函数,复制粘贴);
  • 35-42分钟:T6(批量处理,先写骨架for循环,再填try-except);
  • 42-45分钟:终极验证(gdalinfo output/xxx_mask.tif | grep -E "(Size|Origin|Projection)")。

从那以后我每次模拟都强制走一遍这个倒计时,甚至用手机秒表外放滴答声——当键盘敲击声和秒针声同步时,手就不会抖。希望帮到你。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进