Я смог построить и отобразить свои растровые изображения с помощью 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):
Спасибо за любую помощь!