遥感影像处理9.21作业

本文最后更新于:2 months ago

用的是Python语言和gdal库,效果为

数据

来源是USGS,位置是北京,可以看到西北有张家口,我的家乡~

导入gdal包

下载gdal包还有numpy,因为我的python是3.9版本,所以直接pip install numpy老出错,所以,可以直接进入这个网站
然后在win+r中输入cmd打开命令行,再进入Python的包的安装文件夹,pip install +轮子的文件名,可以看到安装成功,如下图所示

GDAL同理

导入包和指定landset文件夹所在位置

1
2
3
4
5
6
7
8
import os
from osgeo import gdal

# 当前所在路径
os.chdir(r'C:\Users\17503\Desktop\landset_data')
band1_fn = 'LC08_L1TP_123032_20220910_20220914_02_T1_B1.TIF'
band2_fn = 'LC08_L1TP_123032_20220910_20220914_02_T1_B2.TIF'
band3_fn = 'LC08_L1TP_123032_20220910_20220914_02_T1_B3.TIF'

打开蓝波段,在创建输出图像之前,你需要此波段对象,因为它具有输出文件所需要的基本信息(尺寸大小、投影、坐标等)。

1
2
in_ds = gdal.Open(band1_fn)
in_band = in_ds.GetRasterBand(1)

使用GeoTIFF驱动程序创建新数据集

1
2
3
4
5
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('nat_color.tif',
in_band.XSize, in_band.YSize, 3, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())

创建数据集

1
driver.Create(filename, xsize, ysize, [bands], [data_type], [options])
  • filename->’nat_color.tif’是要创建的数据集在当前文件夹下的路径。
  • xsize->in_band.XSize是新数据集中的列数。
  • ysize->in_band.YSize是新数据集中的行数。
  • 3是新数据集中的波段数,可以不写。默认值为1。
  • in_band.DataType是将存储在新数据集中的数据类型。默认值为GDT_Byte。

读取第1波段数据

因为Landsat图像的波段1是蓝色波段,你需要将其放入输出图像的第三个波段以获得RGB顺序的波段。接下来要做的是从out_ds获取第三个波段,然后使用WriteArray将in_data数组中的值复制到新数据集的第三个波段

1
2
3
in_data = in_band.ReadAsArray()
out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)

读取第2波段数据

1
2
3
in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())

读取第3波段数据

1
2
out_ds.GetRasterBand(1).WriteArray(
gdal.Open(band3_fn).ReadAsArray())

统计像元

在计算统计数据之前,你必须确保数据已写入磁盘而不是仅缓存在内存中,因此这是对FlushCache的调用。然后循环遍历波段并计算每个波段的统计数据。将False传递给此函数会告诉你需要实际统计信息而不是估计值,它可能来自概览图层(尚不存在)或从像元子集中采样。如果估计值可以接受,那么你可以传递True;这也将使计算更快,因为不是每个像元都需要检查。

1
2
3
4
out_ds.FlushCache()
for i in range(1, 4):
data=out_ds.GetRasterBand(i).ComputeStatistics(False)
print(data)

分别是最小值,最大值,平均值,均方差

抽样构建金字塔

遥感领域许多图片的大小动辄上G。读取、显示这样的图片极为耗时,影响用户体验。金字塔技术在几乎不降低显示效果的前提下,大大降低了图片处理的耗时,改善了用户体验Tutorial

1
2
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
del out_ds