请问landsat8-OLI数据怎么进行正射校正？

Question

数据已经经过辐射定标和大气校正了，但不知道怎么进行正射校正，请问具体步骤和相关参数是什么啊？

Answer 1

引自免费微信小程序：皆我百晓生

正射校正是将卫星图像从地球表面的投影转换为平面坐标系统的过程，这样可以确保图像中的每个像素都对应地球上一个确切的位置。Landsat 8 OLI数据通常提供的是地理编码的产品，即已经包含了某种程度的正射校正，但可能不够精确，特别是对于高精度应用而言。如果你需要更高精度的正射校正，可以使用像ENVI、QGIS或ERDAS Imagine这样的专业遥感软件，或者是开源库如GDAL（Geospatial Data Abstraction Library）。这里我将为你提供使用GDAL的基本命令行操作步骤。

1. 安装GDAL：首先确保你已经安装了GDAL。在Linux上，可以通过包管理器安装：

sudo apt-get install gdal-bin

或者在Windows上，可以从官方网站下载并安装GDAL。

1. 准备数据：确保你有Landsat 8 OLI的数据，包括主文件（例如：LC08_L1TP_047036_20220519_20220525_01_T1_B1.TIF等）、XML文件（例如：LC08_L1TP_047036_20220519_20220525_01_T1_MTL.txt）以及DEM数据（例如SRTM或其他高程数据）。

2. 使用GDAL warp命令进行正射校正：在命令行中执行以下命令，其中input_file是你的Landsat 8 OLI数据文件名，output_file是你想要保存结果的文件名，dem_file是数字高程模型文件，utm_zone是UTM投影带号，可以根据你的地理位置确定，proj_epsg是输出投影的EPSG代码，这里我们选择常用的WGS84 UTM投影（EPSG:32632），pixel_size是你想要的输出像素大小（单位米）。

gdalwarp -t_srs "+proj=utm +zone=${utm_zone},+datum=WGS84" -tr ${pixel_size} ${pixel_size} -r near -co "TILED=YES" -co "COMPRESS=LZW" -dstnodata -9999 input_file output_file -rpc -srcnodata -9999 -multi -overwrite dem_file

注意替换上述命令中的变量值，比如输入文件名、输出文件名、DEM文件名、UTM区号等。另外，如果你的原始数据已经是UTM投影，那么t_srs参数可以去掉。

这个过程可能会产生一些警告，但通常不会影响结果。如果你对结果不满意，可以尝试调整内插方法（如-r选项）或使用其他工具进行进一步优化。

如果你需要在Python脚本中完成这个任务，可以使用GDAL的Python API。这是一个简单的例子：

from osgeo import gdal

# 设置输入输出参数
input_file = 'LC08_L1TP_047036_20220519_20220525_01_T1_B1.TIF'
output_file = 'output.tif'
dem_file = 'srtm_n37_e008.hgt'  # 假设你的SRTM DEM文件名为srtm_n37_e008.hgt
utm_zone = 32
pixel_size = 30
epsg_code = 32632

# 创建输出驱动
driver = gdal.GetDriverByName('GTiff')

# 打开输入文件
in_dataset = gdal.Open(input_file)

# 获取输入文件的元数据
geotransform = in_dataset.GetGeoTransform()
projection = in_dataset.GetProjection()

# 计算新地理变换
new_gt = [geotransform[0], pixel_size, geotransform[2],
          geotransform[3], geotransform[5] - pixel_size, pixel_size]

# 创建输出数据集
out_dataset = driver.Create(output_file, in_dataset.RasterXSize, in_dataset.RasterYSize, in_dataset.RasterCount,
                            in_dataset.GetRasterBand(1).DataType)
out_dataset.SetGeoTransform(new_gt)
out_dataset.SetProjection(projection)

# 进行正射校正
gdal.Warp(output_file, input_file, options='-t_srs EPSG:%d -tr %f %f -r bilinear -co TILED=YES -co COMPRESS=LZW -dstnodata -9999 -rpc -srcnodata -9999 -multi -overwrite %s' % (epsg_code, pixel_size, pixel_size, dem_file))

# 关闭数据集
in_dataset = None
out_dataset = None

这段代码使用了GDAL的Python接口实现同样的功能。记得根据实际情况修改输入和输出文件路径、DEM文件路径以及UTM区号和分辨率等参数。