Prelucrarea datelor MODIS folosind limbajul R - produsele LST (temperatura suprafeței terestre active) și SC (extindere stratului de zăpadă)

Introducere

În acest articol sunt demonstrate capabilitățile limbajului R care poate rezolva, doar prin utilizarea funcțiilor proprii și a bibliotecilor suplimentare, etapele de descărcare, citire și prelucrare a datelor MODIS. Bibliotecile suplimentare care îndeplinesc cerințele menționate mai sus sunt RCurl, rgdal și raster.

Biblioteca RCurl este o implementare în R a aplicației libcurl renumită pentru facilitățile oferite în transferul fișierelor utilizând protocoalele HTTP, HTTPS, FTP, FTPS, etc. Această bibliotecă ne ajută să selectăm imaginile satelitare de interes găzduite pe serverele FTP puse la dispoziție de organizațiile care produc datele MODIS. Proprietățile celorlalte biblioteci (rgdal și raster) au fost descrise în primul articol dedicat limbajului R, cu ajutorul lor realizându-se prelucrarea efectivă a datelor satelitare.

Procedura de achiziție a produselor MODIS se poate realiza fie printr-o interfață web intuitivă (procedură descrisă deja printr-un articol existent pe geo-spatial.org), fie accesând direct datele disponibile pe serverele FTP. În acest tutorial se va descrie a doua procedură de achiziție a imaginilor satelitare, datele referitoare la temperatura suprafeței terestre active putând fi găsite la adresa ftp://e4ftl01.cr.usgs.gov/ iar informația privind extinderea stratului de zăpadă la adresa ftp://n4ftl01u.ecs.nasa.gov/. Prin utilitarul pus la dispoziție de UNH putem vizualiza tile-urile MODIS care acoperă suprafața României. Indecșii de interes sunt determinați de rândul 4 și coloanele 19 și 20 ale gridului sinusoidal MODIS (“h19v04”, “h20v04”), pentru acoperirea integrală a teritoriului României fiind necesară descărcarea a două imagini.

Etapa de selectare și descărcare a imaginilor de interes se poate realiza direct în R prin funcțiile RCurl. Funcțiile disponibile în biblioteca rgdal sunt utilizate în etapa de citire a datelor MODIS, din formatul nativ HDF-EOS, iar biblioteca raster ajută în procedura de mozaicare a celor două tile-uri care compun teritoriul României.

Din păcate varianta Windows a bibliotecii rgdal este compilată fără driver-ul HDF4, deci utilizatorii acestui sistem de operare nu pot folosi direct R-ul pentru citirea datelor MODIS. Utilizatorii MAC OS X trebuie sa instaleze variantă acestei biblioteci disponibilă aici. Este de asemenea necesar ca pentru sistemele cu LINUX și MAC OS X să fie instalată suita GDAL înainte de a se instala biblioteca rgdal pentru R.

Produsul MODIS LST (temperatura suprafeței terestre active)

Datele zilnice MODIS LST pot fi descărcate accesând adresa FTP mai sus menționată, căruia i se adaugă sufixul MOLT/MOD11A1.005 pentru date furnizate de senzorul Terra sau MOLA/MYD11A1.005 pentru senzorul Aqua. Cu ajutorul variabilelor an, luna și zi se va construi adresa FTP de unde vom extrage produsul MODIS pentru ziua și senzorul de interes.

library(RCurl)
library(raster)
library(rgdal)
# variabilele de selectie a produselor MODIS
url <- "ftp://e4ftl01.cr.usgs.gov/MOLT/MOD11A1.005/"
an<-"2007"
luna<-"07"
zi<-"27"
data<-paste(an,luna,zi,sep=".")
# creeaza adresa ftp unde se gasesc imaginile MODIS
ftp<-paste(url,data,"/",sep="")
fileNames <- getURL(ftp)

După ce adresa FTP a fost construită iar numele tuturor fișierelor care sunt găzduite extrase, se poate trece la descărcarea fișierelor care corespund teritoriului României.

#variabila cu indecșii imaginilor pentru Romania
p<-c("h19v04","h20v04")
# descarcarea imaginilor MODIS pentru Romania
for (i in p)
{print(i)
files = strsplit(fileNames,"\r*\n")[[1]]
# selectie imagini pentru Romania
ff<-files[grep(files,pattern=paste(i,"*.hdf",sep="."))]
ff<-ff[nchar(ff)==83]
ff<-unlist(strsplit(ff," "))
ff<-ff[length(ff)]
download.file(paste(url,data, "/", ff,sep=""),destfile=paste(getwd(), ff, sep="/"), mode='wb')
}

