前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >打开栅格数据的正确方式

打开栅格数据的正确方式

作者头像
卡尔曼和玻尔兹曼谁曼
发布2019-01-22 09:41:34
1.1K0
发布2019-01-22 09:41:34
举报
文章被收录于专栏:给永远比拿愉快

以一个简单例子说明如何打开栅格影像

下面的例子打开一副GeoTIFF影像,输出了影像的一些信息,然后遍历了所有波段,输出波段的一些信息

代码语言:javascript
复制
import gdal

# 打开栅格数据集
ds = gdal.Open('example.tif')

# 获得栅格数据的一些重要信息
print(f'投影信息:{ds.GetProjection()}')
print(f'栅格波段数:{ds.RasterCount}')
print(f'栅格列数(宽度):{ds.RasterXSize}')
print(f'栅格行数(高度):{ds.RasterYSize}')

# 获取数据集的元数据信息
metadata = ds.GetMetadata_Dict()
for key, value in metadata.items():
    print(f'{key} -> {value}')


for b in range(ds.RasterCount):
    # 注意GDAL中的band计数是从1开始的
    band = ds.GetRasterBand(b + 1)
    # 波段数据的一些信息
    print(f'数据类型:{gdal.GetDataTypeName(band.DataType)}')  # DataType属性返回的是数字
    print(f'NoData值:{band.GetNoDataValue()}')  # 很多影像都是NoData,我们在做数据处理时要特别对待
    print(f'统计值(最大值最小值):{band.ComputeRasterMinMax()}')  # 有些数据本身就存储了统计信息,有些数据没有需要计算

# 关闭数据集
ds = None

输出如下:

代码语言:javascript
复制
投影信息:PROJCS["WGS 84 / UTM zone 49N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",111],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32649"]]
栅格波段数:3
栅格列数(宽度):4800
栅格行数(高度):4800
AREA_OR_POINT -> Area
数据类型:Int16
NoData值:-28672.0
统计值(最大值最小值):(-435.0, 6134.0)
数据类型:Int16
NoData值:-28672.0
统计值(最大值最小值):(-468.0, 6265.0)
数据类型:Int16
NoData值:-28672.0
统计值(最大值最小值):(21.0, 7267.0)

如何将Dataset转为Numpy的ndarray

当我们得到Band对象以后,如果按照GDAL的C/C++接口惯例,我们可以使用WriteRaster()方法进行数据写入(C/C++接口是WriteBlock()),但是在Python中我们有很强大的ndarray对象,所以我们一般是将Band对象中存储的数据转为ndarray进行处理以后,然后再写回去。

下面介绍几种转换的方法:

  1. Dataset级别进行转换,转换结果是一个三维数组,第一个维度是波段数
  2. Band级别进行转换,转换的结果是一个二维数据
  3. 使用gdal_array模块中的LoadFile()函数直接进行(相当于第一种转换)
代码语言:javascript
复制
import gdal

# 打开栅格数据集
ds = gdal.Open('example.tif')
# 在数据集层面转换
image = ds.ReadAsArray()

print(f'数据的尺寸:{image.shape}')
# 输出结果为:数据的尺寸:(3, 4800, 4800)
# 这说明ReadAsArray方法将每个波段都转换为了一个二维数组

# 获得第一个波段的数据
band1 = image[0]

# 在波段层面的转换
for b in range(ds.RasterCount):
    # 注意GDAL中的band计数是从1开始的
    band = ds.GetRasterBand(b + 1)
    band = band.ReadAsArray()
    print(f'波段大小:{band.shape}')

# 关闭数据集
ds = None

输出结果:

代码语言:javascript
复制
数据的尺寸:(3, 4800, 4800)
波段大小:(4800, 4800)
波段大小:(4800, 4800)
波段大小:(4800, 4800)

使用gdal_array模块

代码语言:javascript
复制
from osgeo import gdal_array
# gdal_array模块
image = gdal_array.LoadFile('example.tif')
print(f'数据的尺寸:{image.shape}')

在GDAL中使用Python的异常对象

代码语言:javascript
复制
import gdal
import sys

# 允许GDAL跑出Python异常
gdal.UseExceptions()

try:
    ds = gdal.Open('example.tif')
except (FileNotFoundError, RuntimeError) as e:
    print('文件打开失败!')
    print(e)
    sys.exit(1)
本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2018年05月17日,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 以一个简单例子说明如何打开栅格影像
  • 如何将Dataset转为Numpy的ndarray
  • 在GDAL中使用Python的异常对象
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档