#realizat : Alexandru Dumitrescu #data : 2012-01-03 #descriere : prelucrarea produsului LST MODIS #email : alexandru.dumitrescu [at] gmail.com library(RCurl) library(raster) library(rgdal) gdal<-"/Library/Frameworks/GDAL.framework/Programs" #setwd("/Users/alexdum/Documents/Intro_Spatial") # 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) #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') } #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="")) #citire prin GDAL system(paste(gdal,'/gdalwarp HDF4_EOS:EOS_GRID:\"',files[1],'\":MODIS_Grid_Daily_1km_LST:LST_Day_1km fil1.tif', sep="")) system(paste(gdal,'/gdalwarp HDF4_EOS:EOS_GRID:\"',files[2],'\":MODIS_Grid_Daily_1km_LST:LST_Day_1km fil2.tif', sep="")) ls1<-readGDAL("fil1.tif") ls1$band1<-ls1$band1*0.02 ls2<-readGDAL("fil2.tif") ls2$band1<-ls2$band1*0.02 # 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", progress="text") #vizualizarea rezultatelor cu frontiere plot(srm,xlim=c(20,30),ylim=c(43.5,48.3),col=rev(heat.colors(16))) download.file('http://earth.unibuc.ro/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 GEOTIFF, 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)