遥感图切割

使用GDAL包可以进行遥感图的处理，使用ENVI工具可以方便查看遥感图像
分割遥感图，保存成tif格式的，代码如下：

# -*- 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  # 栅格矩阵的行数

        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, 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(im_geotrans)  # 写入仿射变换参数
        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


if __name__ == "__main__":
    # os.chdir(r'E:/data')  # 切换路径到待处理图像所在文件夹
    proj, geotrans, data = GRID().read_img('D:/Image/NECAS/20200203/Ortho/Ortho.tif')  # 读数据
    print(proj)
    print(geotrans)
    # print(data)
    print(data.shape)
    channel, width, height = data.shape
    print(width, height)
    for i in range(width // 4000):  # 切割成4000*3000小图
        for j in range(height // 3000):
            cur_image = data[:, i * 4000:(i + 1) * 4000, j * 3000:(j + 1) * 3000]
            # driver = gdal.GetDriverByName('JPEG')
            # dst_ds = driver.CreateCopy('D:/Image/NECAS/20200203/cut2/{}_{}.jpg'.format(i, j), cur_image)
            GRID().write_img('D:/Image/NECAS/20200203/cut/{}_{}.tif'.format(i, j), proj, geotrans, cur_image)  ##写数据

将遥感图分割后，保存成jpg格式的，代码如下：

# -*- coding: utf-8 -*-
import os
import  numpy  as np
from osgeo import gdal
import cv2


def readTif(fileName):
    merge_img = 0
    driver = gdal.GetDriverByName('GTiff')
    driver.Register()

    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName + "掩膜失败，文件无法打开")
        return
    im_width = dataset.RasterXSize  # 栅格矩阵的列数
    print('im_width:', im_width)

    im_height = dataset.RasterYSize  # 栅格矩阵的行数
    print('im_height:', im_height)
    im_bands = dataset.RasterCount  # 波段数
    im_geotrans = dataset.GetGeoTransform()  # 获取仿射矩阵信息
    im_proj = dataset.GetProjection()  # 获取投影信息
    print(im_bands)

    if im_bands == 1:
        band = dataset.GetRasterBand(1)
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 获取数据
        cdata = im_data.astype(np.uint8)
        merge_img = cv2.merge([cdata, cdata, cdata])

        cv2.imwrite('C:/Users/summer/Desktop/a.jpg', merge_img)
    #
    elif im_bands == 4:
        band1=dataset.GetRasterBand(1)
        band2=dataset.GetRasterBand(2)
        band3=dataset.GetRasterBand(3)
        band4=dataset.GetRasterBand(4)
        for i in range(im_width // 4000):  # 切割成4000*3000小图
            for j in range(im_height // 3000):
                data1 = band1.ReadAsArray(i * 4000,  j * 3000,  4000,3000).astype(np.uint8)  # r #获取数据
                data2 = band2.ReadAsArray(i * 4000,  j * 3000, 4000,3000).astype(np.uint8)  # g #获取数据
                data3 = band3.ReadAsArray(i * 4000,  j * 3000, 4000,3000).astype(np.uint8)  # b #获取数据
                data4 = band4.ReadAsArray(i * 4000,  j * 3000,  4000,3000).astype(np.uint8)  # R #获取数据
                # print(data1[1][45])
                output1= cv2.convertScaleAbs(data1)#alpha=(255.0/65535.0)
                # print(output1[1][45])
                output2= cv2.convertScaleAbs(data2)
                output3= cv2.convertScaleAbs(data3)

                merge_img1 = cv2.merge([output3, output2, output1])  # B G R

                cv2.imwrite('D:/Image/NECAS/20200203/cut2/{}_{}.jpg'.format(i, j), merge_img1)
                print("success")

if  __name__=='__main__':
    readTif("D:/Image/NECAS/20200203/Ortho/Ortho.tif")
    print ("0k")

详细怎么使用GDAL 自行百度。
参考链接：http://www.dengb.com/Pythonjc/1318700.html
https://www.osgeo.cn/pygis/gdal-gdalreadata.html
https://blog.csdn.net/xiaoli_nu/article/details/94064529