DBHM-GIS: PARTE 3
16-11-2014
Universidade Federal de Santa Maria - UFSM
Departamento de Física
Prof. Jônatan Tatsch
Aluno Roilan Hernandez Valdes
1. Derivação das direções de fluxo de baixa resolução
1.1. Convenção: rasters de alta e baixa resolução
O modelo DBHM, assim como outros modelos hidrológicos distribuídos, representa a bacia por células quadradas, usualmente da ordem de cerca de 1 km de dimensão na projeção cartográfica LAEA (ou ~0,01° no sistema de projeções geográficas).
Nesse texto os elementos de maior resolução (ou seja, menores dimensões) estão sendo referenciados como pixels, enquanto que as células do modelo hidrológico compõem a grade de baixa resolução. A Figura 4 ilustra os conceitos de pixels e células.
Uma das exigências quanto à definição das resoluções alta e baixa é que as dimensões dos pixels e as dimensões das células do modelo sejam múltiplas entre si, de tal forma que em cada célula esteja contido o mesmo número inteiro de pixels. Por exemplo, usualmente se adota a dimensão ~1 km para a célula do DBHM, e a resolução de 100 m para os rasters de alta resolução. Nesse exemplo, em cada célula existirão 100 pixels internamente a uma célula (upscaling de 10 vezes). Para a mesma dimensão de célula do exemplo, outras resoluções altas possíveis seriam: 500 m (5 vezes), 300 m (3 vezes), 200 m (2 vezes), etc. Para bacias como a Mogu e Vacacai a alta resolução é 100 m e a baixa resolução será de 1000 m. Portanto se um raster de alta resolução conter 282 colunas, em baixa resolução ele terá 28,2 colunas? Foi para evitar esse inconveniente que foi realizado o ajuste do domínio espacial do MDET na seção 4 da Parte 1.
1.2. Descrição
A determinação das direções de fluxo e das áreas de drenagem acumuladas de alta resolução tem como objetivo principal servir como informação para a derivação das direções de fluxo de baixa resolução, isto é, a direção de fluxo de cada célula do modelo hidrológico. Esse processo é conhecido como upscaling de direções de fluxo e tem apresentado resultados muito mais coerentes do que o procedimento mais simples que consistiria em determinar as direções defluxo a partir do MDET (ou MNT: modelo numérico do terreno) reamostrado para a baixa resolução, através de algoritmos como o da rotina MNTAlta4A. No procedimento de upscaling, geram-se direções de fluxo de baixa resolução que se apresentam em boa concordância com as informações de mais alta definição, que são as direções de fluxo e as áreas de drenagem de alta resolução. Os softwares comerciais de geoprocessamento mais conhecidos (Idrisi, ArcView, ArcGis) não realizam operações desse tipo. Paz et al. (2006) desenvolveu um algoritmo de upscaling de direções de fluxo baseado no trabalho de Reed (2003), o qual está implementado na rotina ../softs/Paz2006_routines/DirFluxo4/DirFluxo4.exe.
O algoritmo da rotina DirFluxo4.exe consiste basicamente em identificar um pixel exutório para cada célula, seguir o caminho de fluxo de pixel em pixel a partir do pixel exutório da célula e conforme esse caminho traçado determina para qual célula vizinha a célula analisada drena. O referido algoritmo pode ser resumidamente descrito em três etapas:
- Determinação do pixel exutório de cada célula: para uma determinada célula, escolhe-se inicialmente como pixel exutório aquele que apresenta a maior área de drenagem acumulada dentre todos os pixels contidos na célula; verifica-se o comprimento do curso d’água principal a montante desse pixel dentro da célula; caso esse comprimento seja superior a um valor mínimo pré-definido, o pixel testado é aceito como pixel exutório; tal valor mínimo corresponde ao parâmetro denominado de Caminho Mínimo de Montante ou CMM, cujo valor recomendado é igual a 1/5 da dimensão da célula; caso não seja atendido o critério do CMM, verifica-se se o pixel testado é o que drena a maior porção da célula e, em caso positivo, tal pixel é aceito para pixel exutório; caso contrário, escolhe-se novo pixel para ser testado dentre os demais de acordo com a maior área de drenagem acumulada e repetem-se as verificações subsequentes;
- Atribuição das direções de fluxo: a atribuição da direção de fluxo para cada uma das células é realizada percorrendo-se o caminho do escoamento desde seu pixel exutório; a cadapixel exutório de uma célula vizinha encontrado, verifica-se o incremento naárea de drenagem; caso seja superior a um valor mínimo pré-definido, a célula analisada drena para essa célula vizinha; tal valor mínimo constitui um parâmetro denominado Área Incremental Mínima ou AIM, e tem valor recomendado igual à área da célula; caso não atenda ao critério da AIM, continua-se a percorrer o caminho do escoamento, até encontrar o pixel exutório de uma célula vizinha que satisfaça o critério mencionado ou que saia da vizinhança; nesse último caso, atribui-se a direção para a última célula visitada; situações particulares são tratadas especificamente, como descrito em Paz et al. (2006);
- Correção de cruzamentos: esporadicamente podem ocorrer cruzamentos entre direções de fluxo de duas células, o que é desfeito com a correção da direção da célula cujo pixel exutório apresenta a menor área de drenagem acumulada dentre as duas célulasenvolvidas.
1.3. Preparação
1.3.1. Configurações básicas
Vamos carregar os pacotes necessários e definir algumas variáveis com a lozalização dos diretórios e dados, similar ao realizado nas etapas anteriores.
## carregando pacotes necessários
pcks <- c("raster", "rgdal", "RSAGA", "adehabitat", "maptools", "sp", "fields", "sfsmisc", "doBy", "stringr", "lattice", "rasterVis", "knitr")
## looping no nomes dos pacotes para carregá-los
l_load_pcks <- sapply(X = pcks,
FUN = require,
character.only = T)
l_load_pcks raster rgdal RSAGA adehabitat maptools sp
TRUE TRUE TRUE TRUE TRUE TRUE
fields sfsmisc doBy stringr lattice rasterVis
TRUE TRUE TRUE TRUE TRUE TRUE
knitr
TRUE Definindo opções, nome do usuário, diretório de trabalho (WD) e carregando funções armazenadas em scripts.
## configurando opção para mostrar barra de progresso durante operações com raster
rasterOptions(progress = "text")
## defina seu nome de usuario que definira o caminho até o diretório scripts
#user_name <- "hidrometeorologista"
user_name <- "hidro2roilan"
## 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"## carregando scripts necessários
## script com função para ajust das extents dos raster
## de forma a linhas e colunas serem múltiplos de 10
source("recode_flowdir.R")
## script para gerar plot de drenagem
source("flow_lines.R")
## script para gerar plot de drenagem
source("gdal_polygonizeR.R")
## script verificação das direções de fluxo pós-upscaling
source("check_flowdir.R")
source("show_flowdir_convention.R")
source("plot_zoom_cell.R")
source("plot_zoom_cell_v2.R")
## script para escrita no formato ascci DBHM
source("raster2ascii_dbhm.R")Vamos definir a bacia de interesse através da variável basin_name. A partir dela defineremos todos diretórios de acesso aos dados da bacia. O nome da bacia é usado para carregar o MDET de alta resolução da bacia na etapa anterior. Se alterar a basin_name altere o parâmetro fator_escala (definido adiante) de acordo com o salto do upscaling.
## defina o nome de sua bacia
# basin_name <- "mogu"
# basin_name <- "vacacai"
basin_name <- "pauliceia"
## diretório da bacia dentro da base de dados
database_basin_dir <- gsub("basin", basin_name, "../data_base/basin/")
database_basin_dir[1] "../data_base/pauliceia/"## mdet da bacia de alta resolução (100 m) gerado na etapa anterior
mdet_alta_file <- paste0(database_basin_dir,
gsub(pattern = "basin",
replacement = basin_name,
x = "mdet_basin_laea_int_crop.rst"
)
)
mdet_alta_file[1] "../data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst"mdet_alta <- raster(mdet_alta_file)
mdet_altaclass : RasterLayer
dimensions : 70, 60, 4200 (nrow, ncol, ncell)
resolution : 100, 100 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst
names : mdet_pauliceia_laea_int_crop
values : 572, 751 (min, max)## definir projeção cartográfica
## le raster tif que contem
mdet_basin_laea_file <- paste0(database_basin_dir, gsub("basin", basin_name, "mdet_basin_laea_int.tif"))
laea <- projection(raster(mdet_basin_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"projection(mdet_alta) <- laea
## pojeção geográfica
mdet_basin_file <- paste0(database_basin_dir,gsub("basin",basin_name,"mdet_basin.tif"))
mdet_basin_file[1] "../data_base/pauliceia/mdet_pauliceia.tif"ll <- projection(raster(mdet_basin_file))
ll[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"Vamos começar criando o sub-diretório da bacia em ../softs/Paz2006_routines/DirFluxo4/ e copiando o executável da rotina para este novo diretório.
## caminho para o diretório DirFluxo4
dirflux4_dir <- "../softs/Paz2006_routines/DirFluxo4"
## se não existir o diretório dirflux4_dir ele é criado
if(!"DirFluxo4" %in% dir("../softs/Paz2006_routines")) dir.create(dirflux4_basin_dir)
## caminho para o diretório DirFluxo4/bacin_name
dirflux4_basin_dir <- paste0(dirflux4_dir, "/", basin_name)
dirflux4_basin_dir[1] "../softs/Paz2006_routines/DirFluxo4/pauliceia"## se o diretório dirflux4_basin_dir não existir ele é criado
if(!basin_name %in% dir(dirflux4_dir)) dir.create(dirflux4_basin_dir)
## copiando executável para o dirflux4_basin_dir
file.copy(from = paste0(dirflux4_dir, "/DirFluxo4.exe"),
to = dirflux4_basin_dir)[1] TRUEVamos copiar os dois rasters gerados na etapa anterior:
- direção de fluxo de alta resolução
../softs/Paz2006_routines/MNTAlta/basin/Dir.* - área de drenagem acumulada de alta resolução
../softs/TauDEM-master/basin/area_taudem_km2.*
para o diretório ../softs/Paz2006_routines/DirFluxo4/basin/.
mntalta_basin_dir <- gsub("basin", basin_name, "../softs/Paz2006_routines/MNTAlta/basin")
## copiando arquivo de direções de fluxo
from <- paste0(mntalta_basin_dir, "/Dir.* ")
to <- dirflux4_basin_dir
cmd <- paste0("cp " , from, to)
cmd[1] "cp ../softs/Paz2006_routines/MNTAlta/pauliceia/Dir.* ../softs/Paz2006_routines/DirFluxo4/pauliceia"## executa comando no terminal linux
system(cmd)
## copiando arquivo de área de drenagem
taudem_basin_dir <- gsub("basin", basin_name, "../softs/TauDEM-master/basin")
from <- paste0(taudem_basin_dir,"/area_taudem_alta_km2.* ")
to <- dirflux4_basin_dir
cmd <- paste0("cp " , from, to)
cmd[1] "cp ../softs/TauDEM-master/pauliceia/area_taudem_alta_km2.* ../softs/Paz2006_routines/DirFluxo4/pauliceia"## executa comando no terminal linux
system(cmd)Esses arquivos raster devem ser renomeados para o padrão do nome dos arquivos de entrada da rotina DirFluxo4.exe, i.e., DirAlta.* e AreaAlta.*.
## renomear arquivos raster (*.rst) com extent ajustada para nomes dos arquivos de entrada para rotina
file.copy(from = paste0(dirflux4_basin_dir,c("/Dir.rst","/area_taudem_alta_km2.rst")),
to = paste0(dirflux4_basin_dir,c("/DirAlta.rst","/AreaAlta.rst")),
overwrite = T)[1] TRUE TRUE## renomear arquivos de documentacao (*.rdc)
file.copy(from = paste0(dirflux4_basin_dir,c("/Dir.rdc","/area_taudem_alta_km2.rdc")),
to = paste0(dirflux4_basin_dir,c("/DirAlta.rdc","/AreaAlta.rdc")),
overwrite = T)[1] TRUE TRUE1.3.2. Parâmetros para execução da rotina DirFluxo4
Durante a execução da rotina DirFluxo4.exe são solicitadas algumas informações sobre os atributos espaciais dos rasters de alta e baixa resolução. As expressões abaixo são para obter essas informações. Uma variável que deve ser definida é a fator_escala que indica o fator de escala correpondente ao salto da alta para baixa resolução. Por exemplo para bacia mogu esse salto é de 10 vezes, ou seja de 100 m para 1 km.
Entre as informações estão:
- a resolução alta (em graus)
## DEFINA O FATOR DE ESCALA PARA O UPSCALING
#fator_escala <- 10 # mogu e vacacai
fator_escala <- 5 # Pauliceia
## raster de alta resolução
ralta <- raster(paste0(dirflux4_basin_dir, "/AreaAlta.rst"))
## projetando de laea para lonlat
projection(ralta) <- laea
## raster de baixa resolução (10 x a res alta)
rbaixa <- aggregate(ralta, fact = fator_escala)
## define valores como os valores de cada célula
rbaixa <- setValues(rbaixa, 1:ncell(rbaixa))
rbaixaclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
data source : in memory
names : AreaAlta
values : 1, 168 (min, max)## escreve no formato idrisi
writeRaster(x = rbaixa,
filename = paste0(dirflux4_basin_dir,"/RASTERBAIXA.rst"),
NAflag = -9999,
format='IDRISI',
overwrite = T,
datatype = "INT2S")class : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/Paz2006_routines/DirFluxo4/pauliceia/RASTERBAIXA.rdc
names : layer
values : 1, 168 (min, max)## projeção do domínio em latlon
ralta_ll <- projectExtent(ralta, crs = ll)
## estimando a resolução usando a aproximação de que 1° equivale a 100 km no equador
one_deg <- 100 # km
## alta resolução
res_alta_ll <- unique(res(ralta))/(one_deg*10^3)
res_alta_ll[1] 0.001resolucao_alta <- sprintf("%f km",res_alta_ll)
resolucao_alta[1] "0.001000 km"res(ralta_ll) <- res_alta_ll
ralta_llclass : RasterLayer
dimensions : 63, 58, 3654 (nrow, ncol, ncell)
resolution : 0.001, 0.001 (x, y)
extent : -47.65808, -47.60008, -21.66993, -21.60693 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 - a resolução baixa (em graus) ou seja a resolução após upscaling;
res_baixa_ll <- res_alta_ll * unique(res(rbaixa)/res(ralta))
res_baixa_ll[1] 0.005- Área Incremental Mínima (Km²)
aim <- prod(res(rbaixa)/1000)
aim[1] 0.25- Caminho mínimo à montante (Km)
cmm <- 1/5 * (unique(res(ralta))/1000)
cmm[1] 0.021.4. Operação
A rotina a ser executada é a DirFluxo4.exe, disponível em ../softs/Paz2006_routines/DirFluxo4. Nesse diretório deve estar a pasta com o nome da bacia (i.e mogu ou vacacai). Os arquivos de entrada são aqueles obtidos após ajuste das linhas e colunas (seção 7.2) e que foram renomeados para DirAlta.rst (formato Idrisi/integer/binary) e DirAlta.rdc; e os arquivos para AreaAlta.rst (formato Idrisi/real/binary) e AreaAlta.rdc, que correspondem respectivamente às direções de fluxo e área de drenagem acumuladas de alta resolução. No diretório da sub-bacia deve estar o executável DirFluxo4.exe, que pode ser copiado do diretório precedente.
Ao executar a rotina DirFluxo4.exe, deve-se informar os valores em graus da baixa (0.01 °) e da alta (0.001 °) resolução, e depois checar se o número de linhas e colunas calculados para a baixa resolução apresentado na tela de saída da rotina está correto. Caso esteja correto deve-se prosseguir a execução da rotina, e informar se deseja ou não que seja lido arquivo referente à máscara de baixa resolução. Em nossos exemplos não usamos máscaras, portanto digite 0.
O passo seguinte é informar os valores dos parâmetros Área Incremental Mínima (AIM) e Caminho Mínimo de Montante (CMM). O valor usual de AIM é igual à área da célula (por exemplo, 1 km² para células de 1 x 1 km (~0.01 x 0.01°). Para o parâmetro CMM, o valor recomendado é igual a 1/5 da dimensão da célula (por exemplo, 0.2 km para célula de 0.01 x 0.01º). Atenção para as unidades dos parâmetros AIM (km²) e CMM (km). Os valores do parâmetro AIM e principalmente do parâmetro CMM não foram exaustivamente analisados quanto à influência sobre o processo de upscaling de direções de fluxo.
A principal saída da rotina DirFluxo4.exe é o arquivo raster DIRBAIXA.rst (formato Idrisi/integer/binary), contendo as direções de fluxo dos elementos de baixa resolução (células), seguindo a mesma codificação usada pela rotina MNTAlta4A.exe.
O arquivo de documentação correspondente pode ser gerado a partir do raster RASTERBAIXA.rst gerado acima.
Após é só renomear o arquivo RASTERBAIXA.rdc para DIRBAIXA.rdc, abra-o no Idrisi e atualize a resolução e os valores mínimoe máximo da imagem usando o menu TOOLS. A imagem DirBaixa.rst deve ter aspecto semelhante à imagem DirAlta.rst, mas com elementos de maiores dimensões devido à diferença de resolução.
Dois arquivos de saída da rotina DirFluxo4.exe dizem respeito aos pixels exutórios determinados pelo algoritmo. O arquivo pixelexu.dat (em formato ascii) contém a linha e a coluna do pixel exutório de cada célula, e serve como informação de entrada para a rotina Trechos.exe que será futuramente aplicada. O outro arquivo é denominado PixelExu.rst (formato Idrisi/integer/binary), e trata-se deuma imagem raster de alta resolução onde todos os pixels têm valor 0 exceto aqueles selecionados como pixel exutório de alguma célula, que têm valor 1. Essa imagem é apenas informativa e pode ser usada para visualizar quais os pixels exutórios escolhidos pelo algoritmo, e entender como foi determinada a direção de fluxo de cada célula. É preciso gerar o arquivo de documentação PixelExu.rdc, o que pode ser feito copiando o arquivo DirAlta.rdc e renomeando essa cópia. Em seguida deve-se abrir o arquivo PixelExu.rdc no Idrisi e atualizar a resolução e os valores mínimo e máximo pelo comando TOOLS.
1.5. Verificação da direção de fluxo do upscaling
- imagem raster com as direções de fluxo de baixa resolução (arquivos DirBaixa.rst e DirBaixa.rdc) (formato Idrisi/integer/binary);
- imagem raster de alta resolução (formato Idrisi/integer/binary), onde os pixels selecionados como pixel exutório de alguma célula, tem valor 1 e os demais têm valor 0 (arquivos PixelExu.rst e PixelExu.rdc).
Importando arquivos de saída da rotina DirFluxo4.exe.
## nome completo do arquivo DIRBAIXA.rst
dir_baixa_file <- paste0(dirflux4_basin_dir, "/DIRBAIXA.rst")
dir_baixa <- raster(dir_baixa_file)
projection(dir_baixa) <- laea
dir_baixaclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/Paz2006_routines/DirFluxo4/pauliceia/DIRBAIXA.rst
names : DIRBAIXA
values : 1, 128 (min, max)## extraindo valores das células
db_vals <- getValues(dir_baixa)
## valores únicos de dif de fluxo
unique(db_vals)[1] 64 128 32 16 2 1 4 8## qtos NAs?
sum(is.na(db_vals))[1] 0plot(dir_baixa)
1.6. Correção das direções de fluxo de baixa resolução
Após a execução da rotina de upscaling de direção de fluxo é preciso corrigir alguns casos em que a direção de fluxos foi incorretamente definida.
Vamos carregar os rasters de saída da rotina no R:
## importando area acum. de alta resolução
area_alta_file <- paste0(dirflux4_basin_dir, "/AreaAlta.rst")
area_alta <- raster(area_alta_file)
area_altaclass : RasterLayer
dimensions : 70, 60, 4200 (nrow, ncol, ncell)
resolution : 100, 100 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/Paz2006_routines/DirFluxo4/pauliceia/AreaAlta.rst
names : AreaAlta
values : 0.01, 25.02 (min, max)projection(area_alta) <- laea
## importando localização dos pixels exutório de alta resoliução
pixelexu <- raster(gsub("AreaAlta.rst","pixelexu.rst", area_alta_file))
projection(pixelexu) <- laea
## NA no IDRISI foi definido como 1
## alterando tal que exutório = 1 e não exutório = NA
pixelexu[is.na(pixelexu)] <- 1
pixelexu[pixelexu != 1] <- NA
pixelexuclass : RasterLayer
dimensions : 70, 60, 4200 (nrow, ncol, ncell)
resolution : 100, 100 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
data source : in memory
names : pixelexu
values : 1, 1 (min, max)plot(pixelexu,
col = tim.colors(2),
main = "Pixels exutórios das células da grade do MDET")
Vamos então aplicar uma função para identificar quais direções de fluxo foram definidas incorretamente. O objeto de saída dessa função é um objeto da classe list no R e possui dois elementos acessíveis pelo operador $. O primeiro (fd_check$mat) é uma matriz com valor 0 para células sem problema e 1 para células com problema. O segundo elemento (fd_check$tab) é uma tabela com a linha e coluna das células com problema.
## script para plotar zoom em torno do ponto com problema
fd_check <- check_flowdir(d = dir_baixa)
## pontos a serem corrigidos
fd_check$tab[1] lin col
<0 rows> (or 0-length row.names)Para correção das direções inconsistentes precisamos visualizar a situação de cada célula. A função plot_zoom_cell foi desenvolvida para gerar um gráfico com zoom na célula de interesse contendo os rasters relevantes para a correção da direção de fluxo da célula. A função show_flowdir_convention mostra em uma nova tela gráfica a convenção de direções de fluxo adotada.
O mais comum é a ocorrência de cruzamentos entre direções de fluxo de duas células, o que deve ser em geral desfeito pela correção da direção da célula cujo pixel exutório apresenta a menor área de drenagem acumulada dentre as duas células envolvidas.
Vamos ver o exemplo da primeira célula com problema na bacia. O terceiro gráfico (à direita) mostra o valor da área de drenagem acumulada do pixel exutório (de alta resolução) de cada célula.
cell = 78 
Direcções de fluxo de baixa resoluçao concertadas
cell = 63 
Após a análise das direções do entorno da célula, corrigimos a direção atribuindo um novo valor de direção de fluxo para célula com problema.
#icell <- cellFromRowCol(object = dir_baixa,
# rownr = fd_check$tab[i, "lin"],
# colnr = fd_check$tab[i, "col"])
desloca_lin <- 0
desloca_col <- 0
icell <- cellFromRowCol(dir_baixa,
rownr = fd_check$tab[i, "lin"] + desloca_lin,
colnr = fd_check$tab[i, "col"] + desloca_col)
icell
dir_baixa[icell]
dir_baixa[icell] <- 128 # para fora do raster (para cima, caso mogu)
dir_baixa[icell] Rodando novamente a função check_flowdir deve sumir a linha da tabela correspondente a célula corrigida anteriormente.
## rodamos novamente
fd_check <- check_flowdir(d = dir_baixa)
## pontos a serem corrigidos
fd_check$tabNo meu caso da bacia Vacacaí parece mais conveniente alterar a direção da célula localizada a noroeste da primeira célula identificada (por isso usamos as variáveis desloca_lin = -1 e desloca_col = -1). Vamos defini-la com direção para fora do raster (32, para oeste). Para os casos em que a célula com cruzamento está fora da área da bacia o critério da célula de menor área de drenagem não precisa ser seguido a risca, pois não influenciará a rede drenagem da bacia. Esse procedimento deve ser repetido para todas linhas da tabela fd_check$tab. Maior atenção deve ser dada às células com cruzamento dentro da bacia.
## fechando telas gráficas
dev.off(); dev.off()
## pegando como exemplo o primeiro caso
i <- 1
fd_check$tab[i, "col"]
fd_check$tab[i, "lin"]
## plot na região da célula correspondente a linha e coluna
plot_zoom_cell(fdbaixa = dir_baixa,
aaalta = area_alta,
exutalta = pixelexu,
icol = fd_check$tab[i, "col"],
irow = fd_check$tab[i, "lin"],
n = 10,
newindow = T,
winwidth = 18,
winheight = 6)
## gráfico com a convenção de direção de fluxo
show_flowdir_convention(newindow = T)
## para fechar os gráficos gerados via linha de comando
#
## após analise
desloca_lin <- 0
desloca_col <- 0
icell <- cellFromRowCol(dir_baixa,
rownr = fd_check$tab[i, "lin"] + desloca_lin,
colnr = fd_check$tab[i, "col"] + desloca_col)
icell
## valor de dreição atual
dir_baixa[icell]
## novo valor de direção definido conforme a situação de cada bacia
dir_baixa[icell] <- 128
dir_baixa[icell]
## rodamos novamente check_flowdir
fd_check <- check_flowdir(d = dir_baixa)
## pontos a serem corrigidos
fd_check$tabdir_baixa_file_corr <- gsub("DIRBAIXA", "DIRBAIXA_corr", dir_baixa_file)
writeRaster(x = dir_baixa,
filename = dir_baixa_file_corr,
NAflag = -9999,
format='IDRISI',
overwrite = T,
datatype = "INT2S")class : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/Paz2006_routines/DirFluxo4/pauliceia/DIRBAIXA_corr.rdc
names : layer
values : 1, 128 (min, max)2. Fração de área de cobertura da bacia
Vamos gerar um arquivo raster de baixa resolução que forneça a porcentagem de cada célula dentro da bacia, onde a bacia é definida pelo polígono com contorno da bacia gerado a partir dos rasters de alta resolução. Então precisamos do polígono da bacia de alta resolução e um raster na baixa resolução, por exemplo o dir_baixa. O polígono da bacia será rasterizado para baixa resolução, dividindo a célula de baixa resolução em 100 x 100 pixels e quantificando a porcentagem da célula contida dentro do polígono. O raster com fração de área em porcentagem dentro da bacia é então escrito no formato .asc.
## carregando shapefile de contorno da bacia
shape_bacia_alta <- gsub("basin",basin_name,"../softs/TauDEM-master/basin/mascara_bacia_alta_pol.shp")
mascara_bacia_alta_pol <- shapefile(shape_bacia_alta)
## definindo projeção do raster de direções de fluxo de baixa resolução
dir_baixaclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
data source : in memory
names : DIRBAIXA
values : 1, 128 (min, max)projection(dir_baixa) <- laea
## rasterizando polígono da bacia gerado com raster de alta resolução
rasterOptions(progress = "") # não mostrar barra de progresso
frac <- rasterize(x = mascara_bacia_alta_pol,
y = dir_baixa,
getCover = TRUE,
silent = TRUE)
projection(frac) <- laea
rasterOptions(progress = "text") # configurando para mostrar barra de progresso
## nome para o arquivo frac
frac_file <- gsub("basin", basin_name, "../frac_basin.asc")
## salvando fração de área dentro da bacia em ascii
raster2ascii_dbhm(r = frac, arqname = frac_file)
plot(frac, main = "Fração de área dentro da bacia (%)")
Para usar o raster frac como uma máscara para a ser aplicada sobre os próximos rasters de entrada que serão gerados vamos definir onde a fração de área da bacia for zero %, como NA. Esse novo raster será definido como fracNA.
fracNA <- frac
fracNA[fracNA <= 0 ] <- NA
plot(fracNA, main = "Fração de área dentro da bacia (%) \n Fora da bacia = NA")
Vamos então gerar o raster de direção de fluxo para a bacia, isto é, um raster contendo valores de direção de fluxo somente para as células que estão completa ou parcialmente dentro da bacia (frac > 0).
## direção de fluxo para a bacia
dir_baixa_mask <- mask(dir_baixa, fracNA)
## converte dir_baixa_mask para a convenção do arcgis requerida pelo DBHM
fdir <- recode_flowdir(dir_baixa_mask,
from = "iph",
to = "arcgis")
plot(fdir, main = "Direção de fluxo (ArcGIS)")
## nome completo para o arquivo dir.asc, salvando em:
## /home/user/DBHM/data_process/geo_data
dir_file_ascii <- gsub("basin", basin_name, "../dir_basin.asc")
## gera arquivo ascii
raster2ascii_dbhm(r = fdir, arqname = dir_file_ascii)3. Área acumulada de baixa resolução
A derivação da área de drenagem acumulada será novamente realizada com a rotina aread8 do TauDEM. No entanto vamos utilizar o raster de direção de fluxo de baixa resolução pós-correções (arquivo ../softs/Paz2006_routines/DirFluxo4/vacacai/DIRBAIXA_corr.rst).
Vamos então converter o raster de direção de fluxo de baixa resolução para a convenção usada pelo TauDEM e exportá-los para o formato .tif e .rst. O formato idrisi também é gerado porque eventualmente a rotina aread8 pode não funcionar como esperado com arquivo .tif gerado com a função writeGDAL. Isso ocorreu raramente, mas caso ocorra use o raster no formato Idrisi para gerá-lo via Idrisi. Esses arquivos são salvos no diretório ../softs/Paz2006_routines/DirFluxo4/basin.
## lê raster de dir fluxo corrigida
dir_baixa <- raster(dir_baixa_file_corr)
plot(dir_baixa, main = "Direção de fluxo (IPH) corrigida")
## reclassifica dir. fluxo da convenção do IPH para a do TauDEM
dir_baixa_taudem <- recode_flowdir(dir_baixa, from = "iph", to = "taudem")
plot(dir_baixa_taudem, main = "Direção de fluxo (TauDEM)")
## gera nome do arquivo de dir. fluxo obtido com o MNTAlta em formato tif requerido pelo taudem
dir_baixa_taudem_file <- gsub("DIRBAIXA.rst", "dir_baixa_taudem.tif", dir_baixa_file)
dir_baixa_taudem_file[1] "../softs/Paz2006_routines/DirFluxo4/pauliceia/dir_baixa_taudem.tif"## fdiralta_taudem em dormato SpatialGridDataFrame (pacote sp)
## para escrita em formato requerido pelo taudem
dir_baixa_taudem_sgdf <- as(dir_baixa_taudem, "SpatialGridDataFrame")
## salvar raster de dir. fluxo na convenção do taudem em formato tif (requerido pelo TauDEM)
writeGDAL(dataset = dir_baixa_taudem_sgdf,
fname = dir_baixa_taudem_file,
drivername = "GTiff",
type="Int16",
mvFlag = -32768)
## salvar raster de dir. fluxo na convenção do taudem em formato .rst
writeRaster(x = dir_baixa_taudem,
filename = gsub("\\.tif","\\.rst",dir_baixa_taudem_file),
NAflag = -9999,
format='IDRISI',
overwrite = TRUE,
datatype = "INT2S")class : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/Paz2006_routines/DirFluxo4/pauliceia/dir_baixa_taudem.rdc
names : layer
values : 1, 8 (min, max)A rotina aread8 do TauDEM será aplicada usando o raster de direções de fluxo de baixa resolução ../softs/Paz2006_routines/DirFluxo4/basin/dir_baixa_taudem.tif para obter o raster de área de drenagem acumulada (arquivo que será nomeado como ../softs/TauDEM-master/basin/area_baixa_taudem.tif).
## Area acumulada com TauDEM
area_baixa_taudem_file <- "../softs/TauDEM-master/basin/area_baixa_taudem.tif"
# dir_baixa_taudem_file <- gsub("dir_baixa_taudem","dir_baixa_taudem_idrisi",dir_baixa_taudem_file)
cmd <- paste0("mpiexec -n 4 aread8 -p ",
dir_baixa_taudem_file,
" -ad8 ",
area_baixa_taudem_file)
## atualizando comando de acordo com o sub-diretório da bacia
cmd <- gsub(pattern = "basin",
replacement = basin_name,
x = cmd)
system(cmd)
## importando raster de area acum de baixa resolução
area_baixa_taudem <- raster(gsub(pattern = "basin",
replacement = basin_name,
x = area_baixa_taudem_file)
)
## atualizando valores mínimos e máximos do raster
area_baixa_taudem <- setMinMax(area_baixa_taudem)
|
| | 0%area_baixa_taudemclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/TauDEM-master/pauliceia/area_baixa_taudem.tif
names : area_baixa_taudem
values : 1, 45 (min, max)Para deixar o raster area_baixa_taudem compatível com a área de drenagem do ArcGIS (padrão requerido pelo DBHM) precisamos descontar 1 célula do total de células drenadas até cada ponto. Dessa forma células de cabeceira terão valor 0, ou seja, nenhuma célula drena para elas.
## ajuste do area para manter consistência com a area acum gerada com o arcgis
area_baixa_taudem <- area_baixa_taudem - 1
area_baixa_taudemclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : area_baixa_taudem
values : 0, 44 (min, max)Vamos visualizar a área de drenagem:
## raster de área de drenagem taudem
plot(area_baixa_taudem,
col = tim.colors(32),
main = "Área de drenagem de baixa res. (n° de células)")
plot(mascara_bacia_alta_pol , add = TRUE)
shapefile_exut_baixa_file <- gsub(pattern = "basin",
replacement = basin_name,
x = paste0(dirname(area_baixa_taudem_file),"/exutorio_baixa.shp"))
sp_xy_exut_baixa <- shapefile(x = shapefile_exut_baixa_file)
plot(sp_xy_exut_baixa, add = T, lty = 6, lwd = 3, col = 2)
# shapefile_exut_baixa_file2 <- "../data_base/data_pauliceia/exutorio_baixa.shp"
# sp_xy_exut_baixa2 <- shapefile(x = shapefile_exut_baixa_file2)
# plot(sp_xy_exut_baixa2, add = T, lty = 6, lwd = 3, col = 3)Dar um (ou mais) zoom(s) na região do ponto exutório, para encontrar a célula correspondente ao exutório da bacia.
sp_xy_exut <- shapefile(gsub("basin",basin_name,"../softs/TauDEM-master/basin/exutorio_alta.shp"))
sp_xy_exut
plot(sp_xy_exut , add = TRUE, cex = 2, col = 2)
## zoom na região do ponto exutório
zoom(area_baixa_taudem, col = tim.colors(32))
plot(mascara_bacia_alta_pol , add = TRUE)
plot(sp_xy_exut , add = TRUE, cex = 2, col = 2)Após definir as coordenadas do exutório vamos convertê-las para um objeto do tipo SpatialPoints (objeto do pacote sp, especípico para manipulação de dados spaciais) e salvá-lo num arquivo shapefile.
## clique no ponto exutório para extrair as coordenadas xy
xy_exut_baixa <- locator(n=1)
## converte coordenadas para SpatialPoints
sp_xy_exut_baixa <- SpatialPoints(xy_exut_baixa)
sp_xy_exut_baixa
## projeta para laea
proj4string(sp_xy_exut_baixa) <- laea
## nome do arquivo shapefile como localizacao do exutório
shapefile_exut_baixa_file <- gsub(pattern = "basin",
replacement = basin_name,
x = paste0(dirname(area_baixa_taudem_file),"/exutorio_baixa.shp")
)
## salva shapefile
shapefile(x = sp_xy_exut_baixa,
filename = shapefile_exut_baixa_file,
overwrite = TRUE)Vamos obter a célula do raster correspondente aquelas coordendas do exutório e obter a área de drenagm no exutório para compará-la com o valor de alta resolução.
## determina a célula do exutório
cell_exut_baixa <- cellFromXY(object = area_baixa_taudem,
xy = sp_xy_exut_baixa)
cell_exut_baixa[1] 113## área acumulada da bacia naquela célula
area_bacia_baixa <- raster::extract(area_baixa_taudem, cell_exut_baixa)
area_bacia_baixa * prod(res(rbaixa)/1000) # subestimativa?
11 Vamos gerar a área de drenagem acumulada só para bacia formada a partir do ponto exutório de baixa resolução, informando o shapefile recém gerado no comando de execução da rotina aread8.
## Area acumulada com TauDEM usando o shape do exutório
cmd <- paste0("mpiexec -n 4 aread8 -p ", dir_baixa_taudem_file," -ad8 ../softs/TauDEM-master/basin/area_taudem_baixa_exutorio.tif", " -o ", shapefile_exut_baixa_file)
## atualizando comando de acordo com o sub-diretório da bacia
cmd <- gsub(pattern = "basin",
replacement = basin_name,
x = cmd)
system(cmd)Vamos importar a área de drenagem de baixa resolução obtida com o ponto exutório e deixá-la compatível com a área de drenagem do ArcGIS (padrão requerido pelo DBHM), novamente precisamos descontar 1 célula do total de células drenadas até cada ponto. Dessa forma células de cabeceira terão valor 0, ou seja, nenhuma célula drena para elas. Após esse procedimento produziremos uma máscara da bacia de baixa resolução.
area_taudem_baixa_exutorio_file <- gsub(pattern = "basin",
replacement = basin_name,
x = "../softs/TauDEM-master/basin/area_taudem_baixa_exutorio.tif")
area_taudem_baixa_exutorio <- raster(area_taudem_baixa_exutorio_file)
## atualizando estatísticas do raster
area_taudem_baixa_exutorio <- setMinMax(area_taudem_baixa_exutorio)
|
| | 0%area_taudem_baixa_exutorioclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/softs/TauDEM-master/pauliceia/area_taudem_baixa_exutorio.tif
names : area_taudem_baixa_exutorio
values : 1, 45 (min, max)## descontando 1 célula para consistência com a área do arcgis
area_taudem_baixa_exutorio <- area_taudem_baixa_exutorio - 1
area_taudem_baixa_exutorioclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : area_taudem_baixa_exutorio
values : 0, 44 (min, max)plot(area_taudem_baixa_exutorio,
col = tim.colors(32),
main = "Área de drenagem (TauDEM) usando exutório")
## contorno de alta resolução
plot(mascara_bacia_alta_pol,
border = "red",
add = TRUE)
A partir da área acumulada de baixa resolução determinada a partir do ponto exutório vamos gerar uma máscara com valor 1 dentro da bacia e 0 fora.
## mascara de baixa resolução
mascara_baixa <- area_taudem_baixa_exutorio
mascara_baixa[!is.na(mascara_baixa)] <- 1
plot(mascara_baixa, main = "Máscara da bacia baseada na Área de drenagem")
## salvando em ascii mascara de baixa resolução
mascara_baixa_file <- paste0(dirname(area_taudem_baixa_exutorio_file), "/mascara_baixa.asc")
raster2ascii_dbhm(r = mascara_baixa,
arqname = mascara_baixa_file)Vamos salvar o raster de área de drenagem de baixa resolução para bacia.
## salvar raster da área de drenagem acumulada de baixa resolução (1000 m)
acc_file_ascii <- gsub("basin", basin_name, "../acc_basin.asc")
## a <- mask(acc, fracNA)
acc <- area_taudem_baixa_exutorio
plot(acc,
col = tim.colors(32),
main = "Área de drenagem de baixa resolução \n (n° de células)")
## células de cabeceira
plot(acc == 0, main = "Células de cabeceira")
#zoom(acc, col = tim.colors(32))
#flow_lines(dir_baixa, length = 0.001, col = "blue", lwd = 2)
## gera arquivo ascii
raster2ascii_dbhm(r = acc, arqname = acc_file_ascii)5. Estrutura de diretórios e dados gerados
A figura abaixo mostra a estrutura de diretórios e os arquivos gerados nessa etapa são destacados.
6. Lista de arquivos gerados
- dir.asc direções de fluxo da bacia (valores fora da bacia são definidos como -9999) expressa na convenção do ArcGIS (Figura 2);
fdirclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
data source : in memory
names : arcgis
values : 1, 128 (min, max)dir_file_ascii[1] "../dir_pauliceia.asc"- acc.asc área de drenagem acumulada da bacia expressa como número de células à montante que escoam para a célula, na projeção cartográfica LAEA.
acc_file_ascii[1] "../acc_pauliceia.asc"accclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : area_taudem_baixa_exutorio
values : 0, 44 (min, max)- frac.asc fração de área de cada célula dentro da bacia (delimitada pelo polígono da bacia derivada do MDET de alta resolução) expressa em porcentagem de 0 a 100%;
frac_file[1] "../frac_pauliceia.asc"fracclass : RasterLayer
dimensions : 14, 12, 168 (nrow, ncol, ncell)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs
data source : in memory
names : layer
values : 0, 100 (min, max)7. Referências Citadas
Paz, A.R. MGBgis - Modelo hidrológico distribuído MGB-IPH, ETAPA 1 de preparação das informações de entrada. Manual do Usuário, IPH-UFRGS, Porto Alegre (RS). 2008.
Paz, A. R., W. Collischonn, e A. L. L. Silveira (2006), Improvements in large scale drainage networks derived from digital elevation models, Water Resourc. Res., 42 (8), doi: 10.1029/2005WR004544.
Reed, S. M. (2003), Deriving flow directions for coarse-resolution (1-4 km) gridded hydrologic modeling, Water Resourc. Res., 39(9), 1238, doi:10.1029/2003WR001989.