• Чтение космических снимков в среде 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')
     
    

    Просмотрим информацию о снимке. Данный снимок представляет собой композит, есть информация о привязке, проекции, каналах.
    R_open_remotedata

    Теперь отобразим наш снимок.
     
     
    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ей строке с расположением каналов можно поиграться, тем самым посмотреть, как картинка будет меняться.

    R_open_remotedata2 Можно выбрать определенные группы каналов, используя функцию подмножества, узнать их названия.
     
     
    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) – отображение снимки изменилось под необходимый экстент.

    R_open_remotedata3

    Функциями 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()
     
    

    Опубликовано · Автор