Python gdal チュートリアル: gdal を使用したラスターデータの読み取り-Python チュートリアル-php.cn

Python gdal チュートリアル: gdal を使用したラスターデータの読み取り

黄舟

リリース： 2016-12-24 17:07:04

オリジナル

8482 人が閲覧しました

GDAL は 100 を超えるラスターデータタイプをネイティブにサポートし、

ArcInfo グリッド、ArcSDE ラスター、Imagine、Idrisi、ENVI、GRASS、GeoTIFF

HDF4、HDF5

USGS DOQ、USGS DEM を含むすべての主流の GIS および RS データ形式をカバーします

ECW、MrSID

TIFF、JPEG、JPEG2000、PNG、GIF、BMP

完全なサポートリストについては、http://www.gdal.org/formats_list.html を参照してください

GDAL サポートライブラリのインポート

古いバージョン (1.5 より前): import gdal, gdalconst

新しいバージョン (1.6 以降): from osgeo import gdal, gdalconst

gdal と gdalconst の両方をインポートするのが最善です。gdalconst の定数は接頭辞であり、一貫性を保つように努めます。他のモジュールとの競合は最小限です。したがって、次のように gdalconst を直接インポートできます。 from osgeo.gdalconst import *

GDAL データドライバーは、OGR データドライバーと同様に、最初に特定のタイプのデータドライバーを作成してから、対応するラスターデータセットを作成する必要があります。

すべてのデータドライバーを一度に登録しますが、読み取りのみ可能で書き込みはできません: gdal.AllRegister()

読み取りと書き込みができるように、また新しいデータセットを作成できるように、特定の種類のデータドライバーを個別に登録します。

driver = gdal.GetDriverByName('HFA')

driver.Register()

既存のラスターデータセットを開きます:

fn = 'aster.img'

ds = gdal.Open(fn, GA_ReadOnly)

ds が None の場合:

print 'Could not open ' + fn

sys.exit(1)

x 方向のピクセル数、y 方向のピクセル数、およびバンドの数を読み取ります

cols = ds. Raster ラスターデータセットの地理座標情報を格納するリスト

adfGeoTransform[0] /* 左上 x x 左上隅の座標*/

adfGeoTransform[1] /* w--e ピクセル解像度東西方向のピクセル解像度 Rate*/

adfGeoTransform[2] /* 回転、画像が「北を上」の場合は 0 北が上を向いている場合、地図の回転角度*/

adfGeoTransform[3 ] /* 左上 y 左上隅の y 座標*/

adfGeoTransform[4] /* 回転、画像が「北を上」の場合は 0 北が上を向いている場合、地図の回転角度*/

adfGeoTransform [5] /* n-s ピクセル解像度南北方向のピクセル解像度*/

ラスターデータセットの座標は通常、左上隅に基づくことに注意してください。

次の例は、ラスターデータセットから Geotransform をリストとして抽出し、その中のデータを読み取ることです

geotransform = ds.GetGeoTransform()

OriginX = geotransform[0]

OriginY = geotransform[3]originY = geotransform[3]

xelWidth = geotransform[1]

pixelHeight = geotransform[5]

ある座標に対応するピクセルの相対位置 (ピクセルオフセット)、つまり座標との相対位置を計算します。左上隅のピクセル、ピクセル数によって計算されます。計算式は次のとおりです:

xOffset = int((x –originX)/pixelWidth)

yOffset = int((y –originY)/pixelHeight)

特定のピクセルの値を読み取るには、2 つのステップを分割する必要があります

最初にバンドを読み取ります: GetRasterBand()、パラメータはバンドのインデックス番号です

次に ReadAsArray(, < yoff>、、)、(xoff,yoff) から始まりサイズ (xsize,ysize) の行列を読み取ります。マトリクスサイズが 1X1 に設定されている場合、1 ピクセルが読み取られます。ただし、この方法では、たとえ 1 つのピクセルしか読み取られなかったとしても、読み取ったデータをマトリクスに配置することしかできません。例:

band = ds.GetRasterBand(1)

data =band.ReadAsArray(xOffset, yOffset, 1, 1)

画像全体を一度に読み取りたい場合は、オフセットを 0、サイズを設定します。画像全体のサイズに設定します。例:

data =band.ReadAsArray(0, 0,cols, rows)

ただし、データから特定のピクセルの値を読み取るには、data を使用する必要があることに注意してください。 [yoff、xoff]。逆にやらないように注意してください。数学の行列は [row,col] ですが、ここではその逆です。ここで、row は y 軸に対応し、col は x 軸に対応します。

バンド = なし、データセット = なしなど、適切な場合にメモリを解放することに注意してください。特に画像が非常に大きい場合

ラスターデータをより効率的に読み取るにはどうすればよいですか?明らかに、1 つずつ読み取るのは非常に非効率的です。ラスターデータセット全体を 2 次元配列に詰め込むことは、依然として大量のメモリを消費するため、得策ではありません。より良いアプローチは、ブロック内のデータにアクセスし、必要なブロックのみをメモリに入れることです。今週のサンプルコードには、ブロックサイズを読み取ることができる utils モジュールが含まれています。

例:

import utils

blockSize = utils.GetBlockSize(band)

xBlockSize = blockSize[0]

yBlockSize = blockSize[1]

タイル(タイル状)、つまりラスターデータはブロックに格納されます。 GeoTiff などの一部の形式はタイル化されず、1 行がブロックになります。 Erdas Imaging 形式は 64x64 ピクセルでタイル化されます。

行がブロックの場合、行ごとに読み取る方がリソースを節約できます。

タイル化されたデータ構造の場合、一度に 1 つのブロックのみを読み取れるように ReadAsArray() のパラメーター値を設定するのが最も効率的な方法です。例:

rows = 13、cols = 11、xBSize = 5、yBSize = 5

for i in range(0, rows, yBSize):

if i + yBSize < rows:

numRows = yBSize