GDAL库学习笔记(二):GDAL库的一些细节
1. GDAL库的一些细节 1.1. 关于ReadRaster 1.1.1. 缩放 下面是我的一个小例子:我使用了ReadAsArray,返回的是Numeric模块定义的Array,我喜欢它的是因为它排列很工整,看起来很好看。 (这里的数据是用GRASS从spearfish示例数据转出的aspect数据集)...
- 作者:lilin来源:啄木鸟CPUG站|2007年02月08日
1. GDAL库的一些细节
1.1. 关于ReadRaster
1.1.1. 缩放
下面是我的一个小例子:我使用了ReadAsArray,返回的是Numeric模块定义的Array,我喜欢它的是因为它排列很工整,看起来很好看。 (这里的数据是用GRASS从spearfish示例数据转出的aspect数据集)
1
>>
>
import
gdal
2
>>
>
from
gdalconst
import
*
3
>>
>
dataset
=
gdal
.
Open
(
"f:/pyproj/gtif/aspect.tif"
)
4
>>
>
band
=
dataset
.
GetRasterBand
(
1
)
5
>>
>
band
.
ReadAsArray
(
100
,
100
,
5
,
5
,
10
,
10
)
6
array
(
[
[
61
,
61
,
58
,
58
,
90
,
90
,
87
,
87
,
45
,
45
]
,
7
[
61
,
61
,
58
,
58
,
90
,
90
,
87
,
87
,
45
,
45
]
,
8
[
36
,
36
,
59
,
59
,
113
,
113
,
88
,
88
,
37
,
37
]
,
9
[
36
,
36
,
59
,
59
,
113
,
113
,
88
,
88
,
37
,
37
]
,
10
[
255
,
255
,
82
,
82
,
135
,
135
,
72
,
72
,
22
,
22
]
,
11
[
255
,
255
,
82
,
82
,
135
,
135
,
72
,
72
,
22
,
22
]
,
12
[
45
,
45
,
129
,
129
,
144
,
144
,
255
,
255
,
255
,
255
]
,
13
[
45
,
45
,
129
,
129
,
144
,
144
,
255
,
255
,
255
,
255
]
,
14
[
90
,
90
,
110
,
110
,
98
,
98
,
35
,
35
,
45
,
45
]
,
15
[
90
,
90
,
110
,
110
,
98
,
98
,
35
,
35
,
45
,
45
]
]
,
'b'
)
16
>>
>
band
.
ReadAsArray
(
100
,
100
,
10
,
10
,
10
,
10
)
17
array
(
[
[
61
,
58
,
90
,
87
,
45
,
255
,
255
,
117
,
65
,
50
]
,
18
[
36
,
59
,
113
,
88
,
37
,
26
,
63
,
165
,
23
,
74
]
,
19
[
255
,
82
,
135
,
72
,
22
,
29
,
67
,
118
,
77
,
120
]
,
20
[
45
,
129
,
144
,
255
,
255
,
36
,
94
,
108
,
88
,
97
]
,
21
[
90
,
110
,
98
,
35
,
45
,
88
,
109
,
121
,
92
,
73
]
,
22
[
94
,
108
,
85
,
45
,
55
,
97
,
144
,
167
,
135
,
21
]
,
23
[
96
,
106
,
82
,
45
,
45
,
45
,
255
,
230
,
251
,
255
]
,
24
[
246
,
236
,
255
,
255
,
3
,
255
,
255
,
247
,
254
,
255
]
,
25
[
236
,
249
,
255
,
255
,
255
,
255
,
255
,
255
,
255
,
255
]
,
26
[
194
,
212
,
255
,
255
,
255
,
255
,
255
,
255
,
255
,
255
]
]
,
'b'
)
27
>>
>
band
.
ReadAsArray
(
100
,
100
,
20
,
20
,
10
,
10
)
28
array
(
[
[
54
,
95
,
91
,
150
,
53
,
135
,
164
,
139
,
35
,
37
]
,
29
[
128
,
152
,
86
,
97
,
96
,
91
,
116
,
199
,
255
,
200
]
,
30
[
101
,
66
,
71
,
135
,
80
,
152
,
255
,
255
,
255
,
210
]
,
31
[
171
,
159
,
87
,
247
,
254
,
202
,
104
,
255
,
255
,
160
]
,
32
[
223
,
255
,
255
,
255
,
255
,
206
,
147
,
193
,
218
,
121
]
,
33
[
227
,
255
,
255
,
116
,
150
,
238
,
255
,
255
,
137
,
162
]
,
34
[
230
,
255
,
231
,
247
,
145
,
156
,
134
,
30
,
130
,
136
]
,
35
[
155
,
196
,
252
,
230
,
187
,
19
,
134
,
195
,
126
,
144
]
,
36
[
80
,
54
,
85
,
105
,
163
,
140
,
192
,
95
,
154
,
146
]
,
37
[
150
,
83
,
71
,
161
,
173
,
255
,
255
,
120
,
149
,
180
]
]
,
'b'
)
38
>>
>
ReadAsArray的原型
1
>>
>
help
(
band
.
ReadAsArray
)
2
Help
on
method
ReadAsArray
in
module
gdal
:
3
ReadAsArray
(
self
,
xoff
=
0
,
yoff
=
0
,
win_xsize
=
None
,
win_ysize
=
None
,
buf_xsize
=
None
4
,
buf_ysize
=
None
,
buf_obj
=
None
)
method
of
gdal
.
Band
instance
看看函数的几个参数的意义。头两个100是取值窗口的左上角在实际数据中所处象元的xy位置。后两个是取值窗口覆盖的区域大小,再后面 是取值窗口取出数组进行缩放后数组的大小。这里需要注意的是这里的buffer大小是根据参数自动分配的,可以不指定,如果不指定,则和第3,4个参数一致。经过5,6两个参数的设置,可以进行缩放。比如真的取值窗口大小可以是20*20,而取出后数组就可以人工设置大小。让他称为10*10的数组,或者是40*40的数组。如果设置20*40则取出的数组对于真实必威现金回扣像来说有了变形。
1
>>
>
band
.
ReadAsArray
(
100
,
100
,
10
,
10
)
2
array
(
[
[
61
,
58
,
90
,
87
,
45
,
255
,
255
,
117
,
65
,
50
]
,
3
[
36
,
59
,
113
,
88
,
37
,
26
,
63
,
165
,
23
,
74
]
,
4
[
255
,
82
,
135
,
72
,
22
,
29
,
67
,
118
,
77
,
120
]
,
5
[
45
,
129
,
144
,
255
,
255
,
36
,
94
,
108
,
88
,
97
]
,
6
[
90
,
110
,
98
,
35
,
45
,
88
,
109
,
121
,
92
,
73
]
,
7
[
94
,
108
,
85
,
45
,
55
,
97
,
144
,
167
,
135
,
21
]
,
8
[
96
,
106
,
82
,
45
,
45
,
45
,
255
,
230
,
251
,
255
]
,
9
[
246
,
236
,
255
,
255
,
3
,
255
,
255
,
247
,
254
,
255
]
,
10
[
236
,
249
,
255
,
255
,
255
,
255
,
255
,
255
,
255
,
255
]
,
11
[
194
,
212
,
255
,
255
,
255
,
255
,
255
,
255
,
255
,
255
]
]
,
'b'
)
通 过几个例子可以看到,读取的4个size参数的作用就是把硬盘上指定区域(前四个参数定义)的数据按比例缩放到用户指定区域(后两个定义)内,必要的时候 进行缩放。这里需要注意的是缩的情况(缩的时候是取几个周围点的平均值)如果是调色板数据就可能引发问题(是否可以设置重采样方式我还要再研究一下)。
补:在发现了一个帖子,发现RasterIO用的都是最临近(都是最临近?),而要设置重采样方式只能在BuildPyramid的时候设置了。现在看来BuildPyramid还是有些用的。
1.1.2. 范围
现在还有一个问题。就是取值窗口超过实际数据最大的范围怎么办?
1
>>
>
band
.
XSize
2
634
3
>>
>
band
.
YSize
4
478
5
>>
>
band
.
ReadAsArray
(
630
,
100
,
10
,
10
)
6
ERROR
5
:
Access
window
out
of
range
in
RasterIO
(
)
.
Requested
7
(
630
,
100
)
of
size
10
x10
on
raster
of
634
x478
.
8
Traceback
(
most
recent
call
last
)
:
9
File
""
,
line
1
,
in
?
10
File
"D:\Python24\lib\gdal.py"
,
line
837
,
in
ReadAsArray
11
buf_xsize
,
buf_ysize
,
buf_obj
)
12
File
"D:\Python24\lib\gdalnumeric.py"
,
line
178
,
in
BandReadAsArray
13
buf_xsize
,
buf_ysize
,
datatype
,
buf_obj
)
14
TypeError
:
Access
window
out
of
range
in
RasterIO
(
)
.
Requested
15
(
630
,
100
)
of
size
10
x10
on
raster
of
634
x478
.
16
>>
>
ERROR: EOF in multi-line statement
可以看到,出错了。获取数据的时候不能越界,很不好。还要调用的时候去判断越界没。还好用python的Numeric模块可以很方便地扩展矩阵。
1.1.3. 效率
对于gtif来说,从横向读取和重纵向读取的效率会相差十分之大(那不是一点的大)! 可以写一个python文件来验证一下
1
# gtif.py
2
import
gdal
3
import
time
4
dataset
=
gdal
.
Open
(
"f:/gisdata/gtif/zz.tif"
)
5
band
=
dataset
.
GetRasterBand
(
1
)
6
width
=
dataset
.
RasterXSize
7
height
=
dataset
.
RasterYSize
8
bw
=
128
9
bh
=
128
10
bxsize
=
width
/
128
11
bysize
=
width
/
128
12
start
=
time
.
time
(
)
13
band
.
ReadAsArray
(
0
,
0
,
width
,
height
)
14
print
time
.
time
(
)
-
start
15
start2
=
time
.
time
(
)
16
for
i
in
range
(
bysize
)
:
17
for
j
in
range
(
bxsize
)
:
18
band
.
ReadAsArray
(
bw
*
j
,
bh
*
i
,
bw
,
bh
)
19
print
time
.
time
(
)
-
start2
运行一下得到结果
F:\pyproj>gtif.py 2.35899996758 0.766000032425
然后把最后一段的循环的两个for位置调换一下(缩进要注意),变成:
1
for
j
in
range
(
bxsize
)
:
2
for
i
in
range
(
bysize
)
:
3
band
.
ReadAsArray
(
bw
*
j
,
bh
*
i
,
bw
,
bh
)
运行结果是:
F:\pyproj>gtif.py 2.48400020599 24.140999794
天!运行时间是第一次的N倍!切注意!切注意!
另外,我们可以看到,如果把必威现金回扣像分块读取,比不分块读取同样大小的必威现金回扣像效率明显要高出许多。