DBHM-GIS: PARTE 8
05-11-2014
Universidade Federal de Santa Maria - UFSM
Departamento de Física
Prof. Jônatan Tatsch
Aluno Roilan Hernandez Valdes
Download dos dados de LAI e FPAR.
Os dados serão obtidos do site do MODIS, usando as funções desenvolvidas para esse fim.
#instalando pacotes necessários
install.packages("RCurl", dependencies = TRUE)
install.packages("rgdal", dependencies = TRUE)
install.packages("climstats", repos="http://R-Forge.R-project.org",dependencies = TRUE)
install.packages("snow",dependencies = T)
install.packages("lmom", dependencies = T)require(RCurl)
require(raster)
require(knitr)
require(rgdal)
#require(zoo)
# require(R.utils)
require(climstats)## configurando opção para mostrar barra de progresso durante operações com raster
rasterOptions(progress = "text")
## definindo-se o seu nome de usuario que establecer o caminho até o diretório scripts
user_name <- system("echo $USER", intern = TRUE )
## caminho para o diretório do DBHM
## substitui a string user pela valor da variável user_name
dbhm_path <- gsub(pattern = "user",
replacement = user_name,
x = "/home/user/DBHM/")
dbhm_path## [1] "/home/hidro2roilan/DBHM/"## definindo diretório de trabalho (combina my_path e caminho até subdir scripts)
script_dir <- paste0(dbhm_path, "data_process/geo_data/scripts")
script_dir## [1] "/home/hidro2roilan/DBHM/data_process/geo_data/scripts"## define diretório de trabalho
setwd(script_dir)
## verifica diretório de trabalho configurado
getwd()## [1] "/home/hidro2roilan/DBHM/data_process/geo_data/scripts"## script para escrita no formato ascii DBHM
source("raster2ascii_dbhm.R")
source("brickstack_to_saved_rasters2.R")Instalação de sofware do MRTools
Existe o software MRTools disponível multiplataforma, aquí só se explica o procedimento para o SO Linux 64-bits. Precisa-se acessar ao site da NASA, fazer o cadastro e escolher o .zip em dependência de seu sistema. Os comandos de ajuda podem ser:
$ unzip MRT_download_Linux64.zip
$ cd MRT_download_Linux64/
$ ./mrt_installA continuação começa a instalação interativa do programa que perguntará por várias informações:
Do you wish to proceed with the MRT installation? [y/n]
y (enter)
Enter the MRT directory path:
/home/user/.MRT (enter)
Directory does not exist. Create /home/user/.MRT? [y/n]
y (enter)Ao final da instalação na terminal sai a declaração das variaveis de entorno, essas quatro linhas devem ser copiadas no arquivo .bashrc no $HOME do usuário
MRT_HOME=/home/user/.MRT
PATH=$PATH:/home/user/.MRT/bin
MRT_DATA_DIR=/home/user/.MRT/data
export MRT_HOME PATH MRT_DATA_DIRApós copiar digitar source ~/.bashrc na terminal.
MODIS
A configuração dos tiles do MODIS é mostrada na siguiente figura
Uma documentação mais detalhada para o download específico do produtos do MODIS pode-se encontrar no site de r-gis essa operação é seguida aquí. Accessar ao site ou baixar os siguientes arquivos: ModisDownload.R e ModisLP.RData e guardar na pasta do directório de trabalho.
source("ModisDownload.R")
load("ModisLP.RData")
# Para mosrar a tabela com os produtos disponíveis
tablemodis <- modisProducts()
tablemodis[,c(1,3,5,6)]## product name resolution temporal
## 1 MOD13Q1 Vegetation Indices 250m 16 day
## 2 MYD14A1 Thermal Anomalies & Fire 1000m Daily
## 3 MOD14A1 Thermal Anomalies & Fire 1000m Daily
## 4 MYD14A2 Thermal Anomalies & Fire 1000m 8 day
## 5 MOD14A2 Thermal Anomalies & Fire 1000m 8 day
## 6 MYD14 Thermal Anomalies & Fire 1000m 5 min
## 7 MYD09CMG Surface Reflectance Bands 1\x967 5600m Daily
## 8 MOD09CMG Surface Reflectance Bands 1\x967 5600m Daily
## 9 MYD09GQ Surface Reflectance Bands 1\x962 250m Daily
## 10 MOD09GQ Surface Reflectance Bands 1\x962 250m Daily
## 11 MYD09GA Surface Reflectance Bands 1\x967 500/1000m Daily
## 12 MYD09A1 Surface Reflectance Bands 1\x967 500m 8 day
## 13 MYD17A2 Gross Primary Productivity 1000m 8 day
## 14 MOD09A1 Surface Reflectance Bands 1\x962 500m 8 day
## 15 MYD09Q1 Surface Reflectance Bands 1\x962 250m 8 day
## 16 MOD09Q1 Surface Reflectance Bands 1\x962 250m 8 day
## 17 MCD43B4 Nadir BRDF-Adjusted Reflectance 1000m 16 day
## 18 MCD43A4 Nadir BRDF-Adjusted Reflectance 500m 16 day
## 19 MCD43C4 Nadir BRDF-Adjusted Reflectance 5600m 16 day
## 20 MYD15A2 Leaf Area Index - FPAR 1000m 8 day
## 21 MOD15A2 Leaf Area Index - FPAR 1000m 8 day
## 22 MCD15A3 Leaf Area Index - FPAR 1000m 4 day
## 23 MOD44W Land Water Mask Derived 250m None
## 24 MYD11C3 Land Surface Temperature & Emissivity 5600m Monthly
## 25 MOD11C3 Land Surface Temperature & Emissivity 5600m Monthly
## 26 MYD11B1 Land Surface Temperature & Emissivity 5600m Daily
## 27 MOD11B1 Land Surface Temperature & Emissivity 5600m Daily
## 28 MYD11A1 Land Surface Temperature & Emissivity 1000m Daily
## 29 MYD11C1 Land Surface Temperature & Emissivity 5600m Daily
## 30 MOD11C1 Land Surface Temperature & Emissivity 5600m Daily
## 31 MOD11A2 Land Surface Temperature & Emissivity 1000m 8 day
## 32 MYD11C2 Land Surface Temperature & Emissivity 5600m 8 day
## 33 MOD11C2 Land Surface Temperature & Emissivity 5600m 8 day
## 34 MYD11_L2 Land Surface Temperature & Emissivity 1000m 5 min
## 35 MCD12Q1 Land Cover Type 500m Yearly
## 36 MCD12C1 Land Cover Type 5600m Yearly
## 37 MCD12Q2 Land Cover Dynamics 500m Yearly
## 38 MOD17A2 Gross Primary Productivity 1000m 8 day
## 39 MCD43C2 BRDF-Albedo Snow-free Quality 5600m 16 day
## 40 MCD43A2 BRDF-Albedo Quality 500m 16 day
## 41 MCD43A1 BRDF-Albedo Model Parameters 500m 16 day
## 42 MCD43B1 BRDF-Albedo Model Parameters 1000m 16 day
## 43 MCD43C1 BRDF-Albedo Model Parameters 5600m 16 day
## 44 MCD43B3 Albedo 1000m 16 day
## 45 MCD43C3 Albedo 5600m 16 day
## 46 MYD13C2 Vegetation Indices 5600m Monthly
## 47 MOD13C2 Vegetation Indices 5600m Monthly
## 48 MYD13C1 Vegetation Indices 5600m 16 day
## 49 MOD13C1 Vegetation Indices 5600m 16 day
## 50 MYD13A3 Vegetation Indices 1000m Monthly
## 51 MOD13A3 Vegetation Indices 1000m Monthly
## 52 MYD13Q1 Vegetation Indices 250m 16 day
## 53 MYD13A2 Vegetation Indices 1000m 16 day
## 54 MOD13A2 Vegetation Indices 1000m 16 day
## 55 MYD13A1 Vegetation Indices 500m 16 day
## 56 MCD43B2 BRDF-Albedo Quality 1000m 16 day
## 57 MYD11A2 Land Surface Temperature & Emissivity 1000m 8 day
## 58 MCD45A1 Burned Area 500m Monthly
## 59 MOD11A1 Land Surface Temperature & Emissivity 1000m Daily
## 60 MOD44B Vegetation Continuous Fields 250m Yearly
## 61 MOD14 Thermal Anomalies & Fire 1000m 5 min
## 62 MOD17A3 Gross Primary Productivity 1000m Yearly
## 63 MCD15A2 Leaf Area Index - FPAR 1000m 8 day
## 64 MOD13A1 Vegetation Indices 500m 16 day
## 65 MCD43A3 Albedo 500m 16 day
## 66 MOD11_L2 Land Surface Temperature & Emissivity 1000m 5 min
## 67 MOD09GA Surface Reflectance Bands 1\x967 500/1000m DailyPara os interesses do DBHM são necessários os produtos LAI (MCD15A2 subset 2 (/6)), FPAR (MCD15A2 subset 1 (/6)), NDVI (MYD13A2 MOD13A2 subset 1 (/11))
Procurando a informação dos arquivos da bacia.
#Definindo a bacia de trabalho:
basin_name <- "pauliceia"
#O arquivo de Area Acumulada pode ser usado para obter a informação necessária
wsdem <- raster(gsub("basin",basin_name,"../wsdem_basin.asc"))
plot(wsdem)
# é precisso saber as coordenadas dos ponto Superior-Esquerdo (UL) e Inferior-Dereito (LR)
UL <- c(wsdem@extent@xmin,wsdem@extent@ymax);UL## [1] -2901.520 3677.145 LR <- c(wsdem@extent@xmax, wsdem@extent@ymin);LR## [1] 3098.480 -3322.855# Procurando a informação da projeção.
laea_file <- gsub("basin",basin_name,"../data_base/basin/mdet_basin_laea_int.tif")
laea <- projection(raster(laea_file))
laea## [1] "+proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs" split_laea <- strsplit(x = laea, split = c(" "))[[1]]
split_laea## [1] "+proj=laea" "+lat_0=-21.64" "+lon_0=-47.63" "+x_0=0"
## [5] "+y_0=0" "+a=6370997" "+b=6370997" "+units=m"
## [9] "+no_defs" lat_laea <-as.numeric(strsplit(x = split_laea[2], split= "=")[[1]][2])
lat_laea## [1] -21.64 lon_laea <-as.numeric(strsplit(x = split_laea[3], split= "=")[[1]][2])
lon_laea ## [1] -47.63 raio_laea <-as.numeric(strsplit(x = split_laea[6], split= "=")[[1]][2])
raio_laea## [1] 6370997 proj_params <- paste(as.character(raio_laea),"0.0 0.0 0.0",as.character(lon_laea), as.character(lat_laea),
"0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0",
sep = " ")
proj_params## [1] "6370997 0.0 0.0 0.0 -47.63 -21.64 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0"Chunk para fazer o download, reprojeção e o recorte da área.
modis_down_path <- "../data_base/MODIS/"
dir.create(modis_down_path )
setwd(modis_down_path )
ModisDownload(x = 'MCD15A2', #Código do produto
h = 13, #Coordenada horizontal do tile
v = 11, #Coordenada vertical do tile
dates =c('2003.01.01','2015.04.30'), #Data de inicial e final
mosaic=FALSE, #Flag para fazer mosaico se tem mais de um tile
MRTpath = '/home/hidro2roilan/.MRT/bin', #Diretorio onde o MRTools foi instalado
bands_subset = "1 1 1 1 1 1", #Bandas a ser processadas no reprojeção
del= TRUE, #Flag paa deletar os tiles descarregados
proj= TRUE, #Flag para fazer a reprojeção
resample_type="NEAREST_NEIGHBOR", #Tipo de reamostragem
proj_type="LA", #Tipo de projeção
proj_params=proj_params, #parámetros da projeção
UL = UL, #Coordenada Upper-Left do dominio
LR = LR, #Coordenada Low-Right do dominio
datum='NODATUM', #Datum ??
pixel_size= 500 #Resolução de saída
)
ModisDownload(x = 'MOD13A2',
h = 13,
v = 11,
dates =c('2003.01.01','2015.04.30'),
mosaic=FALSE,
MRTpath = '/home/hidro2roilan/.MRT/bin',
bands_subset = "1 1 1 1 1 1 1 1 1 1 1 1",
del= TRUE,
proj= TRUE,
resample_type="NEAREST_NEIGHBOR",
proj_type="LA",
proj_params=proj_params,
UL = UL,
LR = LR,
datum='NODATUM',
pixel_size= 500
)
ModisDownload(x = 'MYD13A2',
h = 13,
v = 11,
dates =c('2003.01.01','2015.04.30'),
mosaic=FALSE,
MRTpath = '/home/hidro2roilan/.MRT/bin',
bands_subset = "1 1 1 1 1 1 1 1 1 1 1 1",
del= TRUE,
proj= TRUE,
resample_type="NEAREST_NEIGHBOR",
proj_type="LA",
proj_params=proj_params,
UL = UL,
LR = LR,
datum='NODATUM',
pixel_size= 500
)
setwd(script_dir)Chunk para fazer a reprojeção y o recorte da área de interesse.
Este passo considera que os arquivos do MODIS necessários já foram obtidos, e estão no formato HDF.
## Usando as funções reprojeção do MRTools:
## listando todos os arquivos dos produtos
modis_products_path <- paste(script_dir, "../data_base/MODIS_download", sep="/")
MCD15A2_list <- list.files(path = modis_products_path, pattern = "^MCD15A2.*hdf$")
MOD13A2_list <- list.files(path = modis_products_path, pattern = "^MOD13A2.*hdf$")
MYD13A2_list <- list.files(path = modis_products_path, pattern = "^MYD13A2.*hdf$")
all_MODIS <- c(MCD15A2_list,MYD13A2_list,MOD13A2_list)
#Criando e deslocando à pasta de processamento dos arquivos para a bacia
modis_process_path <- "../data_base/MODIS_process/"
dir.create(modis_process_path)
setwd(modis_process_path)
#Fazendo enlaces simbólicos dos arquivos originais no diretório de trabalho
lapply(all_MODIS, function(x) system(paste("ln -s ",modis_products_path,"/",x," ",getwd() ,sep="")) )
#Função para o nome dos arquivos de saída (seja igual a saídas da função ModisDownload() )
filename <- function(x) {
pref <- strsplit(x,'/')[[1]]
pref <- pref[length(pref)]
pref <- strsplit(pref,'\\.')[[1]]
tero <- pref[2]
date <- as.POSIXct(paste(substr(tero,2,5),"-01-01",sep=""))
day <- as.numeric(substr(tero,6,8))
date_name <- as.Date(day-1, origin=as.Date(date))
pref <- paste(pref[1],"_",pref[3],sep="")
filename=paste(pref,'_',date_name,'.tif',sep='')
}
#Fazendo a reprojeção, interpolação e extração dos subprodutos
for( i in 1:length(all_MODIS)){
reprojectHDF(all_MODIS[i],
filename = filename(all_MODIS[i]),
MRTpath = '/home/hidro2roilan/.MRT/bin',
resample_type="NEAREST_NEIGHBOR",
proj_type="LA",
proj_params=proj_params,
UL = UL,
LR = LR,
datum='NODATUM',
pixel_size= 500)
}
#Deletando os enlaces simbólicos dos arquivos HDF
lapply(all_MODIS, function(x) system(paste("rm ",x) ) )
Lai.mm.yyyy <- list.files(pattern = "^M.*.2005-05-.*.Lai_1km.tif$")
lapply(Lai.mm.yyyy, function(x) plot(raster(x), main=x))
setwd(script_dir)Exemplo dos rasters de saída combinado com a área acumulada como máscara.
setwd(script_dir)
modis_process_path <- "../data_base/MODIS_process/"
ifelse(dir.exists(modis_process_path),getwd(),dir.create(modis_process_path) )## [1] "/home/hidro2roilan/DBHM/data_process/geo_data/scripts" setwd(modis_process_path)
plot(stack(list.files(pattern = "^M.*.2005-05-.*.Lai_1km.tif$")))
#Listando os arquivos LAI
lai_lst_files <- list.files(pattern = ".Lai_1km.tif$")
head(lai_lst_files)## [1] "MCD15A2_h13v11_2003-01-01.Lai_1km.tif"
## [2] "MCD15A2_h13v11_2003-01-09.Lai_1km.tif"
## [3] "MCD15A2_h13v11_2003-01-17.Lai_1km.tif"
## [4] "MCD15A2_h13v11_2003-01-25.Lai_1km.tif"
## [5] "MCD15A2_h13v11_2003-02-02.Lai_1km.tif"
## [6] "MCD15A2_h13v11_2003-02-10.Lai_1km.tif" #Determinando as datas dos arquivos
dates <- as.Date(substr(lai_lst_files,16,25))
#Aplicando máscara e juntando os arquivos
raster_object_from_files <- mask(stack(lai_lst_files[order(dates)]),wsdem)
# Eliminando dados espúrios e aplicando o gain do produto
b <- calc(raster_object_from_files, function(x){ x[x > 101] <- NA; return(x)})
# NAvalue(b) <- -9999
gain(b) <- 0.1
# NAvalue(b) <- -9999
names(b) <- gsub("-","", dates)
plot(b,1)
# raster_list_from_b <- brickstack_to_saved_rasters2(b, output_basename = tempdir)
#
# raster_object_from_b <- stack(raster_list_from_b)
#
raster_object_from_b <- setZ(b,dates,"Date/time")
inds <- as.factor(format(getZ(raster_object_from_b), "%Y-%m-15"))
result <- stackApply(x = raster_object_from_b,
indices = inds,
fun = max, na.rm = TRUE)##
|
| | 0%
|
|====================== | 33%
|
|=========================================== | 67%
|
|=================================================================| 100%
## v_result <- getValues(result)
v_result[is.infinite(v_result)] <- NA
result <- setValues(result, v_result)
avg_result <- cellStats(result, stat = mean)
result_mean <- stackApply(x = raster_object_from_b,
indices = inds,
fun = mean, na.rm = TRUE)##
|
| | 0%
|
|====================== | 33%
|
|=========================================== | 67%
|
|=================================================================| 100%
## v_result_mean <- getValues(result_mean)
v_result_mean[is.infinite(v_result_mean)] <- NA
result_mean <- setValues(result_mean, v_result_mean)
avg_result_mean <- cellStats(result_mean, stat = mean)
dts <- as.Date(as.character(inds))
plot(unique(dts),avg_result, type = "l")
# lines(unique(dts),avg_result_mean, type = "l", col = "blue")
names(result) <- gsub("-","",substr(as.character(unique(dts)),3,7))
# ## lista dos raster de maximos mensais
# ## acrescentar prefixo "lai", "fpar" e "max" para distincao
raster_list_from_maxlai <- brickstack_to_saved_rasters2(result,
output_basename = "LAI/lai_",
use_layernames=F)
for( i in 1:nlayers(result)){
YYMM <- format.Date(unique(dts)[i],format = "%y%m")
raster_lai_file <-list.files("LAI", pattern = gsub("YYMM",YYMM,".XYYMM.grd$"), full.names = T)
output_path <- paste(script_dir,"/../LAI_",basin_name, sep = "")
ifelse(dir.exists(output_path), getwd(), dir.create(output_path))
bin_lai_file <- paste(output_path,gsub("YYMM",YYMM,"/Lamb_YYMM_LAI.bin"), sep= "")
writeBin(getValues(raster(raster_lai_file)),bin_lai_file, size = 4)
} setwd(script_dir)
modis_process_path <- "../data_base/MODIS_process/"
ifelse(dir.exists(modis_process_path),getwd(),dir.create(modis_process_path) )## [1] "/home/hidro2roilan/DBHM/data_process/geo_data/scripts" setwd(modis_process_path)
plot(stack(list.files(pattern = "^M.*.2005-05-.*.Fpar_1km.tif$")))
#Listando os arquivos LAI
fpar_lst_files <- list.files(pattern = ".Fpar_1km.tif$")
head(fpar_lst_files)## [1] "MCD15A2_h13v11_2003-01-01.Fpar_1km.tif"
## [2] "MCD15A2_h13v11_2003-01-09.Fpar_1km.tif"
## [3] "MCD15A2_h13v11_2003-01-17.Fpar_1km.tif"
## [4] "MCD15A2_h13v11_2003-01-25.Fpar_1km.tif"
## [5] "MCD15A2_h13v11_2003-02-02.Fpar_1km.tif"
## [6] "MCD15A2_h13v11_2003-02-10.Fpar_1km.tif" dates <- as.Date(substr(fpar_lst_files,16,25))
#Aplicando máscara e juntando os arquivos
raster_object_from_files <- mask(stack(fpar_lst_files[order(dates)]),wsdem)
# Eliminando dados espúrios e aplicando o gain do produto
c <- calc(raster_object_from_files, function(x){ x[x > 101] <- NA; return(x)})
gain(c) <- 0.01
names(c) <- gsub("-","", dates)
plot(c,1)
raster_object_from_c <- setZ(c,dates,"Date/time")
inds <- as.factor(format(getZ(raster_object_from_c), "%Y-%m-15"))
result <- stackApply(x = raster_object_from_c,
indices = inds,
fun = max, na.rm = TRUE)##
|
| | 0%
|
|====================== | 33%
|
|=========================================== | 67%
|
|=================================================================| 100%
## v_result <- getValues(result)
v_result[is.infinite(v_result)] <- NA
result <- setValues(result, v_result)
avg_result <- cellStats(result, stat = mean)
result_mean <- stackApply(x = raster_object_from_c,
indices = inds,
fun = mean, na.rm = TRUE)##
|
| | 0%
|
|====================== | 33%
|
|=========================================== | 67%
|
|=================================================================| 100%
## v_result_mean <- getValues(result_mean)
v_result_mean[is.infinite(v_result_mean)] <- NA
result_mean <- setValues(result_mean, v_result_mean)
avg_result_mean <- cellStats(result_mean, stat = mean)
#
# dts <- as.Date(as.character(inds))
plot(unique(dts),avg_result, type = "l")
# lines(unique(dts),avg_result_mean, type = "l", col = "blue")
names(result) <- gsub("-","",substr(as.character(unique(dts)),3,7))
# ## lista dos raster de maximos mensais
# ## acrescentar prefixo "lai", "fpar" e "max" para distincao
ifelse(dir.exists("FPAR"),getwd(),dir.create("FPAR") )## [1] "/home/hidro2roilan/DBHM/data_process/geo_data/data_base/MODIS_process" raster_list_from_maxlai <- brickstack_to_saved_rasters2(result,
output_basename = "FPAR/fpar_",
use_layernames=F)
for( i in 1:nlayers(result)){
YYMM <- format.Date(unique(dts)[i],format = "%y%m")
raster_fpar_file <-list.files("FPAR", pattern = gsub("YYMM",YYMM,".XYYMM.grd$"), full.names = T)
output_path <- paste(script_dir,"/../FPAR_",basin_name, sep = "")
ifelse(dir.exists(output_path), getwd(), dir.create(output_path))
bin_fpar_file <- paste(output_path,gsub("YYMM",YYMM,"/Lamb_YYMM_FPAR.bin"), sep= "")
writeBin(getValues(raster(raster_fpar_file)),bin_fpar_file, size = 4)
}