Załadujmy potrzebne pakiety.
library(WicPomiarSatelitarny)
Pobierzmy zdjęcia od 29 czerwca do dziś, ale tylko takie, dla których zachmurzenie nie przekraczało 60%.
(pozycja Puław to z grubsza 21.95 E, 51.5 N)
zdjPulawy = rbind( pobierzZdjecia(21.95, 51.5, "2016-06-29", "2016-06-29", 60, "dane/Pulawy", "wget"), pobierzZdjecia(21.95, 51.5, "2016-08-08", "2016-08-08", 60, "dane/Pulawy", "wget") )
# stwórz listę obiektów R na podstawie zdjęć zapisanych na dysku zdjecia = wczytajZdjecia("dane/Pulawy", "2016-08-08") # skomponuj obraz kolorowy z trzech wskazanych pasm eksportujKolor(zdjecia[["B04"]], zdjecia[["B03"]], zdjecia[["B02"]], 'dane/Pulawy/rgb.tif', 0, 5000) # pokoloruj klasyfikację scen w zadany sposób eksportujMono(zdjecia[["SCL"]], "dane/Pulawy/scl.tif", 0, 11, c('red', 'red', 'white', 'white', 'green', 'brown', 'blue', 'white', 'white', 'white', 'white', 'white')) # pokoloruj pasmo wegetacji na zielono eksportujMono(zdjecia[["B8A"]], "dane/Pulawy/B8A_kolor.tif", 0, 5000, c('black', 'green'))
Policzmy wskaźnik NDVI.
zdjecia[["NDVI"]] = (zdjecia[["B8A"]] - zdjecia[["B04"]]) / (zdjecia[["B8A"]] + zdjecia[["B04"]]) eksportujMono(zdjecia[["NDVI"]], "dane/Pulawy/NDVI_2016-08-08.tif", -1, 1, c('blue', 'black', 'green'))
Analogicznie policzmy NDVI dla 2016-06-29
zdjecia2 = wczytajZdjecia("dane/Pulawy", "2016-06-29") zdjecia2[["NDVI"]] = (zdjecia2[["B8A"]] - zdjecia2[["B04"]]) / (zdjecia2[["B8A"]] + zdjecia2[["B04"]]) eksportujMono(zdjecia2[["NDVI"]], "dane/Pulawy/NDVI_2016-06-29.tif", -1, 1, c('blue', 'black', 'green'))
Ustawmy na braki danych wszystkie piksele niesklasyfikowane jako "roślinność", "goła ziemia" ani "woda" (liczenie NDVI dla chmur, braków danych, itp. nie ma sensu).
zdjecia[["NDVI"]][! zdjecia[["SCL"]] == 4] = NA zdjecia2[["NDVI"]][! zdjecia2[["SCL"]] %in% c(4, 5, 6)] = NA eksportujMono(zdjecia[["NDVI"]], "dane/Pulawy/NDVI_bez_chmur_2016-08-08.tif", -1, 1, c('blue', 'black', 'green')) eksportujMono(zdjecia2[["NDVI"]], "dane/Pulawy/NDVI_bez_chmur_2016-06-29.tif", -1, 1, c('blue', 'black', 'green'))
Oznaczmy miejsca, wg różnicy wartości NDVI pomiędzy 29 czerwca i 8 sierpnia.
roznica = zdjecia[["NDVI"]] - zdjecia2[["NDVI"]] eksportujMono(roznica, "dane/Pulawy/roznica.tif", -1, 1, c('red', 'black', 'green'))
W maju 2016 r. w okolicach miasta McMurray w Kanadzie miał miejsce olbrzymi pożar lasu (trwał ponad miesiąc). Porównując NDVI sprzed i po pożarze spróbujemy odnaleźć pogorzeliska.
Współrzędne miasteczka to z grubsza 111.37 W, 56.7 N. Pobierzemy dane z 12 kwietnia i 21 czerwca.
zdjMcMurray = rbind( pobierzZdjecia(-111.37, 56.7, "2016-04-12", "2016-04-12", 60, "dane/Pulawy", "wget"), pobierzZdjecia(-111.37, 56.7, "2016-06-21", "2016-06-21", 60, "dane/Pulawy", "wget") )
Wczytajmy obrazy.
przed = wczytajZdjecia("dane/McMurray/", "2016-04-12") po = wczytajZdjecia("dane/McMurray/", "2016-06-21")
Policzmy NDVI i jego różnicę.
ndviPrzed = (przed[["B8A"]] - przed[["B04"]]) / (przed[["B8A"]] + przed[["B04"]]) ndviPo = (po[["B8A"]] - po[["B04"]]) / (po[["B8A"]] + po[["B04"]]) roznica = ndviPo - ndviPrzed
Za pomocą klasyfikacji scen zamieńmy na braki danych piksele przysłonięte chmurami na dowolnym z obrazów (różnica wartości NDVI dla tych pikseli nic nam nie powie).
roznica2[! przed[["SCL"]] %in% 4:6 | ! po[["SCL"]] %in% 4:6] = NA
Wyeksportujmy obrazki.
Obszary oznaczone na czerwono na obrazie "roznica2.tiff" w północnej części zdjęcia to prawdopodobnie pogorzeliska.
eksportujMono(roznica, 'dane/McMurray/roznica.tiff', -1, 1, c('red', 'black', 'green')) eksportujMono(roznica2, 'dane/McMurray/roznica2.tiff', -1, 1, c('red', 'black', 'green')) eksportujMono(ndviPrzed, 'dane/McMurray/ndvi_przed.tiff', -1, 1, c('blue', 'black', 'green')) eksportujMono(ndviPo, 'dane/McMurray/ndvi_po.tiff', -1, 1, c('blue', 'black', 'green')) eksportujMono(przed[["SCL"]], 'dane/McMurray/scl_przed.tiff', 0, 11, c('red', 'red', 'white', 'white', 'green', 'brown', 'blue', 'white', 'white', 'white', 'white', 'white')) eksportujMono(po[["SCL"]], 'dane/McMurray/scl_po.tiff', 0, 11, c('red', 'red', 'white', 'white', 'green', 'brown', 'blue', 'white', 'white', 'white', 'white', 'white')) eksportujKolor(przed[["B04"]], przed[["B03"]], przed[["B02"]], 'dane/McMurray/przed_rgb.tiff', 0, 5000) eksportujKolor(po[["B04"]], po[["B03"]], po[["B02"]], 'dane/McMurray/po_rgb.tiff', 0, 5000)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.