Использование библиотеки растров Python Rasterio для создания подмножества изображения TIFF, а затем его отображения и сохранения? - PullRequest
0 голосов
/ 18 ноября 2018

У меня есть два растровых изображения, одно из группы 4 с B4 в конце и другое из группы 5 с B5 в конце.Я хочу настроить растр B5 на 800x600, затем отобразить его и сохранить как GeoTiff.Затем я хочу вычислить NDVI (я предполагаю, что для этого мне понадобятся и B4, и B5, но я не уверен).Затем я хочу отобразить подмножество NDVI растра B5.Отобразите его и сохраните как GeoTiff.

Как бы я создал что-то вроде подмножества растрового изображения TIFF размером 800 x 600 пикселей?Я также хочу взять этот TIFF и сгенерировать NDVI-образ для этого подмножества.

ПРИМЕЧАНИЕ. Я работаю с изображением Landsat.На изображении в конце заголовка файла стоит буква B5.

Что я сделал до сих пор:

import rasterio
from rasterio.windows import Window 
import matplotlib.pyplot plt # for later use

with rasterio.open('MyRasterImage.tif') as src:
    w = src.read(1, window=Window(0, 0, 800, 600))

Я хочу отобразить это, используя записные книжки Spyder или Jupyter.Поэтому я подумал использовать matplotlib и выполнил текущий код:

# Plot
plt.imshow(w)
plt.show()

При этом создается окно matplotlib размером 800x600, но оно все фиолетовое, не знаю, зачем оно это создает.

Теперь я хочучтобы иметь возможность отображать это изображение 800x600.Затем, после этого, я хочу предварительно сформировать NDVI для этого изображения подмножества 800x600.Затем отобразите подмножество 800x600 с отображением NDVI.

Я знаю, что формула: NDVI = (NIR - красный) / (NIR + красный)

Но как мне извлечь NIR и красныйвот из этого единственного изображения Landsat?

Моя попытка сделать это:

band1 = dataset.read(1)
band2 = dataset.read(2)
band3 = dataset.read(3)
print(band[2])

Когда я запускаю этот код для полос, я получаю ошибку:

rasterio indexerror: band index 2 out of range (not in (1,))

КогдаЯ запускаю этот код:

print(w.count)

Возвращает '1'.

Значит, это означает, что изображение Landsat имеет только одну полосу?Но для того, чтобы сделать NDVI, мне не нужны 3 полосы?

Я думаю написать некоторый код, подобный этому, чтобы получить NDVI из этого растра.Но я не уверен, что делать с извлечением полос:

# We handle the connections with "with"
with rasterio.open(bands[0]) as src:
    b3 = src.read(1)

with rasterio.open(bands[1]) as src:
    b4 = src.read(1)

# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')

# Calculate NDVI
ndvi = (b4.astype(float) - b3.astype(float)) / (b4 + b3)

Этот код не работает, потому что полосы не определены как что-либо, поэтому я не знаю, как определить полосы для получения NDVI.

После этого я не уверен, что делать с отображением и сохранением отрендеренного изображения.

...