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)require(RCurl)## Loading required package: RCurl
## Loading required package: bitopsrequire(raster)## Loading required package: raster
## Loading required package: sprequire(knitr)## Loading required package: knitr## 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")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
acc <- raster(gsub("basin",basin_name,"../acc_basin.asc"))## rgdal: version: 0.9-2, (SVN revision 526)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected) plot(acc)
# é precisso saber as coordenadas dos ponto Superior-Esquerdo (UL) e Inferior-Dereito (LR)
UL <- c(acc@extent@xmin,acc@extent@ymax);UL## [1] -2901.520 3677.145 LR <- c(acc@extent@xmax, acc@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.Nov.2004 <- list.files(pattern = "^M.*.2005-05-.*.Lai_1km.tif$")
lapply(Lai.Nov.2004, function(x) plot(raster(x), main=x))
setwd(script_dir) modis_process_path <- "../data_base/MODIS_process/"
dir.create(modis_process_path)## Warning in dir.create(modis_process_path): '../data_base/MODIS_process'
## already exists setwd(modis_process_path)
Lai.Nov.2004 <- list.files(pattern = "^M.*.2005-05-.*.Lai_1km.tif$")
lapply(Lai.Nov.2004, function(x) plot(mask(raster(x)*0.1,acc), main=x))
## [[1]]
## NULL
##
## [[2]]
## NULL
##
## [[3]]
## NULL
##
## [[4]]
## NULLsetwd(script_dir)