Cu ajutorul funcției GDALinfo() putem vizualiza proprietățile fișierelor *.hdf descărcate (numărul de benzi, numele fiecărei benzi, data de producere a imaginilor, etc). În cele ce urmează vom procesa informația LST extrasă ziua (MODIS_Grid_Daily_1km_LST:LST_Day_1km).

#lista cu fisierele *.hdf descarcate
files<-list.files(pattern=".hdf$",full.names=T)
#GDALinfo(files[1])
#citirea datelor MODiS in R
ls1<-readGDAL(paste('HDF4_EOS:EOS_GRID:\"',files[1],'\":MODIS_Grid_Daily_1km_LST:LST_Day_1km',sep=""))
ls2<-readGDAL(paste('HDF4_EOS:EOS_GRID:\"',files[2],'\":MODIS_Grid_Daily_1km_LST:LST_Day_1km',sep=""))
# mozaicarea datelor MODIS
rm<-merge(raster(ls1),raster(ls2))
rm<-crop(rm,extent(1400000,2500000,4600000,5500000))
#transforma din grade Kelvin in grade Celsius
rm<-rm-273.15
#proiecteaza din proiectie Sinusoidala in WGS84
srm<-projectRaster(rm, crs="+init=epsg:4326", method="ngb")

Pentru vizualizarea imaginilor mozaicate și reproiectate în proiecție Geografică WGS84 se folosește funcția plot(), iar prin writeGDAL() se poate exporta informației LST în format raster GeoTIFF. Dacă nu se dorește arhivarea imaginilor originale *.hdf se poate utiliza funcția unlink().

#vizualizarea rezultatelor cu frontiere
plot(srm,xlim=c(20,30),ylim=c(43.5,48.3),col=rev(heat.colors(16)))
download.file('http://www.geo-spatial.org/vechi/file_download/1251','fro.zip',mode = "wb")
unzip('fro.zip')
ro<-readOGR(dsn=getwd(), layer="frontiera_ro")
plot(ro,add=T)
# scrie datele MODIS in format GEOTIF, cu numele corespunzator zilei
srm<-crop(srm,extent(20,30,43.5,48.3))
writeRaster(srm, paste("LST_",gsub("\\.","_",data),".tif",sep=""),overwrite=T)
##stergerea fisierelor HDF4 descarcate
unlink(files)

Dacă se dorește aplicarea liniilor de cod R prezentate mai sus pentru procesarea altui produs, trebuie schimbat doar sufixul din variabila ftp. De asemenea trebuie aflat numele benzii care se vrea a fi prelucrată (lucru posibil prin funcția GDALinfo()), această informație fiind necesară în procedura de citire a fișierelor *.hdf . O descriere detaliată a produselor MODIS disponibile pe serverul ftp://e4ftl01.cr.usgs.gov/ este oferită de https://lpdaac.usgs.gov/lpdaac/products/modis_products_table.

Figura 1. MODIS LST – 27 iulie 2007

Produsul MODIS SC (extinderea stratului de zăpadă)

Prelucrarea produsului SC se realizează după aceeași metodologie aplicată cazului LST, diferențele constând doar în setarea variabilei url și a numelui benzii care este citită în R.

Produsele referitoare la stratul de zăpadă pot fi găsite la adresa ftp://n4ftl01u.ecs.nasa.gov/, o descriere a datelor MODIS găzduite pe acest server fiind disponibilă pe portalul http://nsidc.org/data/modis/data_summaries/index.html. Numele benzii (MOD_Grid_Snow_500m:Snow_Cover_Daily_Tile) care va fi prelucrată este identificat utilizând funcția GDALinfo(), pentru acest material fiind ales produsul MYD10A1 din 26 decembrie 2011. Timpul de procesare în cazul datelor MODIS SC este mai mare fapt explicabil prin rezoluția spațială mai bună a acestui produs (500×500 m).

Figura 2. MODIS SC – 26 decembrie 2011

Concluzii

Cele prezentate mai sus demonstrează eficiența limbajului R care poate rezolva doar prin câteva linii de script achiziția și prelucrarea produselor MODIS. Doar prin modificarea a trei variabile (data, adresă FTP și numele benzii) scriptul prezentat în acest tutorial poate fi introdus în flux operativ și aplicat oricărui produs level 3 MODIS.

Chiar dacă în Windows momentan acest script nu se poate folosi integral, se poate rula doar etapa de achiziție iar pentru partea de procesare pot fi utilizate funcțiile din GDAL sau MRT (MODIS Repreojection Tool). Un articol foarte bun în care sunt prezentate câteva posibilități de utilizare a limbajului R, sub Windows, în tandem cu cele două suite (GDAL și MRT) poate fi consultat aici.

Comenzile prezentate mai sus pot fi rulate linie cu linie prin încărcare scripturilor în mediul R după cum urmează:

• pentru produsul LST – MODIS_temperatura.R
• pentru produsul SC – MODIS_zapada.R