【python】【gdal】读取栅格数据转换成数组，数组转换栅格数据(tif文件)_数组转栅格

every blog every motto: The shortest answer is doing.

0. 前言

本节主要讲解有关利用python、gdal读取栅格数据(tif)，再将读取后的转换成数组。
说明： 添加了直接读取dataset，更加简洁（2021.5.21）

1. 正文

1.1. 读取tif图片,并转换成数组

1.1.1 逐通道读取

读取栅格数据转换成数组，其中需要注意的是，栅格数据波段1开始，数组通道从0开始。

 def load_img_to_array(self, img_file_path):
     """
     读取栅格数据，将其转换成对应数组
     :param: img_file_path: 栅格数据路径
     :return: 返回投影，几何信息，和转换后的数组(5888,5888,10)
     """
     dataset = gdal.Open(img_file_path)  # 读取栅格数据
     print('处理图像的栅格波段数总共有：', dataset.RasterCount)

     # 判断是否读取到数据
     if dataset is None:
         print('Unable to open *.tif')
         sys.exit(1)  # 退出

     projection = dataset.GetProjection()  # 投影
     transform = dataset.GetGeoTransform()  # 几何信息

     ''' 栅格图片转成数组,依次读取每个波段，然后转换成数组的每个通道 '''

     array_channel = 0  # 数组从通道从0开始
     print('处理其中的波段数：', self.Bands)

     # 读取Bands列表中指定的波段
     for band in self.Bands:  # 图片波段从1开始
         print('正在将其中的{}波段转换成数组'.format(band))
         srcband = dataset.GetRasterBand(band)  # 获取栅格数据集的波段
         if srcband is None:
             print('WARN: srcband is None: ' + str(band) + img_file_path)
             continue

         #  一个通道转成一个数组（5888,5888）
         arr = srcband.ReadAsArray()

         if array_channel == 0:  # 如果是第一个通道，要新建数组（5888,5888,10）
             H = arr.shape[0]  # 图像的长 5888
             W = arr.shape[1]  # 图像的宽 5888
             img_array = np.zeros((H, W, len(self.Bands)), dtype=np.float32)

         img_array[:, :, array_channel] = np.float32(arr)

         array_channel += 1  # 进入下一个波段-> 数组

     return projection, transform, img_array
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

1.1.2 一次性读取

 def load_img_to_array(self, img_file_path):
     """
     读取栅格数据，将其转换成对应数组
     :param: img_file_path: 栅格数据路径
     :return: 返回投影，几何信息，和转换后的数组(5888,5888,10)
     """
     dataset = gdal.Open(img_file_path)  # 读取栅格数据
     print('处理图像的栅格波段数总共有：', dataset.RasterCount)

     # 判断是否读取到数据 
     if dataset is None:
         print('Unable to open *.tif')
         sys.exit(1)  # 退出

     projection = dataset.GetProjection()  # 投影
     transform = dataset.GetGeoTransform()  # 几何信息

    # 直接读取dataset
    img_array = dataset.ReadAsArray()

     return projection, transform, img_array
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

1.2 将数组转换成tif图片

依次读取数组的每个通道，将每个通道写成波段，最后保存在硬盘中。
说明： 下方代码是因为对图片进行了相关处理，几何信息发生变换，所以要对其进行调整。如果图片没有发生变化不需要进行调整

    def predit_to_tif(self, mat, projection, tran, mapfile):
        """
        将数组转成tif文件写入硬盘
        :param mat: 数组
        :param projection: 投影信息
        :param tran: 几何信息
        :param mapfile: 文件路径
        :return:
        """

        row = mat.shape[0]  # 矩阵的行数
        columns = mat.shape[1]  # 矩阵的列数

        if self.isSeg == 1:
            x_res = tran[1] * 1  # x方向分辨率
            y_res = tran[5] * 1  # y方向分辨率

            geo_transform = (tran[0] + x_res * (1 - 1) / 2.0, x_res, 0, tran[3] + y_res * (1 - 1) / 2.0, 0, y_res)
        else:
            x_res = tran[1] * self.step
            y_res = tran[5] * self.step

            geo_transform = (
                tran[0] + tran[1] * (self.dim_x - 1) / 2.0, x_res, 0, tran[3] + tran[5] * (self.dim_y - 1) / 2.0, 0,
                y_res)

        print(geo_transform)

        dim_z = mat.shape[2]  # 通道数

        driver = gdal.GetDriverByName('GTiff')  # 创建驱动
        # 创建文件
        dst_ds = driver.Create(mapfile, columns, row, dim_z, gdal.GDT_UInt16)
        dst_ds.SetGeoTransform(geo_transform)  # 设置几何信息
        dst_ds.SetProjection(projection)  # 设置投影

        # 将数组的各通道写入tif图片
        for channel in np.arange(dim_z):
            map = mat[:, :, channel]
            dst_ds.GetRasterBand(int(channel + 1)).WriteArray(map)

        dst_ds.FlushCache()  # 写入硬盘
        dst_ds = None
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

简洁版：

def to_tif(arr, tran=(2037.0, 10.0, 0.0, 11238.000000000002, 0.0, -10.0)):
    """ 将数组保存成tif"""

    # 数组shape
    row = arr.shape[0]  # 行数
    columns = arr.shape[1]  # 列数
    dim = 1  # 通道数

    # 创建驱动
    driver = gdal.GetDriverByName('GTiff')
    # 创建文件
    dst_ds = driver.Create('./data/2.tif', columns, row, dim, gdal.GDT_UInt16)
    # 设置几何信息
    dst_ds.SetGeoTransform(tran)

    # 将数组写入
    dst_ds.GetRasterBand(1).WriteArray(arr)
    # 写入硬盘
    dst_ds.FlushCache()
    dst_ds = None

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

结果：

1.3 数据类型对照

参考文献

[1] http://blog.sina.com.cn/s/blog_132ecc3670102xzhr.html
[2] https://www.cnblogs.com/Romi/archive/2012/03/14/2396627.html
[3] https://gdal.org/tutorials/raster_api_tut.html
[4] http://blog.sina.com.cn/s/blog_b36001e10102x7bu.html
[5] https://blog.csdn.net/X_Cosmic/article/details/107182061