1.GDAL
GDAL(Geospatial Data Abstraction Library )是非常受欢迎的、强大的栅格文件读写库。GDAL库是开源的,但是有宽松的授权,所以许多商业软件包都使用它。现有的大部分GIS或者遥感平台,不论是商业软件ArcGIS,ENVI还是开源软件GRASS,QGIS,都使用了GDAL作为底层构建库。
GDAL库以其能读写很多不同格式而被人熟知,但是也包含一些数据处理功能,如邻近分析。NumPy模块是专门处理大矩阵数据,并且可以使用GDAL直接将数据读进NumPy数组。在操作数据之后,可以使用NumPy或其他模块来处理这些矩阵,可以将数组写回磁盘作为一个栅格数据集。
可以查看GDAL光栅驱动:
每个驱动程序读写特定的数据格式。GDAL的数据结构大体上和栅格数据集匹配,每个数据集包含一个或多个波段,每个波段又包含像素数据和可能的概视图。地理参考信息包含在数据集中,因为所有的波段都使用相同的地理参考信息。
美国陆地卫星(LANDSAT)系列卫星由美国航空航天局(NASA)和美国地质调查局(USGS)共同管理。自1972年起,LANDSAT系列卫星陆续发射8颗(第6颗发射失败)。Landsat1-5目前均已退役,landsat7则自2003年6月以来,该传感器的扫描线校正器(SLC)出现故障导致的数据间隙数据,但目前仍在运行。Landsat8 于2013年2月11日发射,经过100天测试运行后开始获取影像,目前仍在运。
美国地质调查局将Landsat影像处理成 GeoTIFFs 影像集。除了波段6(热红外)和波段8(全色),每一个波段都是30米的分辨率,因为来自于同一个Landsat场景,具有相同的尺寸。波段可以直接叠置在另一个上面,而不需要修改。
1.将独立的栅格波段合成一张图像:
# Script to stack three 1-band rasters into a a single 3-band image.
import os
from osgeo import gdal
os.chdir(r'E:\gis with python\Landsat\Washington')
band1_fn = 'p047r027_7t20000730_z10_nn10.tif'
band2_fn = 'p047r027_7t20000730_z10_nn20.tif'
band3_fn = 'p047r027_7t20000730_z10_nn30.tif'
# 打开GeoTIFF波段1,索引为1
# 波段索引从1开始,而不是0
in_ds = gdal.Open(band1_fn)
in_band = in_ds.GetRasterBand(1)
# 利用和波段1相同的属性创建三波段GeoTIFF
gtiff_driver = gdal.GetDriverByName('GTiff') # 使用驱动对象创建新的数据集
out_ds = gtiff_driver.Create('nat_color.tif',
in_band.XSize, in_band.YSize, 3, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
# 从输入波段复制像素到输出波段3
in_data = in_band.ReadAsArray()
out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)
# 从数据集,而不是波段复制像素
in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())
# 为每个波段计算统计值
out_ds.GetRasterBand(1).WriteArray(
gdal.Open(band3_fn).ReadAsArray())
out_ds.FlushCache()
for i in range(1, 4):
out_ds.GetRasterBand(i).ComputeStatistics(False)
# 创建概视图或金字塔图层
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
del out_ds
RGB图像:
driver.Create(filename, xsize, ysize, [bands], [data_type], [options]):
- filename:创建的数据集的路径
- ysize:新数据集的列数
- bands:新数据集里的波段数,默认值为1
- data_type:将要储存在新数据集中的数据类型,默认为GDT_Byte
- option:创建选项字符串列表。值是基于所创建的数据集类型。
GDAL数据类型:
常量 | 数据类型 |
---|---|
GDT_Unknown | Unknown |
GDT_Byte | Unsigned 8-bit integer (byte) |
GDT_UInt16 | Unsigned 16-bit integer |
GDT_Int16 | Signed 16-bit integer |
GDT_UInt32 | Unsigned 32-bit integer |
GDT_Int32 | Signed 32-bit integer |
GDT_Float32 | 32-bit floating point |
GDT_Float64 | 64-bit floating point |
GDT_CInt16 | 16-bit complex integer |
GDT_CInt32 | 32-bit complex integer |
GDT_CFloat32 | 32-bit complex floating point |
GDT_CFloat64 | 64-bit complex floating point |
GDT_TypeCount | Number of available data types |
从输入波段中获取此信息,但是这些图像使用了GDT_Byte默认类型。
1. 查看空的三波段数据集的SRS:
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
从输入数据集中得到SRS并复制到新数据集,再对 geotransform 做同样的操作。geotransform 提供原始坐标和像素大小,并伴随着旋转值,可使图像顶部朝北。如果需要将数据集放置到正确的空间位置,数据集和像素大小是非常重要的参数。
2. 创建数据集后,就要添加像素值,将读取的波段1的像素值读进 NumPy 数组。不给 Readasarray() 函数提供任何参数,所有的像素值就会返回在一个和栅格具有相同尺寸的二维数组里。in_data变量保存了一个二维数组像素值:
in_data = in_band.ReadAsArray()
3. 因为波段1是蓝色波段,所以必须放到输出图层的第3个波段位置上,从而得到RGB波段顺序。
然后,从out_data获取第3个额伯顿啊,并用WriteArray()函数将in_data数组里的像素值复制到新数据集的第3波段。
out_band = out_ds.GetRasterBand(3)
out_band.WriteArray(in_data)
4. 将红色和绿色波段添加到数据集中,并打开第2波段的GeoTIFF。直接从数据集中读取像素数据。
in_ds = gdal.Open(band2_fn)
out_band = out_ds.GetRasterBand(2)
out_band.WriteArray(in_ds.ReadAsArray())
当数据集上调用ReadAsArray函数时,如果正在读取的数据集有多个波段,将会得到一个三维数组。因为Landsat文件只有一个波段,所以返回的是二维数组。
5. 对红色像素值做同样的事情,但是压缩为更少的代码:
out_ds.GetRasterBand(1).WriteArray(gdal.Open(band3_fn).ReadAsArray())
6. 计算数据集中每个波段的统计数据,能够使软件更容易显示。统计数据包括波段的平均值、最小值、最大值和标准查。使得GIS软件可以使用这些信息来拉伸屏幕数据,使得看起来更好。
统计数据之前,必须确保数据已经写入磁盘,计算每个波段的统计数据。传递False给这个函数,告诉它你需要的实际的统计数据,而不是估计值,因为他可能从概视图中或采样像素的子集中得到。若估计值可以接受,就传递True给函数;使得计算速度更快,因为不是每个像素都需要检查:
out_ds.FlushCache()
for i in range(1, 4):
out_ds.GetRasterBand(i).ComputeStatistics(False)
7. 建立数据集的概视图
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
可以在rasterio上查看,此模块依赖于GDAL,实质上是使用GDAL来读写数据,但rasterio模块尝试使栅格数据的处理工作更加轻松。
imageio模块不依赖于GDAL,是纯用python编写的模块,不专注于地理空间数据,可以读取和写入许多不同的栅格格式,包括视频格式。
2.读取部分数据集
访问子集,而不是访问整个图像。ReadAsAsrray()函数有几个可选参数,会根据你是否在用数据集或波段而变化。
band.ReadAsArray([xoff], [yoff], [win_xsize], [win_ysize], [buf_xsize], [buf_ysize], [buf_obj])
- xoff:列读取的起点,默认值为0
- yoff:行读取的起点,默认值为0
- win_xsize:要读取的列数,默认读取所有列
- win_ysize:要读取的行数,默认读取所有行
- buf_xsize:输出数组里的列数,默认为win_xsize的值,如果值不同于win_xsize,数据将会重采样
- buf_ysize:输出数组里的行数,默认为win_xsize的值,如果值不同于win_ysize,数据将会重采样
- buf_obj:用于存放数组,而不是创建数组的NumPy数组。如果需要,数据将会重采样来适应数组。值将会转换为这种数组的类型。
xoff和yoff参数分别用于指定开始读取的列和行的偏移值。默认从第一行和第一列开始。
win_xsize和win_ysize参数用来显示读取了多少行和列,默认读取所有行和列
如果提供的buf_xsize和buf_ysize和该数组的大小不匹配,就会得到一个错误,但是没必要给出buf_xsize和buf_ysize的值,因为这用数组本身确定。
1. 读取从第6000行第1400列开始的3行和6列数据:
import os
import numpy as np
from osgeo import gdal
data_dir = r'E:\gis with python'
# 打开Landsat波段
os.chdir(os.path.join(data_dir, 'Landsat', 'Washington'))
ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
band = ds.GetRasterBand(1)
# 读取指定范围数据
data = band.ReadAsArray(1400, 6000, 6, 3)
print(data)
2. 数据类型转换:
# 使用numpy将数据转换为浮点数
data = band.ReadAsArray(1400, 6000, 6, 3).astype(float)
print(data)
3. 读取数据时,用GDAL来转换:
首先,创建一个浮点数组,然后把它作为一个buf_obj参数传递给readasarray,同时确保创建的数组于正在读取的数据大小相同。
# 通过将它们读入浮点数组,来转换为浮点数
import numpy as np
data = np.empty((3, 6), dtype=float) # dtype参数用来指定数组的数据类型
band.ReadAsArray(1400, 6000, 6, 3, buf_obj=data)
print(data)
4. 将数据数组写入其他数据集的指定位置,给函数传递偏移值。它将从提供的偏移值开始,写入传递给 WriteArray() 函数的数组里的所有数据:
band2.WriteArray(data, 1400, 6000)
5.处理一个超过内存大小的大数据集: 一次只处理一块。
注意: 栅格数据在磁盘上按块储存,所以在这些块中处理图像效率非常高。如图所示:
将DEM的高程从米制转换成英尺(一次只处理一块):
import os
import numpy as np
from osgeo import gdal
os.chdir(r'E:\gis with python\Washington\dem')
# 打开输入光栅并获得其尺寸
in_ds = gdal.Open('gt30w140n90.tif')
in_band = in_ds.GetRasterBand(1)
xsize = in_band.XSize
ysize = in_band.YSize
# 获取块大小和NoData值
block_xsize, block_ysize = in_band.GetBlockSize()
nodata = in_band.GetNoDataValue()
# 创建具有相同维度和数据类型的输出文件
out_ds = in_ds.GetDriver().Create(
'dem_feet.tif', xsize, ysize, 1, in_band.DataType)
out_ds.SetProjection(in_ds.GetProjection())
out_ds.SetGeoTransform(in_ds.GetGeoTransform())
out_band = out_ds.GetRasterBand(1)
# 在x方向循环这些块
for x in range(0, xsize, block_xsize):
# 获取要读取的列数
if x + block_xsize < xsize:
cols = block_xsize
else:
cols = xsize - x
# 在y方向循环blocks
for y in range(0, ysize, block_ysize):
# 获取要读取的行数
if y + block_ysize < ysize:
rows = block_ysize
else:
rows = ysize - y
# 读取一个块的数据值,将其转换为英尺
# 然后将结果写入输出中相同的块位置
data = in_band.ReadAsArray(x, y, cols, rows)
data = np.where(data == nodata, nodata, data * 3.28084)
out_band.WriteArray(data, x, y)
# 刷新缓存
# 设置NoData值后进行统计
out_band.FlushCache()
out_band.SetNoDataValue(nodata)
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32])
del out_ds
这种方法比一次性读取整个波段来说更加复杂,但是节约内存。
3.现实世界坐标
前提:真实世界坐标和栅格数据坐标使用相同的SRS
需要的数据:坐标原点、像素大小(单个像元所代表的实际地物大小)和旋转值,都储存在从数据集间复制过来的 geotransform(地理变换)里。地理变换是一个元组,包含六个值:
索引 | 描述 |
---|---|
0 | 原点的x坐标 |
1 | 像素宽度 |
2 | x像素旋转(图像朝北为0°)(默认为0) |
3 | 原点的y坐标 |
4 | y像素旋转(图像0°为朝北)(默认为0) |
5 | 像素高度(负值) |
GDAL提供的ApplyGeoTransform函数来进行仿射变换,需要geotransform、x值和y值3个参数。当使用一个数据集的geotransform变换时,这个函数能将图像坐标(偏移)转换为现实世界坐标。
GDAL提供了一个ApplyGeoTransform函数进行仿射变换,parameter:geotransform、x值、y值。
另一种思路(数据集的逆变换):
1. GDAL 1.X版本: InvGeoTransform函数
# 现在得到逆变换
# 原始文件可用于将偏移量转换为真实世界的坐标
# 而逆文件可用于将真实世界的坐标转换为偏移量
import os
import numpy as np
from osgeo import gdal
data_dir = r'D:\geodata'
#从陆地卫星的一个波段中获取geotransform
os.chdir(os.path.join(data_dir, 'Landsat', 'Washington'))
ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
band = ds.GetRasterBand(1)
gt = ds.GetGeoTransform()
print(gt)
# GDAL 1.x: 得到一个成功标志和地理变换
success, inv_gt = gdal.InvGeoTransform(gt)
print(success, inv_gt)
成功返回1,不成功返回0。 0:表示仿射变换不可逆。
2. GDAL 2.X版本(我的版本为3.2.1): InvGeoTransform函数,使用此函数将现实世界坐标转换为图像坐标。
inv_gt = gdal.InvGeoTransform(gt)
print(inv_gt)
(-12060.5, 0.03508771929824561, 0.0, 188406.5, 0.0, -0.03508771929824561)
对照上面的表格进行解释。
3. 坐标转换实例 1: 获取(465200, 5296000)处的像素值
# 使用inverset geotransform从真实世界的UTM坐标获得像素偏移量
# 输出为浮点型
offsets = gdal.ApplyGeoTransform(inv_gt, 465200, 5296000)
print(offsets)
# 将偏移量转换为整型
xoff, yoff = map(int, offsets)
print(xoff, yoff)
[4262.307017543859, 2581.9385964912362]
4262 2581
# 读取像素值
# ReadAsArray函数需要使用整型
# ReadAsArray函数返回一个二维数组
value = band.ReadAsArray(xoff, yoff, 1, 1)[0,0] # [0,0]表示数组里的第一行,第一列的值
print(value)
最终,(465200, 5296000)处的像素值为:
62
存在的问题: 单个像素值的提取效率低,不适用于整个波段的像素值。
改进的方法:
# 先读取所有像素值
# 再取出所需要的像素值
data = band.ReadAsArray()
x, y = map(int, gdal.ApplyGeoTransform(inv_gt, 465200, 5296000))
value = data[yoff, xoff] # 这里需要注意,这里使用的是Numpy读取的
print(value)
62
注意: 使用Numpy数组的[行,列]偏移与GDAL使用的方式相反。
4. 实例 2: 提取图像空间子集,并保存(感兴趣区域已知)。即将图像的现实世界坐标转换为像素偏移值。
数据说明: 一幅真彩色Landsat图像(Landsat_color.tif)
import os
from osgeo import gdal
# 提取的边框的坐标
interested_ulx, interested_uly = 517296, 5296820
interested_lrx, interested_lry = 540315, 5267134
os.chdir(r'D:\geodata\Landsat\Washington')
in_ds = gdal.Open('Landsat_color.tif')
# 将现实世界的坐标转换为像素偏移量
in_gt = in_ds.GetGeoTransform()
inv_gt = gdal.InvGeoTransform(in_gt)
if gdal.VersionInfo()[0] == '1':
if inv_gt[0] == 1:
inv_gt = inv_gt[1]
else:
raise RuntimeError('Inverse geotransform failed')
elif inv_gt is None:
raise RuntimeError('Inverse geotransform failed')
# 获取与边框角坐标对应的偏移量
# 计算左上角和右下角偏移
offsets_ul = gdal.ApplyGeoTransform(
inv_gt, interested_ulx, interested_uly)
offsets_lr = gdal.ApplyGeoTransform(
inv_gt, interested_lrx, interested_lry)
off_ulx, off_uly = map(int, offsets_ul)
off_lrx, off_lry = map(int, offsets_lr)
# 计算需要提取的行列数
rows = off_lry - off_uly
columns = off_lrx - off_ulx
# 创建具有正确行列数的输出光栅
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('interested.tif', columns, rows, 3)
out_ds.SetProjection(in_ds.GetProjection())
# 将新的原点坐标放入地理变换中
subset_ulx, subset_uly = gdal.ApplyGeoTransform(
in_gt, off_ulx, off_uly)
out_gt = list(in_gt)
out_gt[0] = subset_ulx
out_gt[3] = subset_uly
out_ds.SetGeoTransform(out_gt)
# 遍历三个波段
for i in range(1, 4):
in_band = in_ds.GetRasterBand(i)
out_band = out_ds.GetRasterBand(i)
# 利用计算值读入数据
data = in_band.ReadAsArray(
off_ulx, off_uly, columns, rows)
# 从原点处写入数据
out_band.WriteArray(data)
del out_ds
结果展示:
![](https://i-blog.csdnimg.cn/blog_migrate/34a341b1908e1c2737906aa3052bf962.png)
![](https://i-blog.csdnimg.cn/blog_migrate/0450eb37f18bc4444af33fb02135cdc5.png)
图像属性信息:
![](https://i-blog.csdnimg.cn/blog_migrate/042947cd2a91fb6f48c9c0e52ab0dc7a.png)
![](https://i-blog.csdnimg.cn/blog_migrate/cc280f20c4c6987adb6d9160b7266a91.png)