Используя matplotlib и Rasterio, я пытаюсь сохранить растр как GeoTIFF, а также переназначить его?

Question

Я смог построить и отобразить свои растровые изображения с помощью matplotlib. Эта часть успешна. Часть, на которой я застрял, может каким-то образом сохранить этот сюжет. Для растерио я нашел два полезных урока:

https://rasterio.readthedocs.io/en/latest/topics/windowed-rw.html

а также

https://www.earthdatascience.org/courses/earth-analytics-python/multispectral-remote-sensing-in-python/export-numpy-array-to-geotiff-in-python/

Я получил вычисление для функции с именем NDVI, и с помощью matplotlib я могу отобразить ее так, как я хочу, с помощью следующего кода. Но когда я сохраняю файл в формате GeoTIFF, изображение на рабочем столе становится черным. Я также планирую перепроектировать данные, и я закомментировал этот код.

Вот мой код:

import rasterio
import matplotlib.pyplot as plt
import numpy as np


nirband = r"LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF"

redband =r"LC08_L1TP_015033_20170822_20170912_01_T1_B4.TIF"


#rasterio.windows.Window(col_off, row_off, width, height)
window = rasterio.windows.Window(2000,2000,800,600)

with rasterio.open(nirband) as src:
    subset = src.read(1, window=window)

fig, ax = plt.subplots(figsize=(12,6))
plt.imshow(subset)
plt.title(f'Band 5 Subset')





with rasterio.open(nirband) as src:
    nir = src.read(1, window=window)

with rasterio.open(redband) as src:
    red = src.read(1, window=window)

red = red.astype(float)
nir = nir.astype(float)
np.seterr(divide='ignore', invalid='ignore')

ndvi = np.empty(nir.shape, dtype=rasterio.float32)
check = np.logical_or ( red > 0, nir > 0 )
naip_ndvi = np.where ( check,  (1.0*(nir - red )) / (1.0*( nir + red )),-2 )


fig, ax = plt.subplots(figsize=(12,6))
ndvi = ax.imshow(naip_ndvi)
ax.set(title="NDVI")



with rasterio.open("LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF") as src:
    naip_data_ras = src.read()
    naip_meta = src.profile


with rasterio.open('MyExample.tif', 'w',**naip_meta) as dst:
    dst.write(naip_ndvi, window=window)


# =============================================================================
# with rasterio.open('example.tif') as dataset:
# 
#     # Read the dataset's valid data mask as a ndarray.
#     mask = dataset.dataset_mask()
# 
#     # Extract feature shapes and values from the array.
#     for geom, val in rasterio.features.shapes(
#             mask, transform=dataset.transform):
# 
#         # Transform shapes from the dataset's own coordinate
#         # reference system to CRS84 (EPSG:4326).
#         geom = rasterio.warp.transform_geom(
#             dataset.crs, 'EPSG:4326', geom, precision=6)
# 
#         # Print GeoJSON shapes to stdout.
#         print(geom)
# =============================================================================

Вот как выглядит NDVI, когда я использую matplotlib (я хочу сохранить это на своем рабочем столе в виде файла GeoTIFF):

NDVI

Спасибо за любую помощь!

Answer 1

Ваш src.profile (из «LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF»), вероятно, имеет тип dtype int32 или что-то подобное, в то время как ваш naip_ndvi содержит значения с плавающей запятой.

Вы должны установить dtype для вашего dst на float32 (или float64)

      with rasterio.open("LC08_L1TP_015033_20170822_20170912_01_T1_B5.TIF") as src:
    naip_data_ras = src.read()
    naip_meta = src.profile
    
    # set dtype to float
    naip_meta["dtype"] = rasterio.float32

Answer 2

Как вы смотрите на выходное изображение? В средстве просмотра изображений, в ГИС или программном обеспечении дистанционного зондирования, которые могут добавить контрастное растяжение к файлу? Значения NDVI варьируются от -1 до 1 - возможно, диапазон значений слишком мал для вашего программного обеспечения для автоматического отображения. Недавно у меня была похожая проблема с изменением изображений PlanetScope - они отображались, как и ожидалось, с помощью matplotlib, но tiff выглядел черным.

Вы можете попробовать масштабировать вывод, умножив значения ячеек на 100 - это может помочь с проблемами отображения. Вы также можете проверить значения выходного изображения с помощью программного обеспечения, которое может применить растягивание контраста к изображению (QGIS, продукты esri, ImageJ или программное обеспечение для обработки изображений).