Чтение космических снимков в среде R
В данной заметке рассмотрены возможности среды R для работы с космическими снимками. В качестве примера будем использовать снимки, полученные аппаратом Landsat 8. К задачам будут относится такие операции, как синтез каналов снимка, визуализация композитов. В качестве тренировки используются снимки, покрывающие территорию дельты Волги и зону западных подстепных ильменей.
В качестве теоретической базы для данной работы использовались такие зарубежные материалы, как:
- Analyzing Remote Sensing Data in R: The landsatPackage
- Introduction to remote sensing data analysis using R
Для начала нам необходимо убедиться, что у нас в среде RStudio установлены библиотеки, такие как «raster» и «cluster».
Итак, подключив необходимые пакеты, подключаемся к базовой папке с файлами. Например, у меня было так:
library(raster)
r <- brick('D:/delta/LC81690282016223LGN00_B432.tif')
Просмотрим информацию о снимке. Данный снимок представляет собой композит, есть информация о привязке, проекции, каналах.
nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE) plotRGB(r, r = 3, g = 2, b = 1, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")
plotRGB(r, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Landsat Customize Color Composite")
Во второй строке видим как отображены каналы. Данную расстановку менять не будет, так как при визуализации это будет градация «естественные цвета». А вот в 3ей строке с расположением каналов можно поиграться, тем самым посмотреть, как картинка будет меняться.
Можно выбрать определенные группы каналов, используя функцию подмножества, узнать их названия.
rsub <- subset(r, 2:3)
nlayers(r)
nlayers(rsub)
names(r)
Мы можем задать необходимый экстент снимка, изменить необходимую территорию для визуализации.
extent(r) e <- extent(702750, 800815, 5100150, 5158583)
rr <- crop(r, e)
nf <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(rr, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Landsat Subset")
Обратите внимание, что после того как мы узнали настоящий экстент, изменили его – e <- extent(702750, 800815, 5100150, 5158583)
– отображение снимки изменилось под необходимый экстент.
Функциями aggregate и disaggreagte можно изменять пространственное разрешение или размер пикселя.
rra <- aggregate(rr, fact=3)
nfd <- layout(matrix(c(1,0,2), 1, 3, byrow = TRUE), width = c(1,0.2,1), respect = TRUE)
plotRGB(rr, r = 4, g = 3, b = 2, axes = TRUE, stretch = "lin", main = "Landsat 30m")
dev.off()
Опубликовано · Автор Victor