遥感影像分块显示，用于样本制作

遥感影像分块显示，用于样本制作

首先引入gdal和numpy库，这里不过多介绍。自行安装
下面展示一些 python相关代码。

# -*- coding: utf-8 -*-
import os
import numpy
from osgeo import gdal
 
class GRID:
    #读图像文件
    def read_img(self,filename):
        dataset=gdal.Open(filename)       #打开文件
 
        im_width = dataset.RasterXSize    #栅格矩阵的列数
        im_height = dataset.RasterYSize   #栅格矩阵的行数
        print(im_width)
        im_geotrans = dataset.GetGeoTransform()  #仿射矩阵
        im_proj = dataset.GetProjection() #地图投影信息
        im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵
 
        del dataset 
        return im_proj,im_geotrans,im_data
 
    #写文件，以写成tif为例
    def write_img(self,filename,im_proj,im_geotrans,origin_x, origin_y,pixel_width,pixel_height,im_data):
        #gdal数据类型包括
        #gdal.GDT_Byte, 
        #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        #gdal.GDT_Float32, gdal.GDT_Float64
 
        #判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32
 
        #判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1,im_data.shape 
 
        #创建文件
        driver = gdal.GetDriverByName("GTiff")            #数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
        dataset.SetGeoTransform((origin_x, pixel_width, 0, origin_y, 0, pixel_height))              #写入仿射变换参数
        dataset.SetProjection(im_proj)                    #写入投影
 
        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  #写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i+1).WriteArray(im_data[i])
 
        del dataset

# 计算某行列下像元经纬度
def calcLonLat(dataset, x, y):
    minx, xres, xskew, maxy, yskew, yres = dataset.GetGeoTransform()
    lon = minx + xres * x
    lat = maxy +xres * y
    return lon, lat

 
if __name__ == "__main__":

    
    os.chdir(r'H:/chuli')                        #切换路径到待处理图像所在文件夹
    dataset = gdal.Open('image/0305.tif')
    minx, xres, xskew, maxy, yskew, yres = dataset.GetGeoTransform()
    proj,geotrans,data = GRID().read_img('image/0305-3.tif')        #读数据
    print(proj)
    print(geotrans)
    #print(data)
    print(data.shape)
    channel,width,height = data.shape
    for i in range(width//256):#切割成256*256小图
    	for j in range(height//256):
            cur_image=data[:,i*256:(i+1)*256,j*256:(j+1)*256]
            lon=minx + xres * 256 * j
            lat=maxy + yres * (i * 256)
            GRID().write_img('H:/chuli/data/raw1/{}_{}.tif'.format(i, j),proj,geotrans,lon,lat,xres,yres,cur_image)  ##写数据

处理后的每块数据包含原始影像的投影信息以及坐标信息。且能对多波段数据进行处理。处理结果如下图所示：
6个波段的tif影像

真彩色影像分块。

为了便于标注将真彩色tif影像转成png格式。png格式支持RGB三通道或者RGBα四通道。转换代码如下：
下面展示一些 内联代码片。

import os
from osgeo import gdal
 
open_path = "H:/chuli/data/raw3/"
save_path = "H:/chuli/data/png2/"
 
images = os.listdir(open_path)
for image in images:
    im=gdal.Open(os.path.join(open_path,image))
    driver=gdal.GetDriverByName('PNG')
    dst_ds = driver.CreateCopy(os.path.join(save_path,image.split('.')[0]+".png"), im)

结果如下：

本文参考文档：
https://blog.csdn.net/xiaoli_nu/article/details/94064529.
https://blog.csdn.net/weixin_44265800/article/details/109891821.