DBHM-GIS: PARTE 4
17-11-2014
Universidade Federal de Santa Maria - UFSM
Departamento de Física
Prof. Jônatan Tatsch
Aluno Roilan Hernandez Valdes
1. Extração de comprimentos e declividades dos trechos de rio
1.1. Descrição
Os modelos hidrológicos distribuídos utilizam métodos de propagação do escoamento ao longo da rede de drenagem, isto é, entre as células do modelo. Esse métodos necessitam do comprimento e da declividade do trecho de rio associado à ligação entre duas células, ao longo do qual é feita a propagação.
Para a propagação do escoamento, a maior simplificação quanto ao comprimento do trecho de rio seria considerar todos os trechos com comprimento igual à dimensão da célula, quando o trecho é entre duas células localizadas ortogonalmente entre si, ou igual a 1.4142 vezes a dimensão da célula, no caso de células em diagonal. Entretanto, isso significa uma enorme perda na qualidade da representação do processo físico, que pode conduzir ao mau ajuste do modelo ou forçar o ajuste dos parâmetros de forma compensatória. Por outro lado, informações de comprimentos e declividades dos trechos de rio, principalmente em grandes bacias, raramente estão disponíveis e a extração manual de tais informações a partir de mapas impressos representa tarefa inviável face às dimensões das bacias.
Nós utilizaremos o algoritmo desenvolvido por Paz & Collischonn (2007) para extrair automaticamente os comprimentos e declividades dos trechos de rio. O algoritmo se baseia em informações de alta resolução geradas em etapas anteriores (MNT, MNT com depressões removidas, áreas de drenagem acumuladas) e foi desenvolvido baseado no algoritmo de upscaling de direções de fluxo.
1.2 Entradas e saídas
Os arquivos de entrada dessa rotina são:
DirAlta.r*: imagem raster com direções de fluxo de alta resolução, formato do Idrisi (integer/binary), gerada na etapa 5.1.2;
AreaAlta.r*: imagem raster com áreas acumuladas de drenagem de alta resolução, formato do Idrisi (real/binary), gerada na etapa 6.2.1;
MNTfill.r*: imagem raster com MNT de alta resolução com as depressões removidas (modificado pelo algoritmo de upscaling de direções de fluxo), formato do Idrisi (integer/binary), gerada na etapa 5.1.2;
pixelexu.dat: arquivo em formato ascii contendo informação do número da linha e da coluna do pixel exutório de cada célula, gerado na etapa 7.6;
DirBaixa_corr.r*: imagem raster com direções de fluxo de baixa resolução, formato do Idrisi (integer/binary), gerada na etapa 7.7, isto é, as direções de baixa corrigidas.
Os produtos principais são as imagens raster de baixa resolução contendo os comprimentos e as declividades dos trechos de rio associados a cada célula, os quais servirão efetivamente como entrada no DBHM:
trechosTOT.r*: imagem raster de baixa resolução (formato Idrisi/real/binary), onde a cada elemento (célula do modelo) está associado um valor referente ao comprimento do trecho de rio, em km;
decliv.r*: imagem raster de baixa resolução (formato Idrisi/real/binary), onde a cada elemento (célula do modelo) está associado um valor referente à declividade do trecho de rio, em m/km.
Para as células de cabeceira, a rotina atribui um valor simbólico de declividade igual a -0,99. À exceção de tais células, todas as demais devem ter declividade não nula. Isso pode e deve ser checado antes de usar a informação gerada. Recomenda-se gerar outra imagem fazendo uma reclassificação das declividades, de forma que todos os valores -0,99 se tornem 0. Nós faremos essa reclassificação no R.
Diversos arquivos complementares são gerados:
TrechosMONT.r*: imagem raster de baixa resolução (formato Idrisi/real/binary), onde a cada elemento está associado um valor referente ao comprimento do sub-trecho de rio de montante, em km;
TrechosJUS.r*: imagem raster de baixa resolução (formato Idrisi/real/binary), onde a cada elemento está associado um valor referente ao comprimento do sub-trecho de rio de jusante, em km (para cada elemento, a soma do sub-trecho de montante e do sub-trecho de jusante é igual ao comprimento total que consta no arquivo TrechosTOT.rst);
caminhoTOT.r*: imagem raster de alta resolução (formato Idrisi/integer/binary), onde os pixels assumem o valor 0 ou 1. Os pixels com valor 1 fazem parte do trecho de rio associado a alguma célula;
caminhoMONT.r*: imagem raster de alta resolução (formato Idrisi/integer/binary), onde os pixels assumem o valor 0 ou 1. Os pixels com valor 1 fazem parte do sub-trecho de rio de montante associado a alguma célula;
ptostrechos.r*: imagem raster de alta resolução (formato Idrisi/integer/binary), onde os pixels assumem o valor 0 ou 1. Os pixels com valor são aqueles correspondentes ao ponto inicial ou final do trecho de rio associado a alguma célula;
Novos arquivos raster adicionados ao código Fortran Trechos5b_jdt.f:
cotaMon.r*: imagem raster de baixa resolução (formato Idrisi/real/binary) com valores da maior altitude dos pixels a montante do ponto exutório da célula;
cotaJus.r*: imagem raster de baixa resolução (formato Idrisi/real/binary) com valores da menor altitude dos pixels a jusante do ponto exutório da célula;
1.3 Preparação
1.3.1. Configurações básicas
Vamos carregar os pacotes necessários e definir algumas variáveis com a localizaçã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", "doParallel")
## 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 doParallel
TRUE TRUE Definindo opções, nome do usuário, diretório de trabalho (WD) e carregando funções armazenadas em scripts.
NOTA: altere o parâmetro
usernamede acordo com a sua configuração
## 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 <- system("echo $USER" , intern = T )
## 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
source("flowLength.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.
NOTA: altere o parâmetro
basin_namede acordo com a bacia de estudo
## defina o nome de sua bacia
# basin_name <- "mogu"
# basin_name <- "vacacai"
basin_name <- "pauliceia"1.3.2. Dados de entrada e diretórios
Vamos começar criando as variáveis que definem o sub-diretório da bacia que deve existir nos diretórios de upscaling de direções de fluxo e de determinação do comprimento e declividade dos trechos de rio.
## caminho para o diretório DirFluxo4
dirflux4_dir <- "../softs/Paz2006_routines/DirFluxo4"
trechos5b_dir <- "../softs/Paz2006_routines/Trechos5b"
## se não existir o diretório trechos5b_dir ele é criado
if(!"Trechos5b" %in% dir("../softs/Paz2006_routines")) dir.create(trechos5b_dir)
## COPIAR O CÓDIGO FORTRAN PARA ESSE DIRETÓRIO
##
## caminho para o diretório DirFluxo4/basin_name
dirflux4_basin_dir <- paste0(dirflux4_dir, "/", basin_name)
dirflux4_basin_dir[1] "../softs/Paz2006_routines/DirFluxo4/pauliceia"## caminho para o diretório DirFluxo4/basin_name
trechos5b_basin_dir <- paste0(trechos5b_dir, "/", basin_name)
trechos5b_basin_dir[1] "../softs/Paz2006_routines/Trechos5b/pauliceia"## se o diretório dirflux4_basin_dir não existir ele é criado
if(!basin_name %in% dir(trechos5b_dir)) dir.create(trechos5b_basin_dir)
## copiando código fortran para o dirflux4_basin_dir
fortran_file <- "trechos5b_jdt.f90"
file.copy(from = paste0(trechos5b_dir, "/", fortran_file),
to = trechos5b_basin_dir)[1] TRUEVamos copiar os arquivos de entrada da rotina Trechos5b para o diretório ../softs/Paz2006_routines/Trechos5b/basin. Os arquivos podem ser copiados pelas expressões abaixo.
## localização dos arquivos de entrada da rotina que devem ser copiados
files_from <- list.files(path = dirflux4_basin_dir,
pattern = "Alta|BAIXA_corr|pixelexu.*[rst,rdc,dat]$",
full.names = T)
files_from[1] "../softs/Paz2006_routines/DirFluxo4/pauliceia/AreaAlta.rdc"
[2] "../softs/Paz2006_routines/DirFluxo4/pauliceia/AreaAlta.rst"
[3] "../softs/Paz2006_routines/DirFluxo4/pauliceia/DirAlta.rdc"
[4] "../softs/Paz2006_routines/DirFluxo4/pauliceia/DirAlta.rst"
[5] "../softs/Paz2006_routines/DirFluxo4/pauliceia/DIRBAIXA_corr.rdc"
[6] "../softs/Paz2006_routines/DirFluxo4/pauliceia/DIRBAIXA_corr.rst"
[7] "../softs/Paz2006_routines/DirFluxo4/pauliceia/pixelexu.dat"
[8] "../softs/Paz2006_routines/DirFluxo4/pauliceia/pixelexu.rdc"
[9] "../softs/Paz2006_routines/DirFluxo4/pauliceia/pixelexu.rst" ## localização do destino dos arquivos que devem ser copiados
files_to <- gsub("DirFluxo4","Trechos5b", files_from)
files_to <- gsub("DIRBAIXA_corr","DirBaixa", files_to)
files_to[1] "../softs/Paz2006_routines/Trechos5b/pauliceia/AreaAlta.rdc"
[2] "../softs/Paz2006_routines/Trechos5b/pauliceia/AreaAlta.rst"
[3] "../softs/Paz2006_routines/Trechos5b/pauliceia/DirAlta.rdc"
[4] "../softs/Paz2006_routines/Trechos5b/pauliceia/DirAlta.rst"
[5] "../softs/Paz2006_routines/Trechos5b/pauliceia/DirBaixa.rdc"
[6] "../softs/Paz2006_routines/Trechos5b/pauliceia/DirBaixa.rst"
[7] "../softs/Paz2006_routines/Trechos5b/pauliceia/pixelexu.dat"
[8] "../softs/Paz2006_routines/Trechos5b/pauliceia/pixelexu.rdc"
[9] "../softs/Paz2006_routines/Trechos5b/pauliceia/pixelexu.rst" file.copy(from = files_from, to = files_to, overwrite = T)[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEO raster MNTfill.rst gerado na DBHM-GIS Parte 2 também é uma entrada, logo deve ser copiado para o sub-diretório da bacia no diretório ../softs/Paz2006_routines/Trechos5b/basin.
mdet_alta_fill_file <- gsub(pattern = "basin",
replacement = basin_name,
x = "../softs/Paz2006_routines/MNTAlta/basin/MNTfill.rst")
## localização do destino dos arquivos que devem ser copiados
file_to <- gsub("MNTAlta","Trechos5b", mdet_alta_fill_file)
file_to <- gsub("basin",basin_name, file_to)
file_to[1] "../softs/Paz2006_routines/Trechos5b/pauliceia/MNTfill.rst"file.copy(from = mdet_alta_fill_file,
to = file_to,
overwrite = T)[1] TRUEfile.copy(from = gsub("\\.rst","\\.rdc", mdet_alta_fill_file),
to = gsub("\\.rst","\\.rdc", file_to),
overwrite = T)[1] TRUEOs comprimentos e declividades dos trechos são determinados com a rotina trechos5b_jdt.f, disponível no diretório ../softs/Paz2006_routines/Trechos5b. Durante sua execução devem ser informados alguns atributos espaciais do domínio da bacia em alta e baixa resolução, tais como:
- a resolução alta do MDET na projeção geográfica;
## obtendo string com parâmetros da projeção laea
database_basin_dir <- gsub("basin", basin_name, "../data_base/basin/")
database_basin_dir[1] "../data_base/pauliceia/"mdet_basin_laea_file <- paste0(database_basin_dir,
gsub("basin", basin_name, "mdet_basin_laea_int.tif")
)
mdet_basin_laea_file[1] "../data_base/pauliceia/mdet_pauliceia_laea_int.tif"mdet_alta <- raster(mdet_basin_laea_file)
mdet_altaclass : RasterLayer
dimensions : 79, 65, 5135 (nrow, ncol, ncell)
resolution : 100, 100 (x, y)
extent : -3101.52, 3398.48, -3822.855, 4077.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/data_base/pauliceia/mdet_pauliceia_laea_int.tif
names : mdet_pauliceia_laea_int
values : 566, 751 (min, max)laea <- projection(mdet_alta)
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"## obtendo string com parâmetros da 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"## nome completo do arquivo DIRBAIXA.rst
dir_baixa_file <- grep("DirBaixa.rst", files_to, value = TRUE)
dir_baixa_file[1] "../softs/Paz2006_routines/Trechos5b/pauliceia/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/Trechos5b/pauliceia/DirBaixa.rst
names : DirBaixa
values : 1, 128 (min, max)## projetando de laea para lonlat
## projeção do domínio em latlon
p_mdet_ll <- projectExtent(dir_baixa, crs = ll)
## resolução em latlon
# res(p_mdet_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(mdet_alta))/(one_deg*10^3)
## alta resolução
sprintf("%f",res_alta_ll)[1] "0.001000"- a resolução baixa (após upscaling) do MDET na projeção geográfica;
res(p_mdet_ll) <- res_alta_ll
p_mdet_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 ## resolução baixa, após upscaling
res_baixa_ll <- res_alta_ll * unique(res(dir_baixa)/res(mdet_alta))
sprintf("%f",res_baixa_ll)[1] "0.005000"- as coordenadas do canto superior direito do domínio espacial da bacia na projeção geográfica;
min_x <- xmin(p_mdet_ll)
min_x[1] -47.65808max_y <- ymax(p_mdet_ll)
max_y[1] -21.606931.4 Execução
Diferente dos procedimentos anteriores nós iremos compilar e executar o código fortran da rotina terchos5b_jdt.f. Abra um terminal linux (ctrl+alt+t) e digite:
$ cd ../softs/Paz2006_routines/Trechos5b/basin
$ ifort trechos5b_jdt.f90 -assume byterecl -o trechos5b_jdt; ./trechos5b_jdtComo ocorre nas rotinas anteriormente descritas, os arquivos de documentação (.rdc) correspondentes a cada imagem raster (.rst) devem ser criados manualmente. Isso é feito a partir de cópia de arquivo .rdc equivalente, renomeando-o, abrindo-o no Idrisi e atualizando o valor da resolução e do máximo e mínimo da imagem através da opção TOOLS/ Calculate Resolution e Calculate Min/Max.
Para criar um arquivo .rdc referente a uma imagem de baixa resolução, deve-se usar como base (copiar) um arquivo .rdc correspondente a uma imagem também de baixa resolução (p.ex.: ../softs/Paz2006_routines/DirFluxo4/pauliceia/RASTERBAIXA.rdc). E assim analogamente para gerar arquivos .rdc de imagens de alta resolução. Também deve-se ter atenção quanto ao tipo de dado da imagem, se inteiro ou real. Ao criar o arquivo .rdc e abri-lo no Idrisi, deve-se ajustar o formato manualmente antes de efetuar a atualização da resolução e dos valores mínimo e máximo.
Os arquivos cotaMon.rst e cotaJus.rst são gerados no formato real e contém a elevação do pixel à montante e à jusante do pixel exutório da célula de baixa resolução. Esses rasters serão usados posteriormente para derivação do MDET de baixa resolução.
1.5. Pós-Processamento dos rasters de declividade e comprimento dos trechos
Vamos copiar o arquivo .rdc que servirá como base para atualização dos rasters de baixa resolução através do Idrisi. A atualização dos rasters de alta resolução podem ser feitas da cópia de arquivos como o MNTFill.rst.
file.copy(from = paste0(dirflux4_basin_dir, "/RASTERBAIXA.rdc"),
to = paste0(trechos5b_basin_dir, "/RASTERBAIXA.rdc"),
overwrite = TRUE)[1] TRUEFaça atualização através do Idrisi dos arquivos raster:
- decliv.rst
- trechosTOT.rst
- cotaMon.rst
- cotaJus.rst
Vamos importar os dois principais rasters de saída da rotina trechos5b e ajustá-los para os padrões requeridos pelo DBHM.
## Declividade dos trechos de rio
slope_iph_file <- paste0(trechos5b_basin_dir, "/", "decliv.rst")
slope_iph <- raster(slope_iph_file)
slope_iphclass : 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/Trechos5b/pauliceia/decliv.rst
names : decliv
values : 0, 175.8582 (min, max)plot(slope_iph)
A declividade dos trechos de rio está em m/km. Valores negativos e muito elevados (superiores à 600, geralmente ocorrendo nas bordas do domínio) devem ser substituídos por NA. O raster de fração de área dentro da bacia (frac.asc) para aplicá-lo como máscara sobre o raster de declividade.
## nome do arquivo raster de fração de área dentro da bacia
frac_file <- gsub("basin",basin_name,"../frac_basin.asc")
frac <- raster(frac_file)
## raster frac com NA para células completamente fora da bacia
fracNA <- frac
fracNA[fracNA <= 0 ] <- NA
## por precaução eliminamos dados de declividade negativa
slope_iph[slope_iph < 0] <- NA
## aplicar máscara
slope_iph_mask <- mask(slope_iph, fracNA)
plot(slope_iph_mask,
col = tim.colors(32),
main = "Declividade (m/Km)")
Vamos converter a declividade para m/m e substituir valores inferiores a 0.0001 m/m (~10 cm/km) por 0.0001, para evitar problemas numéricos no modelo de propagação do DBHM.
## declividade em metros/metros
rise_run <- slope_iph_mask/1000
plot(rise_run,
col = tim.colors(32),
main = "Declividade (m/m)")
## limite inferior de declividade
slp_lim_inf <- 0.0001
## plot com localizações de baixa declividade
plot(rise_run <= slp_lim_inf,
main = paste0("Declividade < ", sprintf("%1.4f",slp_lim_inf)," m/m"))
rise_run[rise_run <= slp_lim_inf] <- slp_lim_inf
plot(rise_run,
col = tim.colors(32),
main = paste0("Declividade restrita ao \n mínimo de ", sprintf("%1.4f",slp_lim_inf), " m/m"))
A devidade em % é determinada multiplicando-se a declividade em m/m por 100. Esse raster de declividade (%) do leito do rio usado na rotina de propagação da água na rede de drenagem será salvo como slope_iph.asc.
## declividade em porcentagem
slope_iph_pctg <- round(rise_run*100, 6)
slope_iph_pctgclass : 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 : layer
values : 0.01, 10.09227 (min, max)## preparação de plot
## vetor de valores de slope_iph_pctg_mask
v <- getValues(slope_iph_pctg)
## posições das marcas na barra de cores
position_scale <- hist(v, plot=F)$breaks
## lista para legenda de cores para levelplot
myColorkey <- list(at = position_scale, ## where the colors change
labels = list(at = position_scale)
)
## plot usando função levelplot do pacote rasterVis
levelplot(slope_iph_pctg,
col.regions = rev(terrain.colors(length(position_scale)+1)),
colorkey = myColorkey,
main = "Declividade (%)"
)
## nome do arquivo de declividade
slope_iph_pctg_file <- gsub("basin", basin_name, "../slope_iph_basin.asc")
# writeRaster(slope_iph_pctg, filename = slope_ip_pctg_file, overwrite = T)
raster2ascii_dbhm(r = slope_iph_pctg, arqname = slope_iph_pctg_file)O ascii slope_iph.asc foi gerado pela função raster2ascii_dbhm disponível no script raster2ascii_dbhm.R. A declividade também é usada pelo modelo de superfície SiB2 porém sua unidade é em radianos. Essa declividade será salva como slope_sib.asc.
## seno da declividade para a bacia em radianos
slpsib_iph <- round(sin(atan(rise_run)), 6)
plot(slpsib_iph, main = "Declividade, sin(ang)")
slpsib_iph_file <- gsub("basin", basin_name, "../slpsib_iph_basin.asc")
#writeRaster(sina, filename = slope_sib_file, overwrite = T)
raster2ascii_dbhm(r = slpsib_iph, arqname = slpsib_iph_file)Vamos importar o raster de comprimento dos trechos de rios, converter para metros, restringir os valores a um limite inferior de 0.6*resolução espacial, e salvar como ascii.
reaches_file <- paste0(trechos5b_basin_dir, "/trechosTOT.rst")
reaches <- raster(reaches_file)
reachesclass : 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/Trechos5b/pauliceia/trechosTOT.rst
names : trechosTOT
values : 0.0995038, 1.298764 (min, max)## conversão para metros
reaches <- round(reaches * 1000, 2)
reachesclass : 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 : layer
values : 99.5, 1298.76 (min, max)reaches_mask <- mask(reaches, fracNA)
plot(reaches_mask, main = "Comprimento dos trechos de rios (m)")
#plot(reaches_mask)
## porcentagem da área da bacia com trechos inferiores ao limiar de 0.6*res()
## limiar de comprimento de trecho de rio
limiar_comp <- mean(res(reaches_mask))*0.6
limiar_comp[1] 300## histograma dos comprimentos de trechos de rio
histogram(reaches_mask,
main = "histograma de Comprimento de trechos",
xlab = "comprimento (m)",
ylab = "frequencia de ocorrência (%)")
## % da bacia com comprim. inferiores a limiar_comp
prop.table(table(reaches_mask[] < limiar_comp))*100
FALSE TRUE
71.875 28.125 ## substituindo valores inferiores a limiar_comp por limiar_comp
reaches_mask[reaches_mask < limiar_comp] <- limiar_comp
reaches_maskclass : 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 : layer
values : 300, 1298.76 (min, max)## preparação de plot
## vetor de valores de slope_iph_pctg_mask
v_reaches <- getValues(reaches_mask)
## posições das marcas na barra de cores
position_scale <- hist(v_reaches, plot=F)$breaks
## lista para legenda de cores para levelplot
myColorkey <- list(at = position_scale,
labels = list(at = position_scale))
## plot
levelplot(reaches_mask,
col.regions = rev(terrain.colors(length(position_scale)+1)),
colorkey = myColorkey
#main = "Comprimento de trechos limitados ao mínimo de 600 m"
)
## nome do arquivo raster de comprimento dos trechos
rivlen_sib_file <- gsub("basin", basin_name, "../rivlen_basin.asc")
# writeRaster(reaches, filename = rivlen_sib_file, overwrite = T)
raster2ascii_dbhm(r = reaches_mask, arqname = rivlen_sib_file)Vamos comparar o comprimento dos trechos de rio com a aproximação frequentemente usada que assume valores de 0.96194*resolução horizontal para células cardinais e 1.36039*resolução horizontal para células diagonais, conforme De Smith (2004). Para isso usamos o raster de direções de fluxo de baixa resolução na convenção do IPH para definir os trechos diagonais e cardinais.
## trechos de rio estimados por aproximação mais comum
trechos_baixa <- setRiverLength(r = dir_baixa,
lado = 0.96194,
diagonal = 1.36039)
## gráfico dos trechos assumindo apenas 2 valores
levelplot(trechos_baixa)
## diferença
plot(trechos_baixa - reaches_mask,
col = tim.colors(32),
main = "Diferença entre trechos de rio pelo\n método tradicional e via upscaling")
## histograma da diferença entre as declividades dos dois métodos
update(histogram(trechos_baixa - reaches_mask),
scales = list(cex = 1.2),
main = "histograma da diferença entre os trechos de rio")
2. MDET de baixa resolução
Para gerar o MDET de baixa resolução utilizaremos:
- MDET de alta resolução (
mdet_alta_fill) - a altitude máxima do trecho de rio a montante dada pelo raster
../softs/Paz2006_routines/Trechos5b/basin/cotaMon.rst - a altitude mínima do trecho de rio a jusante dada pelo raster
../softs/Paz2006_routines/Trechos5b/basin/cotaJus.rst
Lembre-se que é necessário gerar o arquivo .rdc com o Idrisi, a partir da cópia de um arquivo de baixa resolução que pode ser necessário ter seu formato alterado para real no Idrisi antes da atualização das estatísticas e da resolução do raster. A mediana desses 3 rasters será usada como aproximação do MDET da baixa resolução.
Vamos extrair as altitudes do mdet_alta_fill para os pixels exutórios das células e gerar um raster de baixa resolução com esses valores.
## MDET de alta resolução preenchido e recortado para bacia
mdet_alta_fill <- raster(mdet_alta_fill_file)
mdet_alta_fillclass : 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/MNTAlta/pauliceia/MNTfill.rst
names : MNTfill
values : 572, 751 (min, max)## raster com pixels exutórios (alta resolução) de cada célula
pixelexu_file <- paste0(trechos5b_basin_dir,"/pixelexu.rst")
pixelexu <- raster(pixelexu_file)
## substituindo pixels com valor zero por NA
pixelexu[pixelexu==0] <- 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. : NA
data source : in memory
names : pixelexu
values : 1, 1 (min, max)## extraindo altitude para os pontos exutórios de cada célula
elev_pixelexu <- mdet_alta_fill * pixelexu
plot(elev_pixelexu,
col = tim.colors(32),
main = "Elevação dos pixels exutórios das células")
## fator de salto entre as resoluções alta e baixa
fator <- unique(res(trechos_baixa))/unique(res(mdet_alta_fill))
fator[1] 5## gerando um raster dessas elevações na resolução do mdet de baixa resolução (1000 m)
## agregando píxels em células
elev_pixelexu_baixa <- aggregate(elev_pixelexu,
fact = fator,
FUN = "mean",
na.rm = TRUE)
## gráfico da elevação dos pixels exutórios
plot(elev_pixelexu_baixa,
col = tim.colors(32),
main = "raster gerado da elevação dos pixels exutórios")
Vamos importar os raster das altitudes dos pixels dentro do trecho de rio com maior altitude à montante e menor altitude à jusante (só é possível importar esses rasters após atualizar arquivos .rdc com Idrisi).
## raster de cota montante
cota_mon_file <- paste0(trechos5b_basin_dir, "/cotaMon.rst")
cota_mon <- raster(cota_mon_file)
## jusante
cota_jus_file <- paste0(trechos5b_basin_dir, "/cotaJus.rst")
cota_jus <- raster(cota_jus_file)
## diferença entre as cotas mon e jus deve ser positiva
plot(cota_mon-cota_jus)
A diferença entre as altitudes dos rasters cota_mon e cota_jus deve resultar em valores somente positivos. Após essa verificação, vamos calcular a mediana da altitude dos três rasters.
## empilhando os três rasters em um stack
s <- stack(cota_jus, elev_pixelexu_baixa, cota_mon)
#s <- stack(cota_jus, cota_mon)
sclass : RasterStack
dimensions : 14, 12, 168, 3 (nrow, ncol, ncell, nlayers)
resolution : 500, 500 (x, y)
extent : -2901.52, 3098.48, -3322.855, 3677.145 (xmin, xmax, ymin, ymax)
coord. ref. : NA
names : cotaJus, layer, cotaMon
min values : 572, 572, 576
max values : 734, 745, 745 ## calculando a mediana dos 3 rasters
wsdem <- calc(s, median)
## aplicando máscara
wsdem <- mask(wsdem, fracNA)
wsdemclass : 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 : layer
values : 589, 736 (min, max)## projetando wsdem
projection(wsdem) <- laea
plot(wsdem, main = "MDET da bacia em baixa resolução")
Salvando em ascii o MDET de baixa resolução só para bacia.
## salvando como ascii
wsdem_file_ascii <- gsub("basin", basin_name, "../wsdem_basin.asc")
raster2ascii_dbhm(r = wsdem, arqname = wsdem_file_ascii)Com MDET de baixa resolução para a bacia (wsdem) podemos calcular a declividade usando o método tradicional, usado por exemplo no software ArcGis, dada pela máxima diferença entre as cotas das 8 células vizinhas.
## calculando decividade pelo método do ArcGis
slope_baixa <- terrain(x = wsdem,
opt = "slope",
unit = "radians",
neighbors = 8)
rise_run_baixa <- tan(mask(slope_baixa, fracNA))
rise_run_baixa[rise_run_baixa < slp_lim_inf] <- slp_lim_inf
slope_baixa_pctg <- rise_run_baixa*100
plot(slope_baixa_pctg, main = "Declividade (%) usando o método tradicional")
## plot da diferença entre as declividades dos dois métodos
plot(slope_baixa_pctg - slope_iph_pctg,
main = "Diferença entre declividades: tradicional - upscaling")
## histograma da diferença entre as declividades dos dois métodos
update(histogram(slope_baixa_pctg - slope_iph_pctg, main = ""),
scales = list(cex = 1.2),
xlim = c(-20,20))
## salvando ascii
slope_baixa_pctg_file <- gsub("basin", basin_name, "../slope_basin.asc")
raster2ascii_dbhm(r = round(slope_baixa_pctg, 6),
arqname = slope_baixa_pctg_file)
## gerando slope para o sib2 baseado no slope_baixa_pctg
slpsib_baixa <- round(sin(atan(slope_baixa_pctg/100)), 6)
plot(slpsib_baixa, main = "Declividade SiB2 (radianos) \n usando o método tradicional")
## salvando slope sib2
slpsib_file <- gsub("basin", basin_name, "../slpsib_basin.asc")
raster2ascii_dbhm(r = round(slpsib_baixa, 6), arqname = slpsib_file)3. Distância a foz da bacia
Uma informação necessária para o DBHM é a distância de cada célula ao exutório da bacia. Esse raster pode ser obtido com a função flowLength disponível no script flowLength.R. Essa função tem como um de seus argumento o parâmetro mc_cores que deve ser definido de acordo com a configuração de sua máquina. No exemplo a máquina utilizada possui 8 núcleos, o mc_core = 7 significa que o looping em cada célula (dentro da função flowLength) é dividido, ou seja, paralelizado, entre os 7 núcleos. Altere o valor da variável ncores de acordo com a disponibilidade em sua máquina. Esta operação pode levar alguns minuto dependendo da máquina.
NOTA: altere o parâmetro
ncoresde acordo com a disponibilidade em sua máquina
## o número de núcleos de acordo com a disponibilidade de sua máquina (htop)
ncores <- 7
## gerando direção de fluxo baixa para bacia apenas
dir_baixa_mask <- mask(dir_baixa, fracNA)
## Garbage Collection (limpeza de memória)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 2435430 130.1 3886542 207.6 3886542 207.6
Vcells 2784687 21.3 4701432 35.9 3851194 29.4## distância ao exutório usando distância fixas
## (lado = 0.96194, diagonal = 1.36039, usando operadores de distancia, como comentado por De Smith (2004))
## AVISO: pode levar alguns minutos, ALTERE mc_cores de acordo com a disponibilidade de sua maquina (htop)Se o arquivo de distância (seguindo a rede de drenagem) ao exutório ../dis_basin.asc já tiver sido gerado previamente ele será lido, do contrário, ele então será calculado e exportado para ascii. Essa condição é para evitar refazer esse cálculo, pois demanda uma quantidade de tempo significativa dependendo do tamanho da bacia. Se for de interesse gerar novamente o ../dis_basin.asc basta apagá-lo.
dis_file <- gsub("basin", basin_name,"../dis_basin.asc")
# if(file.exists(dis_file)) {
# dis <- raster(dis_file)
# } else {
dis <- flowLength(fd = dir_baixa_mask, mc_cores = ncores)elapsed
0.043 ## salvando ascii
raster2ascii_dbhm(r = dis, arqname = dis_file)
# } # end if dis_fileO mesmo procedimento será repetido para a distância ao exutório usando os comprimentos de trechos de rio obtidos dos rasters de alta resolução.
## nome do arquivo dis_iph em formato ascii
dis_iph_file <- gsub("basin", basin_name, "../dis_iph_basin.asc")
# if(file.exists(dis_iph_file)) {
# dis_iph <- raster(dis_iph_file)
# } else {
## usando comprimento de trechos de alta resolução, saída da Trechos5b
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 2439437 130.3 3886542 207.6 3886542 207.6
Vcells 2791646 21.3 4701432 35.9 3851194 29.4 dis_iph <- flowLength(fd = dir_baixa_mask,
fl = reaches_mask ,
mc_cores = ncores)elapsed
0.044 ## salvando ascii
raster2ascii_dbhm(r = dis, arqname = dis_file)
# } # end if dis_fileVamos fazer um gráfico de comparação entre as distâncias à Foz geradas pelas diferentes formas.
Qual o ponto mais distante do exutório? Vamos salvar o caminho do ponto mais distante do exutório em arquivo shapefile. Para isso precisamos dos rasters de direção de fluxo e de área de drenagem acumulada no formato de entrada do DBHM (dir_basin.asc, acc_basin.asc) gerados na etapa anterior.
## area de drenagem acum.
acc_file <- gsub("basin", basin_name, "../acc_basin.asc")
acc <- raster(acc_file)
acc <- setMinMax(acc)
|
| | 0%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 : /home/hidro2roilan/DBHM/data_process/geo_data/acc_pauliceia.asc
names : acc_pauliceia
values : 0, 44 (min, max)## direção de fluxo na convenção requerida pelo DBHM (ArcGis)
dir_file <- gsub("basin", basin_name, "../dir_basin.asc")
fdir <- raster(dir_file)
fdir <- setMinMax(fdir)
|
| | 0%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. : NA
data source : /home/hidro2roilan/DBHM/data_process/geo_data/dir_pauliceia.asc
names : dir_pauliceia
values : 1, 128 (min, max)Vamos determinar a(s) célula(s) mais distante(s) do exutório da bacia.
## distância máxima para cada dis
mx <- cellStats(dis, stat = max)
mx[1] 5407.375mx_iph <- cellStats(dis_iph, stat = max)
mx_iph[1] 5067.67## célula mais distante
cell_dis_max <- Which(dis == mx, cells = T)
cell_dis_max[1] 46## se houver mais de uma célula, escolhemos arbitrariamente a primeira
cell_dis_max <- cell_dis_max[1]E agora podemos determinar a trajetória ao longo da rede de drenagem da célula mais distante do exutório até o exutório. Novamente, devido ao tempo de execução da função flowPath, se o arquivo shapefile com a trajetória já existir de uma rodada prévia ele será usado e o cálculo da trajetória não será necessário. Do contrário o cálculo é refeito. Se for de interesse gerar novamente o arquivo shapefile ../longestPath_pauliceia.shp e os arquivos associados ao shapefile, basta apagá-los.
## para não mostrar barra de progresso
rasterOptions(progress = "")
## nome do arquivo shapefile com trajetória
splines_longest_file <- paste0(dirname(dis_iph_file), "/longestPath_", basin_name,".shp")
## se arquivo existir ele será lido
# if(file.exists(splines_longest_file)){
# ## le arquivo shapefile
# splines_longest <- shapefile(splines_longest_file)
# ## células que intersectam o caminho no wsdem
# fp_cells <- extract(x = wsdem, y = splines_longest, cellnumbers = TRUE )
# fp_cells <- fp_cells[[1]][,"cell"]
# } else {
## caminho a jusante partindo da célula mais longe
fp_cells <- flowPath(x = fdir, p = cell_dis_max)
## coordenadas das fp_cells
coords_fp_cells <- data.frame(xyFromCell(dis, fp_cells),
acc = raster::extract(acc, y = fp_cells))
head(coords_fp_cells) x y acc
1 1848.48 1927.14519 0
2 1348.48 1427.14519 2
3 1348.48 927.14519 4
4 848.48 427.14519 6
5 848.48 -72.85481 15
6 348.48 -72.85481 16 ## coordenadas em ordem da área acumulada
coords_fp_cells_o <- coords_fp_cells[order(coords_fp_cells$acc),]
## criando spatialLines a partir das coordenadas
splines_longest <- SpatialLines(list(Lines(list(Line(coords_fp_cells_o[,1:2])), ID = "1")))
proj4string(splines_longest) <- laea
## salvando como shapefile
shapefile(x = splines_longest, splines_longest_file, overwrite = TRUE)
# } # end if splines_longest_file
## reset barprogress
rasterOptions(progress = "text")Gráfico da trajetória até o ponto da bacia mais distante do exutório.
plot(wsdem, main = "Trajetória de escoamento do ponto\n mais a montante no MDET")
plot(splines_longest, add = TRUE, col = "blue", lwd = 2)
plot(shapefile(x = "../softs/TauDEM-master/pauliceia/exutorio_baixa.shp"), add = T, lty = 6, lwd = 3, col = 2)
Gráfico da elevação ao longo da distância à jusante do ponto mais a montante na bacia.
## Perfil da elevação ao longo do caminho entre o exutório e o ponto mais a montante.
## o "::" é uma forma de acessar a função dentro de um pacote
## útil para situações em que uma função tem mesmo nome em diferentes pacotes
elev_fp <- raster::extract(x = wsdem, y = fp_cells)
## distâncias
dis_fp <- raster::extract(x = dis, y = fp_cells)
dis_iph_fp <- raster::extract(x = dis_iph, y = fp_cells)
x_range <- range(c(dis_iph_fp[], dis_iph[])/1000, na.rm = T)
plot(x = dis_iph_fp/1000,
y = elev_fp,
type = "o",
pch = 20,
xlim = x_range,
xlab = "distância à montante do exutório (Km)",
ylab = "elevação (m)",
las = 1,
main = "Perfil de elevação ao longo da trajetória \n do ponto mais a montante do MDET")
lines(dis_fp/1000, elev_fp, type = "o", pch = 20, col = 2)
eaxis(1, labels = F)
eaxis(2, labels = F)
## legenda
legend("topleft",
legend = c("tradicional", "IPH"),
lty = 1,
pch = 20,
col = c(2,1),
bty = "n")
4. Área da células do MDET
No DBHM é necessário um arquivo raster denominado gridarea.asc com a área de cada célula, que pode ser gerado da seguinte forma.
## área de cada célula da grade
akm2 <- sqrt(prod(res(dir_baixa)))/1000
## estrutura do raster de baixa resolução
gridarea <- raster(dir_baixa)
## atribuindo valores de 1 km2 para cada céula do raster
gridarea <- setValues(gridarea, akm2)
## raster com valores só na área da bacia
gridarea <- mask(gridarea, fracNA)
plot(gridarea, main = "área de cada célula: .90 km²")
## nome do arquivo gridarea
gridarea_file <- gsub("basin", basin_name, "../gridarea_basin.asc")
## salvando como ascii
raster2ascii_dbhm(r = gridarea, arqname = gridarea_file)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 verificação de arquivos
- wsdem.asc MDET da bacia com -9999 (NA) atribuído para células fora da área da bacia;
wsdemclass : 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 : 589, 736 (min, max)wsdem_file_ascii[1] "../wsdem_pauliceia.asc"- gridarea.asc área de cada célula (km²);
gridareaclass : 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.5, 0.5 (min, max)gridarea_file[1] "../gridarea_pauliceia.asc"- slope.asc declividade do terreno (%);
slope_iph_pctgclass : 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 : layer
values : 0.01, 10.09227 (min, max)slope_iph_pctg_file[1] "../slope_iph_pauliceia.asc"slope_iph_pctgclass : 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 : layer
values : 0.01, 10.09227 (min, max)slope_iph_pctg_file[1] "../slope_iph_pauliceia.asc"- slpsib.asc declividade do terreno (radianos) para o SiB2;
slpsib_iphclass : 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 : layer
values : 1e-04, 0.100413 (min, max)slpsib_iph_file[1] "../slpsib_iph_pauliceia.asc"slpsib_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 : layer
values : 0.016505, 0.073833 (min, max)slpsib_file[1] "../slpsib_pauliceia.asc"- dis.asc distância a foz da bacia (km);
disclass : 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 : 480.97, 5407.375 (min, max)dis_file[1] "../dis_pauliceia.asc"dis_iphclass : 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 : 300, 5067.67 (min, max)dis_iph_file[1] "../dis_iph_pauliceia.asc"- rivlen.asc comprimento (m) dos trechos de rio que conectam células vizinhas;
reaches_maskclass : 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 : layer
values : 300, 1298.76 (min, max)rivlen_sib_file[1] "../rivlen_pauliceia.asc"7. Referências
De Smith, M. (2004). Distance transforms as a new tool in spatial analysis, urban planning, and GIS. Environment and Planning B: Planning and Design 31, pp. 85-104.
Fairfield, J., e P. Leymarie (1991), Drainage networks from grid digital elevation models, Water Resourc. Res., 27(5), 709-717, doi:0043-1397/91/90WR-02658.
Jenson, S. K., e Domingue, J.O. (1988), Extracting topographic structure from digital elevation data for geographic information system analysis, Photogramm. Eng. Remote Sens., 54(11), 1593-1600.
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.
Paz, A. R. ; Collischonn, W. River reach length and slope estimates for large-scale hydrological models based on a relatively high-resolution digital elevation model. Journal of Hydrology (Amsterdam), v. 343, p. 127-139, 2007.
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.