Problem z obliczaniem NDVI przy użyciu Rasterio Python

Question

Próbuję obliczyć NDVI za pomocą dwóch przyciętych obrazów rastrowych Landsat 7 (NIR i czerwone pasma przycięte za pomocą pliku maski) przy użyciu następującego kodu:

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

with rio.open(r'D:\clip_test_b3.tif') as src:
    red = src.read(1)   # (Rows, Columns) = (2731, 3660)
    
with rio.open(r'D:\clip_test_b4.tif') as src:
    nir = src.read(1)   # (Rows, Columns) = (2730, 3635)
    
np.seterr(divide = 'ignore', invalid = 'ignore')

ndvi = (nir.astype(float) - red.astype(float))/(nir + red)

plt.imshow(ndvi)

W powyższym kodzie oba paski (czerwone i NIR) mają różne kształty (różne rzędy i kolumny). Po uruchomieniu powyższego kodu otrzymuję komunikat „ValueError: operandy nie mogły być transmitowane razem z kształtami (2730,3635) (2731,3660)”.

Ale kiedy te same obliczenia NDVI, które próbuję wykonać w ArcMap (używając kalkulatora rastrowego), obliczane są NDVI.

Czy ktoś może mi pomóc w rozwiązaniu tego błędu.

Answer 1

Aby rozwiązać problem, musisz upewnić się, że kratki pokrywają ten sam obszar i mają takie same wymiary. Jednym ze sposobów osiągnięcia tego celu jest z reproject_match metodą in rioxarray(geoprzestrzennych przedłużenie xarray zasilany przez rasterio).

import rioxarray

red = rioxarray.open_rasterio("D:\clip_test_b3.tif")
nir_original = rioxarray.open_rasterio("D:\clip_test_b4.tif")

nir = nir_original.rio.reproject_match(red)
ndvi = (nir.astype(float) - red.astype(float))/(nir + red)

ndvi.plot()