gdal python 图像金字塔,针对tif,tiff（生成图像金字塔，展示各层级图像）_python tiff层数

作者：在线问答5 | 2024-07-20 03:38:15

踩

python tiff层数

创建图像金字塔

原创不易
导包

import os
from osgeo import gdal
import cv2
import shutil
import skimage.io as io
from PIL import Image
Image.MAX_IMAGE_PIXELS = None
1
2
3
4
5
6
7

判断文件夹是否存在,存在就删除清空，不存在保留

if os.path.exists(self.dir_path):
       	shutil.rmtree(self.dir_path)
        os.makedirs(self.dir_path)
 else:
     	os.makedirs(self.dir_path)
1
2
3
4
5

检查输入和输出路径


class Image_pyramid:
    def __init__(self,inputfpath):
        self.inputtfpath=inputfpath//输入文件路劲
        # print(self.inputtfpath)
        dir_list = self.inputtfpath.strip(' ').split('.')[0].split("/")[0:-1]
        print(dir_list)
        dir_list_pic=self.inputtfpath.strip(' ').split('.')[0].split("/")[-1]
        dir_path2=""
        for i in range(0,len(dir_list)):
            if len(dir_list[i])==0:
                continue
            else:
                dir_list2=str(os.path.join('/',str(dir_list[i])))
                dir_path2=str(dir_path2)+dir_list2
        self.dir_path=os.path.join(dir_path2,dir_list_pic)
        print(self.dir_path)
        if os.path.exists(self.dir_path):
        	shutil.rmtree(self.dir_path)
             os.makedirs(self.dir_path)
    	else:
        	os.makedirs(self.dir_path)
    print(self.dir_path)
    self.outputtfpath = os.path.join(dir_path2,"10")
    print(self.outputtfpath)
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

生成新的tif图像的名称10.tif

def gdal_pymaid2(self):
      if os.path.exists(self.outputtfpath):
          shutil.rmtree(self.outputtfpath)
          os.makedirs(self.outputtfpath)
      else:
          os.makedirs(self.outputtfpath)
      outputfilename = os.path.join(self.outputtfpath, "10.tif")
      print(333333333333333)
      print(outputfilename)
1
2
3
4
5
6
7
8
9
'运行

获取tif图像的宽度和高度


 inputDataset = gdal.Open(self.inputtfpath)
 print(inputDataset==None)
  inputBand = inputDataset.GetRasterBand(1)
  img_arraytif = inputDataset.ReadAsArray(0, 0,
                                     inputDataset.RasterXSize, inputDataset.RasterYSize)
1
2
3
4
5
6

获取原始卫星图像获取仿射矩阵信息，获取投影信息，创建tiff图像驱动，获取波段数量

im_geotrans = inputDataset.GetGeoTransform()  # 获取仿射矩阵信息
im_proj = inputDataset.GetProjection()  # 获取投影信息
driver = gdal.GetDriverByName("GTiff")，#创建tiff图像驱动
im_bands=int(inputDataset.RasterCount)#获取波段数量
print(im_bands)
1
2
3
4
5

将原始卫星图像的矩阵信息宽度信息，高度信息，反射变化矩阵，投影信息，等写入到新的tif图像中

     if im_bands == 1:
            datasets = driver.Create(outputfilename,int((inputDataset.RasterXSize)),
            int((inputDataset.RasterYSize)), 1, inputBand.DataType)
            datasets.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
            datasets.SetProjection(im_proj)  # 写入投影
      else:
            datasets = driver.Create(outputfilename, int((inputDataset.RasterXSize)), int((inputDataset.RasterYSize)), 3, inputBand.DataType)
            datasets.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
            datasets.SetProjection(im_proj)  # 写入投影
1
2
3
4
5
6
7
8
9

卫星图像可能有一个波段，3个波段，甚至7个波段的信息，将每个波段的信息写入到对应的tif图像中，在这之前需要区分RGB通道，可以利用均值直方图信息，这里将图像一次缩小4倍，16倍，64倍…

 if im_bands >= 3:
            # for i in range(3):
            datasets.GetRasterBand(1).WriteArray(img_arraytif[2])
            datasets.GetRasterBand(2).WriteArray(img_arraytif[1])
            datasets.GetRasterBand(3).WriteArray(img_arraytif[0])
            for i in range(1, 4):
                datasets.GetRasterBand(i).ComputeStatistics(False)
                datasets.BuildOverviews('average', [1,2, 4, 8, 16, 32, 64])
 else:
            datasets.GetRasterBand(1).WriteArray(img_arraytif)
            datasets.GetRasterBand(1).ComputeStatistics(False)
            datasets.BuildOverviews('average', [1,2, 4, 8, 16, 32,64])
 del datasets
1
2
3
4
5
6
7
8
9
10
11
12
13

cv写入不同层数的图像金字塔所对应的图像

def show_pymaid(self):
    self.gdal_pymaid2()

    outputfilename = os.path.join(self.outputtfpath, "10.tif")
    
