使用Rasterio讀取柵格數(shù)據(jù)的實(shí)例講解

更新時(shí)間：2019年11月26日 10:20:38 作者：TheOneGIS

今天小編就為大家分享一篇使用Rasterio讀取柵格數(shù)據(jù)的實(shí)例講解，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧

Rasterio簡(jiǎn)介

有沒(méi)有覺(jué)得用GDAL的Python綁定書(shū)寫(xiě)的代碼很不Pythonic，強(qiáng)迫癥的你可能有些忍受不了。不過(guò)，沒(méi)關(guān)系，MapBox旗下的開(kāi)源庫(kù)Rasterio幫我們解決了這個(gè)痛點(diǎn)。

Rasterio是基于GDAL庫(kù)二次封裝的更加符合Python風(fēng)格的主要用于空間柵格數(shù)據(jù)處理的Python庫(kù)。

Rasterio中柵格數(shù)據(jù)模型基本和GDAL類(lèi)似，需要注意的是：

在Rasterio 1.0以后，對(duì)于GeoTransform的表示棄用了GDAL風(fēng)格的放射變換，而使用了Python放射變換的第三方庫(kù)affine庫(kù)的風(fēng)格。

對(duì)于放射變換

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

GDAL中對(duì)應(yīng)的參數(shù)順序是：(c, a, b, f, d, e)

采用新的放射變換模型的好處是，如果你需要計(jì)算某個(gè)行列號(hào)的地理坐標(biāo)，直接使用行列號(hào)跟給放射變換對(duì)象相乘即可，完全符合數(shù)學(xué)上矩陣乘法的操作，更加直觀和方便。

柵格數(shù)據(jù)讀取代碼示例

下面的示例程序中演示了如何讀取一個(gè)GeoTIFF文件并獲取相關(guān)信息，需要注意的是：

1、rasterio使用rasterio.open()函數(shù)打開(kāi)一個(gè)柵格文件

2、rasterio使用read()函數(shù)可以將數(shù)據(jù)集轉(zhuǎn)為numpy.ndarray，該函數(shù)如果不帶參數(shù)，將把數(shù)據(jù)的所有波段做轉(zhuǎn)換（第一維是波段數(shù)），如果指定波段，則只取得指定波段對(duì)應(yīng)的數(shù)據(jù)（波段索引從1開(kāi)始）

3、數(shù)據(jù)的很多元信息都是以數(shù)據(jù)集的屬性進(jìn)行表示的

import rasterio

with rasterio.open('example.tif') as ds:
 print('該柵格數(shù)據(jù)的基本數(shù)據(jù)集信息（這些信息都是以數(shù)據(jù)集屬性的形式表示的）：')
 print(f'數(shù)據(jù)格式：{ds.driver}')
 print(f'波段數(shù)目：{ds.count}')
 print(f'影像寬度：{ds.width}')
 print(f'影像高度：{ds.height}')
 print(f'地理范圍：{ds.bounds}')
 print(f'反射變換參數(shù)（六參數(shù)模型）：\n {ds.transform}')
 print(f'投影定義：{ds.crs}')
 # 獲取第一個(gè)波段數(shù)據(jù)，跟GDAL一樣索引從1開(kāi)始
 # 直接獲得numpy.ndarray類(lèi)型的二維數(shù)組表示，如果read()函數(shù)不加參數(shù)，則得到所有波段（第一個(gè)維度是波段）
 band1 = ds.read(1)
 print(f'第一波段的最大值：{band1.max()}')
 print(f'第一波段的最小值：{band1.min()}')
 print(f'第一波段的平均值：{band1.mean()}')
 # 根據(jù)地理坐標(biāo)得到行列號(hào)
 x, y = (ds.bounds.left + 300, ds.bounds.top - 300) # 距離左上角東300米，南300米的投影坐標(biāo)
 row, col = ds.index(x, y) # 對(duì)應(yīng)的行列號(hào)
 print(f'(投影坐標(biāo){x}, {y})對(duì)應(yīng)的行列號(hào)是({row}, {col})')
 # 根據(jù)行列號(hào)得到地理坐標(biāo)
 x, y = ds.xy(row, col) # 中心點(diǎn)的坐標(biāo)
 print(f'行列號(hào)({row}, {col})對(duì)應(yīng)的中心投影坐標(biāo)是({x}, {y})')
 # 那么如何得到對(duì)應(yīng)點(diǎn)左上角的信息
 x, y = (row, col) * ds.transform
 print(f'行列號(hào)({row}, {col})對(duì)應(yīng)的左上角投影坐標(biāo)是({x}, {y})')

輸出如下：

該柵格數(shù)據(jù)的基本數(shù)據(jù)集信息（這些信息都是以數(shù)據(jù)集屬性的形式表示的）：
數(shù)據(jù)格式：GTiff
波段數(shù)目：3
影像寬度：4800
影像高度：4800
地理范圍：BoundingBox(left=725385.0, bottom=2648415.0, right=869385.0, top=2792415.0)
反射變換參數(shù)（六參數(shù)模型）：
 | 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
(投影坐標(biāo)725685.0, 2792115.0)對(duì)應(yīng)的行列號(hào)是(10, 10)
行列號(hào)(10, 10)對(duì)應(yīng)的中心投影坐標(biāo)是(725700.0, 2792100.0)
行列號(hào)(10, 10)對(duì)應(yīng)的左上角投影坐標(biāo)是(725685.0, 2792115.0)

以上這篇使用Rasterio讀取柵格數(shù)據(jù)的實(shí)例講解就是小編分享給大家的全部?jī)?nèi)容了，希望能給大家一個(gè)參考，也希望大家多多支持腳本之家。

您可能感興趣的文章: