python整形怎么切片_遥感影像切分切片

遥感影像切片

生活当中，我们可能经常遇到处理一个很大的遥感影像情况，例如做一些逐像元运算，不便于并行处理。因此制作了将影像切分成多个小片的程序，切分后保持原有的数值、数据类型、波段数、投影。

切片

用法为

python 此文件的路径 栅格路径 输出文件夹 并行度 输出图像的行列数

例如：

python ./splitImage.py aa.tif ./output 8 1024

splitImage.py内容如下

import numpy as np

import gdal

import os

import sys

from multiprocessing import Pool

NP2GDAL_CONVERSION = {

"uint8": 1,

"int8": 1,

"uint16": 2,

"int16": 3,

"uint32": 4,

"int32": 5,

"float32": 6,

"float64": 7,

"complex64": 10,

"complex128": 11,

}

def split(rasterpath,outpath,dex,i, j,size,cols,rows,gdaltype,datatype,gt,proj,nodata):

dataset = gdal.Open(rasterpath)

if ((i+1)*size > cols) | ((j+1)*size>rows):

#x向越界

if ((i+1)*size > cols) & ((j+1)*size<=rows):

data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, cols-i*size,size)

#y向越界

elif ((i+1)*size <= cols) & ((j+1)*size>rows):

data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, size,rows-j*size)

#xy方向均越界

else:

print(cols-i*size,rows-j*size)

data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, cols-i*size,rows-j*size)

else:

data0 = dataset.GetRasterBand(1).ReadAsArray(i * size, j * size, size,size)

#如果第一个波段的最大值等于最小值，认为是无效值，不对其创建分片

if data0.max() == data0.min():

return

format = "GTiff"

driver = gdal.GetDriverByName(format)

bandsNum = dataset.RasterCount

basename = os.path.basename(rasterpath)

newbasename = basename[:basename.rfind(".")]

dirpath = os.path.dirname(rasterpath)

dst_ds = driver.Create(outpath + newbasename+'_%s%s.tif' % (dex[i], dex[j]), size, size, bandsNum, gdaltype)

gtnew = (gt[0] + i * size * gt[1], gt[1], gt[2], gt[3] + j * size * gt[5], gt[4], gt[5])

dst_ds.SetGeoTransform(gtnew)

dst_ds.SetProjection(proj)

for k in range(bandsNum):

if ((i + 1) * size > cols) | ((j + 1) * size > rows):

# x向越界

if ((i + 1) * size > cols) & ((j + 1) * size <= rows):

data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, cols - i * size, size)

# y向越界

elif ((i + 1) * size <= cols) & ((j + 1) * size > rows):

data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, size, rows - j * size)

# xy方向均越界

else:

data2 = dataset.GetRasterBand(k+1).ReadAsArray(i * size, j * size, cols - i * size, rows - j * size)

smally,smallx=data2.shape

if nodata==None:

data1=np.zeros((size,size),dtype=datatype)

else:

data1 = np.ones((size, size), dtype=datatype)*nodata

data1[0:smally,0:smallx]=data2

else:

data1=dataset.GetRasterBand(k+1).ReadAsArray(i*size,j*size,size,size)

dst_ds.GetRasterBand(k+1).WriteArray(data1)

if nodata != None:

dst_ds.GetRasterBand(k + 1).SetNoDataValue(nodata)

dataset=None

dst_ds=None

if __name__=="__main__":

args=sys.argv

try:

rasterpath = args[1]

outpath = args[2]

parelle = int(args[3])

# 分片的像元行列数，输出为正方形的

size = int(args[4])

except:

print("""Usage:python 此文件的路径 栅格路径 输出文件夹 并行度 输出图像的行列数\nUsage:python ./splitImage.py aa.tif ./output 8 1024""")

raise

# rasterpath = "H:\\testsplit\\myimage.tif"

# outpath = 'H:\\testsplit\\test4'

# parelle=6

# # 分片的像元行列数，输出为正方形的

# size = 2048

# 区分每块的位置，000到999

dex = ["%.3d" % i for i in range(1000)]

if not os.path.exists(outpath):

os.makedirs(outpath)

if not outpath.endswith(os.path.sep):

outpath = outpath + os.path.sep

dataset=gdal.Open(rasterpath)

proj=dataset.GetProjection()

gt=dataset.GetGeoTransform()

datatype=dataset.GetRasterBand(1).ReadAsArray(0,0,1,1).dtype.name

#获取无效值，认为每个波段的无效值是相同的

nodata=dataset.GetRasterBand(1).GetNoDataValue()

#numpy的数据类型，转换为gdal的数据类型

gdaltype=NP2GDAL_CONVERSION[datatype]

#总行列数

cols,rows=dataset.RasterXSize,dataset.RasterYSize

dataset=None

#获取分片的行数和列数

numx,numy=int(np.ceil(cols/size)),int(np.ceil(rows/size))

pool = Pool(parelle)

for i in range(numx):

for j in range(numy):

pool.apply_async(split, args=(rasterpath,outpath,dex,i, j,size,cols,rows,gdaltype,datatype,gt,proj,nodata))

pool.close()

pool.join()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

切分后的文件列表如下

.

├── myimage_000000.tif

├── myimage_000001.tif

├── myimage_000002.tif

├── myimage_000003.tif

├── myimage_000004.tif

├── myimage_000005.tif

├── myimage_001000.tif

├── myimage_001001.tif

├── myimage_001002.tif

├── myimage_001003.tif

├── myimage_001004.tif

├── myimage_001005.tif

├── myimage_002000.tif

├── myimage_002001.tif

├── myimage_002002.tif

├── myimage_002003.tif

├── myimage_002004.tif

├── myimage_002005.tif

├── myimage_003000.tif

├── myimage_003001.tif

└── myimage_003002.tif

拼接

上一步拆分完了，做一些处理，最后可能还需拼接成一个大的影像。这里比较简单，直接调用gdal自带的gdal_merge.py拼接程序。如果电脑安装了python的gdal库的话就会有这个文件，可以通过搜索查找。

使用前首先将gdal_merge.py复制到当前目录下。

这里将一个文件夹下的所有tif文件拼接成一个文件。 具体可以修改tif变量。

用法为

python 此文件的路径 输入文件夹 输出tif nodata值

例如：

python ./mergeImage.py ./output /out.tif 0

mergeImage.py文件内容

import os

import sys

args=sys.argv

if len(args)!=4:

print("usage:python thispy folder_path outtif nodata")

sys.exit(1)

path=args[1]

outtif=args[2]

nodata=eval(args[3])

tifs=[os.path.join(path,i) for i in os.listdir(path) if i.endswith(".tif")]

#print(tifs)

#将 gdal_merge.py复制到当前目录下

os.system("python gdal_merge.py -init %s -n %s -a_nodata %s -o %s %s"%(nodata,nodata,nodata,outtif," ".join(tifs)))1

2

3

4

5

6

7

8

9

10

11

12

13

14