Python gdal チュートリアル: マップ代数とラスター データの書き込み

NDVI の計算を例に挙げます:

NDVI=(NIR-RED)/(NIR+RED)

ここで、NIR はバンド 3、RED はバンド 2 です

プログラミングのポイントは次のとおりです:

1.バンド 3 を配列 data3 に読み取り、バンド 2 を配列 data2

2 に読み取ります。計算式は次のとおりです: data3 と data2 が両方とも 0 (たとえば、0 は NODATA を表します) の場合、0 による除算エラーが発生します。が発生し、プログラムがクラッシュします。 0の値を削除するには、マスクを選択して使用する必要があります

コードは次のとおりです。行は4行だけです

data2 = Band2.ReadAsArray(0, 0,cols,rows).astype(Numeric.Float16)

data3 =band3.ReadAsArray(0,0,cols,rows).astype(Numeric.Float16)

mask = Numeric.greater(data3 + data2, 0)

ndvi = Numeric.choose(mask, (-99, ( data3 - data2) / (data3 + data2)))

新しいラスター データセットを作成します

計算したばかりのデータを新しいラスター データセットに書き込みます

まず、データ ドライバーをコピーします:

driver = inDataset.GetDriver()

後で新しいデータセットを作成します

Create(, , , [], [])

バンドのデフォルト値は 1 です。および GDALDataType のデフォルト値です。デフォルトのタイプは GDT_Byte です。例:

outDataset = driver.Create(filename,cols,rows, 1, GDT_Float32)

このステートメントの実行中、ストレージ領域はハードウェアに割り当てられています。ディスク

を書き込む前に、最初にバンド オブジェクトを導入する必要もあります

outBand = outDataset.GetRasterBand(1)

バンド オブジェクトは行列への直接書き込みをサポートしており、2 つのパラメーターは x 方向のオフセットと y 方向のオフセットです

outBand.WriteArray(ndvi, 0 , 0)

以下の例は、今回と前回のブロック単位の書き込み方法をまとめたものです

xBlockSize = 64

yBlockSize = 64

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

if i + yBlockSize

numRows = yBlockSize

else:

numRows = rowsnumRows = rows –– ii

for j in range(0,cols, xBlockSize) :

if j + xBlockSize numCols = xBlockSize

outBand.WriteArray(outData, j, i)

NoData 値

outBand.SetNoDataValue(-99)

次のこともできますNoData 値を読み取ります

ND = outBand.GetNoDataValue()

バンドの統計を計算します

最初に FlushCache( ) を使用します

キャッシュされたデータをディスクに書き込みます

それから GetStatistics(, ) を使用します統計を計算します。 around_ok=1 の場合、計算はピラミッドに基づいて行われ、force=0 の場合、グラフ全体を再読み取りする必要がある場合、統計は計算されません。

outBand.FlushCache()

outBand.GetStatistics(0, 1)

新しい画像の地理的参照点を設定します

新しい画像が別の画像と全く同じ地理的参照情報を持っている場合、それは非常に簡単です

geoTransform = inDataset.GetGeoTransform()

outDataset.SetGeoTransform(geoTransform)

proj = inDataset.GetProjection()

outDataset.SetProjection(proj)

ピラミッドを作成する

Imagineスタイルのピラミッドを設定する

gdal.SetConfigオプション(' HFA_USE_RRD', 'YES')

ピラミッドの構築を強制

outDataset.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])

画像のステッチ

1.行数と列数、原点（minX、maxY）、ピクセル長さ、ピクセル幅を読み取り、座標範囲を計算します

maxX1 = minX1 + (cols1 *PixelWidth)

minY1 = maxY1 + (rows1 *xelHeight)

2 . 出力画像の座標範囲を計算します:

minX = min(minX1, minX2, …) maxX = max(maxX1, maxX2, …)

minY = min( minY1, minY2, …) maxY = max(maxY1, maxY2, …)

3. 出力画像の行数と列数を計算します:

cols = int((maxX – minX) /xelWidth)

rows = int((maxY – minY) / abs(pixelHeight)

4. 出力画像を作成して初期化します

5. ステッチする画像ごとに: オフセット値を計算します

xOffset1 = int((minX1 - minX) /xelWidth )

yOffset1 = int((maxY1 - maxY) /xelHeight)

データを読み取り、上記で計算されたオフセットに従って書き込みます

6. 出力画像の場合: 統計を計算し、geotransform を設定します: [minX, ピクセル幅, 0] 、maxY、0、pixelHeight]、投影を設定し、ピラミッドを作成します

上記は Python gdal チュートリアルの内容です: マップ代数とラスター データの書き込み その他の関連コンテンツについては、PHP 中国語 Web サイト (www.php.cn) に注目してください。