专栏首页给永远比拿愉快使用Rasterio读取栅格数据

使用Rasterio读取栅格数据

Rasterio简介

有没有觉得用GDAL的Python绑定书写的代码很不Pythonic,强迫症的你可能有些忍受不了。不过,没关系,MapBox旗下的开源库Rasterio帮我们解决了这个痛点。

Rasterio是基于GDAL库二次封装的更加符合Python风格的主要用于空间栅格数据处理的Python库。

Rasterio中栅格数据模型基本和GDAL类似,需要注意的是:

在Rasterio 1.0以后,对于GeoTransform的表示弃用了GDAL风格的放射变换,而使用了Python放射变换的第三方库affine库的风格。

对于放射变换

affine.Affine(a, b, c,
              d, e, f)

GDAL中对应的参数顺序是:(c, a, b, f, d, e)

采用新的放射变换模型的好处是,如果你需要计算某个行列号的地理坐标,直接使用行列号跟给放射变换对象相乘即可,完全符合数学上矩阵乘法的操作,更加直观和方便。

栅格数据读取代码示例

下面的示例程序中演示了如何读取一个GeoTIFF文件并获取相关信息,需要注意的是:

  1. rasterio使用rasterio.open()函数打开一个栅格文件
  2. rasterio使用read()函数可以将数据集转为numpy.ndarray,该函数如果不带参数,将把数据的所有波段做转换(第一维是波段数),如果指定波段,则只取得指定波段对应的数据(波段索引从1开始)
  3. 数据的很多元信息都是以数据集的属性进行表示的
import rasterio

with rasterio.open('example.tif') as ds:
    print('该栅格数据的基本数据集信息(这些信息都是以数据集属性的形式表示的):')
    print(f'数据格式:{ds.driver}')
    print(f'波段数目:{ds.count}')
    print(f'影像宽度:{ds.width}')
    print(f'影像高度:{ds.height}')
    print(f'地理范围:{ds.bounds}')
    print(f'反射变换参数(六参数模型):\n {ds.transform}')
    print(f'投影定义:{ds.crs}')
    # 获取第一个波段数据,跟GDAL一样索引从1开始
    # 直接获得numpy.ndarray类型的二维数组表示,如果read()函数不加参数,则得到所有波段(第一个维度是波段)
    band1 = ds.read(1)
    print(f'第一波段的最大值:{band1.max()}')
    print(f'第一波段的最小值:{band1.min()}')
    print(f'第一波段的平均值:{band1.mean()}')
    # 根据地理坐标得到行列号
    x, y = (ds.bounds.left + 300, ds.bounds.top - 300)  # 距离左上角东300米,南300米的投影坐标
    row, col = ds.index(x, y)  # 对应的行列号
    print(f'(投影坐标{x}, {y})对应的行列号是({row}, {col})')
    # 根据行列号得到地理坐标
    x, y = ds.xy(row, col)  # 中心点的坐标
    print(f'行列号({row}, {col})对应的中心投影坐标是({x}, {y})')
    # 那么如何得到对应点左上角的信息
    x, y = (row, col) * ds.transform
    print(f'行列号({row}, {col})对应的左上角投影坐标是({x}, {y})')

输出如下:

该栅格数据的基本数据集信息(这些信息都是以数据集属性的形式表示的):
数据格式:GTiff
波段数目:3
影像宽度:4800
影像高度:4800
地理范围:BoundingBox(left=725385.0, bottom=2648415.0, right=869385.0, top=2792415.0)
反射变换参数(六参数模型):
 | 30.00, 0.00, 725385.00|
| 0.00,-30.00, 2792415.00|
| 0.00, 0.00, 1.00|
投影定义:CRS({'init': 'epsg:32649'})
第一波段的最大值:5459
第一波段的最小值:-313
第一波段的平均值:489.80300625
(投影坐标725685.0, 2792115.0)对应的行列号是(10, 10)
行列号(10, 10)对应的中心投影坐标是(725700.0, 2792100.0)
行列号(10, 10)对应的左上角投影坐标是(725685.0, 2792115.0)

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 地图投影

    我们的地球是圆的,而我们的纸张是平面。为了将地球绘制在平面纸张上,我们需要将地球表面投影到平面上。地图投影的实质是建立空间地理坐标和平面直角坐标关系的过程。

    卡尔曼和玻尔兹曼谁曼
  • NumPy中的维度Axis

    NumPy中的维度是一个很重要的概念,很多函数的参数都需要给定维度Axis,如何直观的理解维度呢?我们首先以二维数组为例进行说明,然后推广到多维数组。

    卡尔曼和玻尔兹曼谁曼
  • NumPy中的维度Axis

    NumPy中的维度是一个很重要的概念,很多函数的参数都需要给定维度Axis,如何直观的理解维度呢?我们首先以二维数组为例进行说明,然后推广到多维数组。

    卡尔曼和玻尔兹曼谁曼
  • Numpy统计计算、数组比较,看这篇就够了

    如上述例子所示,axis = 1计算的是行的和,结果以列的形式展示。axis = 0计算的是列的和,结果以行的形式展示。

    华章科技
  • 用MobX管理状态(ES5实例描述)-1.核心概念和基本流程

    江米小枣
  • 学好 Python 的 11 个优秀资源

    Python是目前最流行、最易学最强大的编程语言之一(学习Python的五大理由),无论你是新手还是老鸟,无论是用于机器学习还是web开发(Pinterest就...

    小莹莹
  • 【LEETCODE】模拟面试-46. Permutations

    notice: if ( curList.contains(arr[i]) ){ continue; ...

    杨熹
  • PHP中echo与print和print_r

    echo是PHP语句, print和print_r是函数,语句没有返回值,函数可以有返回值(即便没有用)

    用户7657330
  • Python 打印语句

    首先申明下,本文为笔者学习《Python学习手册》的笔记,并加入笔者自己的理解和归纳总结。

    py3study
  • 大数据周周看 | 谷歌开疆拓土增建数据中心 科达股份斥资32.7亿元加码汽车大数据

    <数据猿导读> 随着“大数据”在各行各业的不断深入,无论是政府还是企业都争相使出浑身解数在大数据领域开疆拓土。下面就请随小编一起看看在本周大数据领域又发生了什么...

    数据猿

扫码关注云+社区

领取腾讯云代金券