python学习——实现不同大小两景影像的合成

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

 直奔主题，这次要实现的是将一个有随机数生成的栅格和一个目标栅格进行合成，最终为一个栅格。

第一景，坐标系为WGS-84，像元大小为1KM，换算成度为0.009009009

第二景，坐标系为WGS-84,像元大小为1KM

 

import gdal,ogr
import numpy as np
def rasterlayer(radom_raster,dust_path):
    dust_raster = dust_path+'\\'+'dust_raster.tif'
    dust_ds = gdal.Open(dust_raster)
    radom_ds = gdal.Open(radom_raster)

    dust_geo_transform = dust_ds.GetGeoTransform()
    dust_cols=dust_ds.RasterXSize
    dust_rows=dust_ds.RasterYSize
    dust_x_min = dust_geo_transform[0] #获取左至坐标
    dust_y_max = dust_geo_transform[3] #获取上至坐标

    #对随机数栅格进行操作
    radom_geo_transform = radom_ds.GetGeoTransform()
    radom_cols=radom_ds.RasterXSize
    radom_rows=radom_ds.RasterYSize
    radom_x_min = radom_geo_transform[0] #获取左至坐标
    radom_y_max = radom_geo_transform[3] #获取上至坐标


    # 进行沙尘栅格和随机数栅格的叠加
    dust_band = dust_ds.GetRasterBand(1)
    dust_data = dust_band.ReadAsArray()
    radom_band = radom_ds.GetRasterBand(1)
    radom_data = radom_band.ReadAsArray()
    # 沙尘栅格在随机数栅格中的下行号
    down_row = round((radom_y_max - dust_y_max) / 0.009009009)
    # 沙尘栅格在随机数栅格中的上行号
    up_row = down_row+dust_rows
    # 沙尘栅格在随机数栅格中的左列号
    left_col = round((dust_x_min - radom_x_min) / 0.009009009)
    # 沙尘栅格在随机数栅格中的右列号
    right_col = left_col+dust_cols
    print(down_row,up_row,left_col,right_col)
    #制作一个掩膜文件
    mask1 = dust_data == 0
    #使用np.choose将大栅格中小栅格所在的所有行和列替换成小栅格的值
    radom_data[down_row:up_row, left_col:right_col] \
        = np.choose(mask1, (dust_data, radom_data[down_row:up_row, left_col:right_col]))
    #之后进行输出。
    radom_proj = radom_ds.GetProjection()
    driver = gdal.GetDriverByName("GTiff")  
    subset = driver.Create(dust_path+'\\'+'dust_random.tif', 6837, 3932, 1, gdal.GDT_Float32)
    ''':type:osgeo.gdal.Dataset'''
    inband = subset.GetRasterBand(1)
    subset_array = subset.GetRasterBand(1)
    subset_array.WriteArray(radom_data)
    subset_array.FlushCache()
    transform = [radom_x_min, radom_geo_transform[1], 0, radom_y_max, 0,
                 -radom_geo_transform[1]]
    subset.SetGeoTransform(transform)
    subset.SetProjection(radom_proj)


if __name__ == '__main__':
    radom_raster = r'E:\JS\Temp\random_raster.dat'#随机数栅格
    dust_path = r'G:\rasterdust'   #行列数较小栅格
    rasterlayer(radom_raster,dust_path)

 结果图：

 

 

 如有错误，欢迎指正，万分感谢！