处理遥感影像通常涉及以下步骤:
环境准备
确保已安装Python和GDAL库。
可以通过`pip`安装GDAL库:`pip install gdal`。
读取遥感影像
使用`rasterio`或`gdal`库读取影像文件。
示例代码:
```python
from osgeo import gdal
file_path = "path_to_your_image.tif"
dataset = gdal.Open(file_path)
im_data = dataset.ReadAsArray()
im_width = dataset.RasterXSize
im_height = dataset.RasterYSize
im_geotrans = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
处理影像
对影像进行裁剪、重采样、滤波等操作。
示例代码:
```python
假设对影像进行简单的阈值处理
threshold = 100
processed_data = im_data > threshold
计算年度平均影像
读取多个影像文件,计算每年的平均影像。
示例代码:
```python
def calculate_yearly_mean(input_folder, output_folder):
yearly_data = {}
for filename in os.listdir(input_folder):
if filename.endswith('.tif'):
year = filename.split('_')[-1].split('.')
data = gdal.Open(os.path.join(input_folder, filename)).ReadAsArray()
yearly_data[year] = data
for year, data in yearly_data.items():
yearly_data[year] = np.mean(data, axis=(0, 1))
for year, data in yearly_data.items():
output_file = os.path.join(output_folder, f"{year}_mean.tif")
driver = gdal.GetDriverByName('GTiff')
new_dataset = driver.Create(output_file, im_width, im_height, 1, gdal.GDT_Float32)
new_dataset.GetRasterBand(1).WriteRaster(0, 0, im_width, im_height, data.tobytes())
影像镶嵌拼接
使用`rasterio`和`gdal`库进行影像的镶嵌拼接。
示例代码:
```python
def mosaic_images(input_folder, output_file):
images = [os.path.join(input_folder, f) for f in os.listdir(input_folder) if f.endswith('.tif')]
datasets = [gdal.Open(img) for img in images]
cols = [ds.RasterXSize for ds in datasets]
rows = [ds.RasterYSize for ds in datasets]
min_cols = min(cols)
min_rows = min(rows)
driver = gdal.GetDriverByName('GTiff')
new_dataset = driver.Create(output_file, min_cols, min_rows, len(images), gdal.GDT_Float32)
for i, ds in enumerate(datasets):
new_dataset.GetRasterBand(i + 1).WriteRaster(0, 0, min_cols, min_rows, ds.ReadAsArray())
for ds in datasets:
ds = None
影像辐射定标
对影像进行辐射定标,例如将像素值除以10000。
示例代码:
```python
def radiometric_calibration(tif_path, output_path):
dataset = gdal.Open(tif_path)
data = dataset.ReadAsArray()
calibrated_data = data / 10000.0
driver = gdal.GetDriverByName('GTiff')
new_dataset = driver.Create(output_path, dataset.RasterXSize, dataset.RasterYSize, 1, gdal.GDT_Float32)
new_dataset.GetRasterBand(1).WriteRaster(0, 0, dataset.RasterXSize, dataset.RasterYSize, calibrated_data.tobytes())
修复遥感影像条带
使用GDAL修复由于卫星传感器故障导致的影像条带问题。
示例代码: