Легко сегментируйте покров и почву с помощью NDVI и Rasterio

В этом посте мы попытаемся сегментировать растительный покров и почву на спутниковых снимках.

Мы будем заимствовать идеи из этой бумаги. Я также буду использовать идеи из моего предыдущего сообщения в блоге по этой теме.

В идеале мы хотим перейти от обычного спутникового снимка:

К этому:

Апельсин - это почва. Красный - это растительность.

NDVI

Как отмечается в документе, нам потребуется извлечь Нормализованный разностный индекс растительности. Это полезный индекс для растительности. Вот пример изображения NDVI:

И вот формула

Вот мой код для получения изображения NDVI в виде массива numpy.

def get_ndvi(red_file):

    nir_file = os.path.dirname(red_file) + '/nir.tif'
    band_red = rasterio.open(red_file)
    band_nir = rasterio.open(nir_sfile)
    red = band_red.read(1).astype('float64')
    nir = band_nir.read(1).astype('float64')

    # ndvi calculation, empty cells or nodata cells are reported as 0
    ndvi=np.where( (nir==0.) | (red ==0.), -255 , np.where((nir+red)==0., 0, (nir-red)/(nir+red)))

    return ndvi

NDVI рассчитывается для определенной области интересов. Это определяется в статье:

Для каждого полевого участка область интереса (ROI) была установлена ​​вручную путем выбора двух центральных строк, и было извлечено среднее значение индексов растительности, соответствующих каждому участку.

Это довольно просто реализовать на Python.

# get average ndvi of center two rows
  def get_region_of_interest(ndvi, multiplier = 1/2):

    # undo the background adjustment
    region = ndvi.copy()
    region = np.where(region==-255, 0, region)

    # mean of center rows
    center_row1 = np.mean(region[int((multiplier)*len(region))])
    center_row2 = np.mean(region[int((multiplier)*len(region))+1])

    # mean of both rows
    mean = (center_row1.copy()+center_row2.copy())/2

    return mean

Фракционный растительный покров

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

THRESHOLD = 0.3

def get_fc(ndvi):

    ndvi_copy = ndvi.copy()

    vegetation = np.where(ndvi_copy > THRESHOLD, 1, 0)
    vegetation_count = np.count_nonzero(vegetation)

    total = ndvi_copy.shape[0]*ndvi_copy.shape[1]
    fractional_cover = vegetation_count/total

    return fractional_cover

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

График доли растительного покрова в сравнении с NDVI

Теперь мы можем воссоздать сюжеты на седьмой странице бумаги. Мы построим график зависимости доли растительности от NDVI для каждого изображения.

Мы также хотим построить линию регрессии по методу наименьших квадратов.

Вот фрагмент кода, который позволяет нам это сделать.

def plot_fc_vs_ndvi(fc, ndvi):

    y = np.array(fc).reshape(1, -1)
    x = np.array(ndvi).reshape(1,-1)

    slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)

    x = np.linspace(min(ndvi),max(ndvi),100)
    f, ax = plt.subplots(1,1,figsize=(10,10))
    ax.plot(x, slope*x+intercept, '-r', label='fc='+str(round(slope, 2))+'*ndvi+'+str(round(intercept, 2)), color='black')
    ax.set_title('Fractional Cover vs NDVI at threshold of '+ str(THRESHOLD))

    scatter = ax.scatter(x=ndvi, y=fc, edgecolors='black')

    ax.set_xlabel('Normalized difference vegetation index (NDVI)')
    ax.set_ylabel('Fractional Cover (fc)')

    ax.text(min(ndvi)+0.8*(max(ndvi)-min(ndvi)), min(fc)+0.2*(max(fc)-min(fc)),s='R^2 = {}'.format(round((r_value**2), 4)), fontdict={'fontsize':14, 'fontweight':'bold'})
    f.savefig('fc_vs_ndvi_plots/fc_vs_ndvi_'+str(self.plot_counter)+'.jpg')
    f.show()

См. Полный код для более подробной информации.

Мы можем изменить установленное ранее значение threshold и посмотреть, как это повлияет на регрессию. В документе отмечается, что мы должны выбрать значение threshold, которое имеет лучшую модель регрессии.

Мы запустим этот код для всех изображений и построим график результатов.

Из этого мы видим, что самый высокий R² связан с порогом 0,45.

Мы можем применить этот порог ко всем изображениям. Любое значение NDVI, превышающее 0,45, относится к растительности. Все, что ниже 0,45, является почвой.

Создать двоичный массив

Используя порог 0,45, мы можем создать двоичный массив. В массиве 1 обозначает растительность. 0 обозначает почву.

Этот фрагмент кода делает именно это. Порог определяется функцией __init__. См. Подробности в полном коде.

def create_mask(self, red_file):

    nir_file = os.path.dirname(red_file) + '/nir.tif'
    band_red = rasterio.open(red_file)
    band_nir = rasterio.open(nir_file)
    red = band_red.read(1).astype('float64')
    nir = band_nir.read(1).astype('float64')

    # get raw ndvi and save as jpg
    self.raw_ndvi = np.where((nir+red)==0., np.nan, (nir-red)/(nir+red))
    
    # create canopy cover mask
    self.canopy_cover = np.where(np.isnan(self.raw_ndvi), np.nan, np.where(self.raw_ndvi<self.ndvi_threshold, 0, 1))
    self.canopy_cover = np.ma.masked_where(np.isnan(self.canopy_cover), self.canopy_cover)

    # show ndvi mask and save it as jpg
    print('canopy cover')
    print(np.unique(self.canopy_cover))

    return self.canopy_cover

Вот результат в более понятной форме. Первый ряд - это пороговое покрытие навеса. Вторая строка - это спутниковый снимок в формате RGB.

Мы видим, что он довольно неплохо отделяет растительный покров от почвы. Это намного лучше, чем моя предыдущая попытка кластеризации K-средних.

Вывод

В этом сообщении блога я описал способ сегментирования спутниковых изображений с помощью NDVI и растерио.

Я выполнял эту работу для небольшого стартапа в Сиднее. Без их помощи я бы не справился. Я многому у них научился.

Полный код находится на Github.

Я надеюсь, что это поможет людям. Если я допустил ошибку, напишите мне в твиттере. Спасибо!

Первоначально опубликовано на https://spiyer99.github.io 31 января 2021 г.