Python遥感开发之GDAL读写遥感影像

前言：主要介绍了使用GDAL读写遥感影像数据的操作，包括读取行、列、投影、值以及数据的简单运算和生成新的tif影像。

---

1 读取tif信息方法一

from osgeo import gdal
import numpy as np

if __name__ == '__main__':
    dataset = gdal.Open("lucc.tif")#读取的是某地的土地类型

    col = dataset.RasterXSize  # 图像长度
    print("col:",col)

    row = dataset.RasterYSize  # 图像宽度
    print("row:", row)

    geotrans = dataset.GetGeoTransform()  # 读取仿射变换
    print("geotrans:", geotrans)

    proj = dataset.GetProjection()  # 读取投影
    print("proj:", proj)

    # num_bands = dataset.RasterCount  # 查看波段个数，单波段默认是1
    # print(num_bands)
    # data_band = dataset.GetRasterBand(1)  # 1波段的具体内容
    # print(data_band.ReadAsArray())

    data = dataset.ReadAsArray()  # 转为numpy格式
    data = data.astype(np.float32)
    a = data[0][0]
    data[data == a] = np.nan
    print("data:", data)

    #遍历每一行像元值
    for i in range(0,row):
        print(i,data[i])

    #遍历读取每一个像元
    for i in range(0,row):
        for j in range(0,col):
            if not np.isnan(data[i][j]):#筛选有效值
                print(data[i][j])
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

2 读取tif信息方法二

from osgeo import gdalnumeric
import numpy as np

if __name__ == '__main__':
    data = gdalnumeric.LoadFile("lucc.tif")
    data = data.astype(np.float32)
    a = data[0][0]
    data[data == a] = np.nan
    #遍历每一行像元
    for d in data:
        print(d)
    #遍历每一个像元
    for d in data:
        for s in d:
            print(s)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

3 自己封装读取tif的方法（推荐）

import numpy as np
from osgeo import gdal,gdalnumeric

def read_tif01(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)#转为float类型
    a = data[0][0]
    data[data == a] = np.nan #原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值
    return [col, row, geotrans, proj, data]

def read_tif02(filepath):
    data = gdalnumeric.LoadFile(filepath)
    data = data.astype(np.float32)
    a = data[0][0]
    data[data == a] = np.nan
    return data

if __name__ == '__main__':
    col, row, geotrans, proj, data = read_tif01("lucc.tif")
    data2 = read_tif02("lucc.tif")
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

4 对读取的tif数据进行简单运算

import numpy as np
from osgeo import gdal

def read_tif01(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)#转为float类型
    a = data[0][0]
    data[data == a] = np.nan #原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值
    return [col, row, geotrans, proj, data]

if __name__ == '__main__':
    col, row, geotrans, proj, data = read_tif01("lucc.tif")
    print(data[1])
    data = data*2 #在原来的data基础上所有值乘以2
    print(data[1])

    #可以进行条件筛选
    data[data==6] = 10#所有像元值为6的重新赋值为10
    print(data[1])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

5 写出tif影像(推荐)

import numpy as np
from osgeo import gdal

def read_tif01(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize#图像长度
    row = dataset.RasterYSize#图像宽度
    geotrans = dataset.GetGeoTransform()#读取仿射变换
    proj = dataset.GetProjection()#读取投影
    data = dataset.ReadAsArray()#转为numpy格式
    data = data.astype(np.float32)#转为float类型
    a = data[0][0]
    data[data == a] = np.nan #原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值
    return [col, row, geotrans, proj, data]

def save_tif(data, file, output):
    ds = gdal.Open(file)
    shape = data.shape
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Float32)#以float类型进行存储
    dataset.SetGeoTransform(ds.GetGeoTransform())
    dataset.SetProjection(ds.GetProjection())
    dataset.GetRasterBand(1).WriteArray(data)

if __name__ == '__main__':
    col, row, geotrans, proj, data = read_tif01("lucc.tif")
    data = data*2 #在原来的data基础上所有值乘以2
    #生成新的tif
    save_tif(data,"lucc.tif","new_lucc.tif")#可以自己指定文件目录
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

6.读取多波段遥感数据

如图所示，此TIF有绿光、蓝光、红光、近红外四个波段，需要同时读取四个波段的数据，代码如下：

import numpy as np
from osgeo import gdal

def read_tif01(filepath):
    dataset = gdal.Open(filepath)
    col = dataset.RasterXSize  # 图像长度
    row = dataset.RasterYSize  # 图像宽度
    geotrans = dataset.GetGeoTransform()  # 读取仿射变换
    proj = dataset.GetProjection()  # 读取投影
    data = dataset.ReadAsArray()  # 转为numpy格式
    data = data.astype(np.float32)  # 转为float类型
    a = data[0][0]
    data[data == a] = np.nan  # 原因：读取某一个行政区的影像图的时候，往往第一行的第一列值为空值
    return [col, row, geotrans, proj, data]


def save_tif(data, file, output):
    ds = gdal.Open(file)
    shape = data.shape
    driver = gdal.GetDriverByName("GTiff")
    dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Float32)  # 以float类型进行存储
    dataset.SetGeoTransform(ds.GetGeoTransform())
    dataset.SetProjection(ds.GetProjection())
    dataset.GetRasterBand(1).WriteArray(data)

if __name__ == '__main__':
    sentinel_tif_path = r"C:\Users\15739\Desktop\随机森林 (2)\测试数据.tif"
    col, row, geotrans, proj, data = read_tif01(sentinel_tif_path)

    b1_list = data[0]
    b2_list = data[1]
    b3_list = data[2]
    b4_list = data[3]
    for i in range(0,row):
        for j in range(0,col):
            if np.isnan(b1_list[i][j]) or np.isnan(b2_list[i][j]) or np.isnan(b3_list[i][j]) or np.isnan(b4_list[i][j]):
                pass
            else:
                list1 = []
                list1.append(b1_list[i][j])
                list1.append(b2_list[i][j])
                list1.append(b3_list[i][j])
                list1.append(b4_list[i][j])
                print(list1)
    print(data.shape)
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
44
45