Использование R для рассчета вегетационных индексов

[_taras_]

Доброго времени!
Намеди пришлось считать вегетационные снимки за несколько лет. В результате сделал некое подобие скрипта для автоматизации подсчетов, удаления облаков и обрезки по полигону, а также сохранению в GeoTIF.
Делался для Landsat 4-7, но трансформировать под другие спутники думаю не составит труда.
Буду благодарен любым конструктивным замечаниям и дополнениям.

[gamm]

единственное замечание общего характера. С данными типа raster R работает своеобразно, гоняя их на диск и обратно, что очень долго (это я проверял, даже когда-то покупал под временные файлы ssd, не сильно помогает). Если памяти хватает, то лучше сразу затащить данные в память, а потом с ними работать.
Типа такого:

Код: Выделить всё

B3 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B3.TIF"))- 0.2
B4 = 0.0000275*(raster ("LT05_L2SP_180028_19940304_20200913_02_T2_SR_B4.TIF"))- 0.2
NIR=B4[]
Red=B3[]
NDVI = (NIR − Red)/(NIR + Red)
r.NDVI=B1 # чтобы не мудрить с созданием
r.NDVI[]=NDVI
writeRaster(r.NDVI,out.fn)

Еще говорят, что brick в отличие от raster держит данные в памяти, но я не проверял.

[Константин Силкин]

Что сказать, Тарас… Раз уж поделились наработками, то молодец! А смысл? Как это может помочь тем, кто не знает что делать для вычисления NDVI? Ведь нельзя ваш скрипт быстро приспособить под свои нужды. Может только одна польза в том, что собрали все знакомые вам индексы в одном месте.
Добавьте тогда уж ещё водорослевый индекс NDAI.
Напишите лучше статью о том как вы писали этот скрипт. Там снабдите всё ссылками, разъяснениями, комментариями…

[_taras_]

В комплекте снимков Landsat-ов 2 уровня имеются снимки с классификацией облаков, теней и т.п. (*SR_CLOUD_QA.TIF). Расшифровку значений и дополнительную информацию можно посмотреть по адресу https://d9-wret.s3.us-west-2.amazonaws. ... ide-v4.pdf стр 13.
Это позволяет автоматизировать создание маски облаков

Код: Выделить всё

#######################ФОРМИРУЕМ РАСТР ОБЛАКОВ###################################
cloud_mask <- raster("LT05_L2SP_181028_19930220_20200914_02_T1_SR_CLOUD_QA.TIF")
cloud_mask[cloud_mask == 0] <- 1
cloud_mask[cloud_mask == 2] <- 0
cloud_mask[cloud_mask == 4] <- 0
cloud_mask[cloud_mask == 8] <- 0
cloud_mask[cloud_mask == 34] <- 0
cloud_mask[cloud_mask != 0] <- 1
#plot(cloud_mask)
writeRaster(cloud_mask, "/home/username/GRIDS/CLOUD_mm_yy.tiff")
CLOUD = raster ("CLOUD_mm_yy.tiff")

Из плюсов меньше тратится времени на обработку и единый подход к удалению облачного покрова.
Минус - плохо сделано распознавание туч и снега на зимних снимках и песчаных и ракушняковых пляжей и облаков на летних. Но эти проблемы в той или иной степени наблюдаются у всех методов выделения облаков