1
2
3
4
5
'运行

读取新创建的卫星图像基本信息，获取概率图深度，即金字塔层数和波段数量

 dataset = gdal.Open(outputfilename)
 # print(dataset)

 band = dataset.GetRasterBand(1)
 overviewNum = band.GetOverviewCount()

 print(dataset.RasterCount)
 print(overviewNum)
 im_bands = int(dataset.RasterCount)
1
2
3
4
5
6
7
8
9

不同波段的卫星图像，图像金字塔所对应的原始图像大小

if im_bands>1:
	rband = dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(6)
	r = roverviewBand.ReadAsArray()
	gband = dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(6)
	g = goverviewBand.ReadAsArray()
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(6)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([r, g, b])#正解
	# data = cv2.merge([b, g, r])
	print(data.shape)
	cv2.imwrite(f'{self.dir_path}/1.bmp', data)
	rband=dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(5)
	r = roverviewBand.ReadAsArray()
	gband=dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(5)
	g = goverviewBand.ReadAsArray()
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(5)
	b = broverviewBand.ReadAsArray()
	data=cv2.merge([r,g,b])
	# data = cv2.merge([b, g, r])
	print(data.shape)
	cv2.imwrite(f'{self.dir_path}/1.bmp',data)
	del data
	
	rband=dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(4)
	r = roverviewBand.ReadAsArray()
	gband=dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(4)
	g = goverviewBand.ReadAsArray()
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(4)
	b = broverviewBand.ReadAsArray()
	data=cv2.merge([r,g,b])
	# data = cv2.merge([b, g, r])
	print(data.shape)
	
	cv2.imwrite(f'{self.dir_path}/2.bmp',data)
	
	del data
	
	
	rband=dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(3)
	r = roverviewBand.ReadAsArray()
	gband=dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(3)
	g = goverviewBand.ReadAsArray()
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(3)
	b = broverviewBand.ReadAsArray()
	data=cv2.merge([r,g,b])
	# data = cv2.merge([b, g, r])
	print(data.shape)
	
	cv2.imwrite(f'{self.dir_path}/3.bmp',data)
	
	del data
	#
	rband=dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(2)
	r = roverviewBand.ReadAsArray()
	gband=dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(2)
	g = goverviewBand.ReadAsArray()
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(2)
	b = broverviewBand.ReadAsArray()
	data=cv2.merge([r,g,b])
	# data = cv2.merge([b, g, r])
	
	cv2.imwrite(f'{self.dir_path}/4.bmp',data)
	
	del data
	#
	#
	rband=dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(1)
	r = roverviewBand.ReadAsArray()
	gband=dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(1)
	g = goverviewBand.ReadAsArray()
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(1)
	b = broverviewBand.ReadAsArray()
	data=cv2.merge([r,g,b])
	# data = cv2.merge([b, g, r])
	print(data.shape)
	
	cv2.imwrite(f'{self.dir_path}/5.bmp',data)
	
	del data
	
	
	rband=dataset.GetRasterBand(3)
	roverviewBand = rband.GetOverview(0)
	r = roverviewBand.ReadAsArray()
	gband=dataset.GetRasterBand(2)
	goverviewBand = gband.GetOverview(0)
	g = goverviewBand.ReadAsArray()
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(0)
	b = broverviewBand.ReadAsArray()
	import numpy as np
	data=cv2.merge([r,g,b])
	# data = cv2.merge([b, g, r])
	print(data.shape)
	# io.imsave(f'{self.dir_path}/6.bmp',data)
	cv2.imwrite(f'{self.dir_path}/6.bmp',data)
	
	
	del data
else:
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(6)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([b, b, b])
	cv2.imwrite(f'{self.dir_path}/1.bmp', data)
	
	del data
	
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(5)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([b, b, b])
	
	cv2.imwrite(f'{self.dir_path}/1.bmp',data)
	
	del data
	
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(4)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([b, b, b])
	print(data.shape)
	
	cv2.imwrite(f'{self.dir_path}/2.bmp',data)
	
	del data
	
	
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(3)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([b, b, b])
	
	cv2.imwrite(f'{self.dir_path}/3.bmp',data)
	
	del data
	
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(2)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([b, b, b])
	
	cv2.imwrite(f'{self.dir_path}/4.bmp',data)
	
	del data
	
	bband = dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(1)
	b = broverviewBand.ReadAsArray()
	data = cv2.merge([b, b, b])
	
	cv2.imwrite(f'{self.dir_path}/5.bmp',data)
	
	del data
	
	bband=dataset.GetRasterBand(1)
	broverviewBand = bband.GetOverview(0)
	b = broverviewBand.ReadAsArray()
	data=cv2.merge([b,b,b])
	print(data.shape)
	
	cv2.imwrite(f'{self.dir_path}/6.bmp',data)
	
	# del data
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
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/在线问答5/article/detail/854962