DBHM-GIS: PARTE 5
17-11-2014
Universidade Federal de Santa Maria - UFSM
Departamento de Física
Prof. Jônatan Tatsch
Aluno Roilan Hernandez Valdes
1. Discretização da bacia hidrográfica
Para a gestão de recursos hídricos é frequentemente necessário entender como mudanças (uso da terra, obras de engenharia) em uma parte da bacia devem afetar o uso da água em outra parte da bacia. O conhecimentos desses efeitos implica na relação topológica entre diferentes partes da bacia, isto é, saber qual parte (sub-bacia) recebe (drena) escoamento de (para) uma dada parte (sub-bacia). A relação topológica dentro da bacia é muito importante para simulação hidrológica para propagar o runoff ao longo da rede drenagem até o exutório da bacia. Esse conhecimento pode ser obtido das informações derivadas do MDET: direção de fluxo, área de drenagem acumulada e distância ao exutório. O esquema de numeração de Pfafstetter é uma metodologia amplamente reconhecida para descrição da topologiada bacia hidrográfica e é usado no DBHM. Ela baseia-se no controle topográfico e na topologia da rede de drenagem. O sistema é fundamentado nos conceitos articulados pioneiramente por Pfafstetter (1989) e documentado por Verdin e Verdin (1999). O esquema de numeração usa um arranjo hierárquico de dígitos decimais e as bacias são delineadas a partir das junções com da rede de drenagem. A hierarquia implica em diferentes níveis de subdivisão das sub-bacias.
Para cada nível, cada sub-bacia é definida por um número inteiro m variando de m = 0 a 9, baseado na sua localização e sua função dentro da rede de drenagem. Os valores de 0 a 9 são tratados de forma binária como ímpares ou pares. O valor ordinário de um dígito indica a posição relativa: montante ou jusante, enquanto que a paridade (ímpar/par) indica a posição na rede de drenagem: dentro ou fora do rio principal.
O nível 0 corresponderia a bacia continental ou alternativamente a bacia que drena para o oceano. Níveis superiores representa progressivamente subdivisões mais detalhadas do nível 0 da bacia.
Para as bacias de nível 1, a procura por sub-bacias começa no exutório, buscando-se junções entre sub-bacias por meio dos rasters de direção de fluxo e área acumulada. Primeiro, detecta-se a célula a montante nas células vizinhas de acordo com a direção de fluxo. Então verifica-se o valor da área acumulada da célula à montante. A célula à montante com maior área é reconhecida como parte do rio principal e as outras células à montante são assumidas como células tributárias. Seguindo ao longo do caminho das direções de fluxo até atingir a borda da bacia, as células do rio principal com as 4 maiores áreas são registradas. Essas células são as junções de sub-bacias e são a partição (nó) entre sub-bacias, codificadas como sub-bacias 2, 4, 6, e 8 de jusante para montante. As bacias entre essas 4 sub-bacias (inter-bacias) são codificadas como 1, 3, 5, 7, e 9 de jusante para montante. A codificada como m = 9 é reservada para bacia de cabeceira.
O procedimento de sub-divisão é realizado iterativamente para obter um conjunto de sub-bacias (ou ottobacias) de nível 2 sobre as bacias de nível 1. Dentro de cada bacia de nível 1, ou de cada interbacia, rio principal e tributários são identificados usando o mesmo método. Os tributários são então ordenados pelo valor de área acumulada e os quatro maiores são númerados 2, 4, 6, e 8, indo de jusante para montante ao longo do rio principal. A Figura 1 ilustra sub-bacias de código Pfafstetter nível 2.
As áreas de drenagem são derivadas usando a direção de fluxo e a todas células são acrescentadas o segundo dígito de Pfafstetter. Teoricamente, o sistema não é limitado ao em um determinado número n de níveis. Na prática esse número de níveis é limitado pela resolução do MDET, n = 6 a 8 níveis é geralmente suficiente. Isso torna-se evidente quando não é mais possível identificar 4 tributários e o número da bacia naquele nível é codificado como 0 no sistema do DBHM.
A discretização da bacia em sub-bacias será realizada com o código fortran priver_tatsch.f com pequenas modificações em relação a versão original de Tang (2006). O código está disponível na página da disciplina no Moodle. Baixe o arquivo e salve no diretório softs. A seguir criamos dentro desse diretório um sub-diretório nomeado de acordo com a bacia de interesse.
2. Preparação
Para executar um código fortran a partir do R, precisamos criar um link simbólico para o ifort. O que poder ser feito por meio do seguinte comando no terminal linux:
$ sudo ln -s /opt/intel/fc/10.1.018/bin/ifort /usr/bin/ifortA seguir, vamos configurar os diretórios, o usuário e carregar os pacotes e scripts necessários.
## carregando pacotes necessários
pcks <- c("raster", "rgdal", "RSAGA", "adehabitat", "maptools", "sp", "fields", "sfsmisc", "doBy", "stringr", "lattice", "rasterVis")
## 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 ## versão de cada pacote
pcks_version <- unlist(lapply(pcks, function(x) as.character(packageVersion(x))))
## combinando informações em um dataframe
pcks_info <- data.frame(pcks, pcks_version)
pcks_info pcks pcks_version
1 raster 2.3.12
2 rgdal 0.9.2
3 RSAGA 0.93.6
4 adehabitat 1.8.15
5 maptools 0.8.34
6 sp 1.0.17
7 fields 8.2.1
8 sfsmisc 1.0.27
9 doBy 4.5.13
10 stringr 0.6.2
11 lattice 0.20.30
12 rasterVis 0.35## configurando opção para mostrar barra de progresso durante operações com raster
rasterOptions(progress = "text")NOTA: altere o parâmetro
usernamede acordo com a sua configuração
## defina seu nome de usuario que definira o caminho até o diretório scripts
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
my_path <- gsub(pattern = "user",
replacement = user_name,
x = "/home/user/DBHM/")
my_path[1] "/home/hidro2roilan/DBHM/"## definindo diretório de trabalho (combina my_path e caminho até subdir scripts)
script_dir <- paste0(my_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"## função para converter raster em polígono
source("gdal_polygonizeR.R")
## converção para ascii no formato dbhm
source("raster2ascii_dbhm.R")Vamos criar o diretório para execução da rotina de discretização da bacia (../softs/priver) onde deve ser salvo o arquivo priver_tatsch.f, fornecido através da página do curso no moodle. Dentro deste diretório criaremos o diretório da bacia ../softs/priver/basin. A string ‘basin’ será substituída pelo nome da bacia, definido pela variável basin_name, semelhante ao realizado nas etapas anteriores.
## define o nome da bacia
##basin_name <- "mogu"
basin_name <- "pauliceia"
#basin_name <- "vacacai"
## diretório para a rotina priver
priver_dir <- "../softs/priver"
## se não existir o diretório priver_dir, ele será criado
if(!basename(priver_dir) %in% dir("../softs")) dir.create(priver_dir)
## define nome do sub-diretório da bacia
basin_dir <- gsub(pattern = "basin",
replacement = basin_name,
x = paste0(priver_dir,"/basin")
)
## se ele não existir, ele será criado
if(!basename(basin_dir) %in% dir(priver_dir)) dir.create(basin_dir)Alguns dos arquivos raster gerados nas etapas anteriores (acc.asc, dir.asc, dis.asc e wsdem.asc) são as entradas da rotina priver_tatsch.f. Vamos copiá-los para o diretório ../softs/priver/pauliceia.
## lista de arquivos (files list)
fl <- list.files("..", pattern = "\\.asc$", full.names = T)
fl <- fl[grep(basin_name, basename(fl))]
fl <- fl[grep("wsdem|dis|dir|acc", basename(fl), value = F)]
fl[1] "../acc_pauliceia.asc" "../dir_pauliceia.asc"
[3] "../dis_pauliceia.asc" "../wsdem_pauliceia.asc"## copiar arquivos raster (*.asc) para diretório priver/basin
dest_dir <- paste0(basin_dir, "/",
gsub(pattern = paste0("_", basin_name),
replacement = "",
x = basename(fl)
)# end gsub
)# end paste0
## copiando arquivos ascii
file.copy(from = fl,
to = dest_dir,
overwrite = T)[1] TRUE TRUE TRUE TRUEVamos carregar o raster de área de drenagem acumulada (acc.asc) para obtermos a estrutura espacial do domínio da bacia, como: extent, resolução, número de linhas, colunas e etc.
## caminho completo até o arquivo asc da bacia (WaterShed DEM)
acc_file <- paste0(basin_dir,"/acc.asc")
file.exists(acc_file)[1] TRUE## lendo raster do modelo digital de elevação do terreno reamostrado de 90 m para 1km
acc <- raster(acc_file)
## informações do raster acc
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/softs/priver/pauliceia/acc.asc
names : acc
values : 0, 44 (min, max)plot(acc)
Para definir a projeção do raster acc será definida a partir da string com os parâmetros da projeção cartográfica Azimutal de área equivalente de Lambert (LAEA) obtida da leitura do arquivo shapefile do polígono da bacia (contorno) determinado da área de drenagem de alta resolução.
## arquivos shapefile no diretório
pol_basin_file <- list.files(path = gsub(pattern = "basin",
replacement = basin_name,
x = "../softs/TauDEM-master/basin"),
pattern = "\\.shp$",
full.names = T)
## seleção do arquivo de interesse
pol_basin_file <- grep("mascara_bacia_alta_pol",pol_basin_file, value = T)
## leitura do shapefile
pol_basin_alta <- shapefile(pol_basin_file)
## extraindo string com parâmetros da projeção
laea <- proj4string(pol_basin_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"## definindo projeção do acc
projection(acc) <- laea
## plot
plot(acc, col = tim.colors(32))
plot(pol_basin_alta, add = T, border = "green", lwd = 2)
Vamos obter algumas informações que definem o domínio espacial do raster acc e armazená-las em um vetor do tipo caractere. A função sprintf imprime o valor numérico com número de dígitos desejado após a vírgula (indicado pelo argumento, fmt = %.7f). Esse valor deve ser alterado conforme a precisão usada no cabeçalho do arquivo acc.asc.
## número de linha e colunas do raster acc
new_row <- as.character(nrow(acc))
new_row[1] "14"new_col <- as.character(ncol(acc))
new_col[1] "12"## Valores das coordenadas do canto inferior esquerdo do domínio espacial
## xmin
xmin_acc <- sprintf(fmt = "%.7f", xmin(acc))
xmin_acc <- paste0("'",xmin_acc,"'")
xmin_acc[1] "'-2901.5200449'"## ymin
ymin_acc <- sprintf("%.7f", ymin(acc))
ymin_acc <- paste0("'",ymin_acc,"'")
ymin_acc[1] "'-3322.8548110'"## resolucao
cell_size <- mean(res(acc))
## vetor do tipo character combinando essas 4 informações
new_strings <- as.character(c(new_col, new_row, xmin_acc, ymin_acc, cell_size))
new_strings[1] "12" "14" "'-2901.5200449'" "'-3322.8548110'"
[5] "500" 3. Aplicação da rotina de discretização da bacia
O código fortran priver_tatsch.f faz a discretização da bacia usando o sistema de Pfafstetter. Nessa rotina as dimensões das variáveis e os loopings nas linhas e colunas são definidos pelo domínio e a resolução espacial do raster da bacia. Por isso é necessário alterar esses valores conforme o domínio da bacia a ser discretizada. As linhas e colunas das matrizes da bacia no código fortran estão definidas pelas strings ‘yyy’ e ‘xxx’, respectivamente. Vamos criar um vetor (do tipo character no R) com essas strings:
old_strings <- c("xxx","yyy", "XLLCORNER", "YLLCORNER", "sizecell")
#old_strings <- c("xxx","yyy")
old_strings[1] "xxx" "yyy" "XLLCORNER" "YLLCORNER" "sizecell" Antes de executar a rotina precisamos substituir essas strings pelos valores correspondentes ao da bacia em questão. Vamos importar o código fortran.
## nome arquivo fortran
fortran_file <- paste0(priver_dir, "/priver_tatsch.f")
fortran_file[1] "../softs/priver/priver_tatsch.f"file.exists(fortran_file)[1] TRUE## leitura do código fortran em linhas
fcode <- readLines(con = fortran_file)Warning in readLines(con = fortran_file): incomplete final line found on
'../softs/priver/priver_tatsch.f'# encontre as dimensões de acc(yyy, xxx)
head(fcode, 20) [1] "c -------------------------------------------------------------------"
[2] "c\tThe Program can code the drainage basin by Pfafstetter coding system (3 levels)."
[3] "c\tInput: acctxt, dirtxt,(dem,dis)"
[4] "c\tOutput: basin0, basin1, basin2, basin3"
[5] "c\toutput 2: mainstem grid, distance-altitue graph"
[6] "c\tContact with: (TANG, Qiuhong) tangqh@iis.u-tokyo.ac.jp"
[7] "c -------------------------------------------------------------------"
[8] "c Atualizacao do codigo para leitura de ascii grid salvo no R"
[9] "c Jonatan Tatsch, 07 jun 2010"
[10] "c Mudanca na leitura do cabecalho"
[11] "c -------------------------------------------------------------------"
[12] ""
[13] " character*5 ncols,nrows"
[14] " character*14 xllcorner,yllcorner"
[15] " character*14 cellsize"
[16] " character*14 nodata"
[17] " integer x0,y0"
[18] " integer s, znodata,nc,nr"
[19] " integer acc(yyy,xxx), dir(yyy,xxx)"
[20] " integer maxtmp, outi, outj, max1i,max1j,surround(8,3),max1" Vamos identificar as linhas do código com as strings xxx, yyy, XLLCORNER, YLLCORNER, sizecell e depois substituí-las pelas strings correspondentes às da bacia: 12, 14, ‘-2901.5200449’, ‘-3322.8548110’, 500.
## definindo regex para as strings procuradas
target_string <- paste(old_strings, collapse="|")
target_string[1] "xxx|yyy|XLLCORNER|YLLCORNER|sizecell" ## procura todas linhas (posições) com row_string ou col_string
pos <- grep(target_string, fcode)
## linhas do código com o padrão procurado
fcode[pos] [1] " integer acc(yyy,xxx), dir(yyy,xxx)"
[2] " integer outlet(2000,2),id,basin(yyy,xxx)"
[3] " integer basin2(yyy,xxx),basin3(yyy,xxx)"
[4] "\treal dis(yyy,xxx),dem(yyy,xxx)"
[5] " write(1,'(A14, A15)') xllcorner,XLLCORNER"
[6] " write(1,'(A14, A15)') yllcorner,YLLCORNER"
[7] " write(1,'(A14, A15)') xllcorner,XLLCORNER"
[8] " write(1,'(A14, A15)') yllcorner,YLLCORNER"
[9] " write(1,'(A14, A15)') xllcorner,XLLCORNER "
[10] " write(1,'(A14, A15)') yllcorner,YLLCORNER "
[11] " write(2,'(A14, A15)') xllcorner,XLLCORNER"
[12] " write(2,'(A14, A15)') yllcorner,YLLCORNER"
[13] "\t integer basin2(yyy,xxx)"
[14] " integer acc(yyy,xxx), dir(yyy,xxx)"
[15] " integer basin(yyy,xxx)"
[16] " integer acc(yyy,xxx), dir(yyy,xxx) "
[17] " integer acc(yyy,xxx), dir(yyy,xxx),mainstem(yyy,xxx) "
[18] "\tinteger dem(yyy,xxx),dis(yyy,xxx),ooi,ooj"
[19] "\tdo i=1,yyy"
[20] "\t\tdo j=1,xxx"
[21] " write(1,'(A, i)') 'ncols',xxx"
[22] " write(1,'(A, i)') 'nrows',yyy"
[23] " write(1,'(A, A20)') 'xllcorner',XLLCORNER"
[24] " write(1,'(A, A20)') 'yllcorner',YLLCORNER"
[25] " write(1,'(A, i)') 'cellsize',sizecell"
[26] "\tdo i=1,yyy"
[27] " \tdo j=1,xxx"
[28] " \tinteger acc(yyy,xxx), dir(yyy,xxx)"
[29] " \tinteger acc(yyy,xxx), dir(yyy,xxx), basin(yyy,xxx)" ## cópia do código original
new_code <- fcode
## lopping em cada old_string para substituir por cada new_string
for(i in 1:length(old_strings)) new_code <- gsub(old_strings[i], new_strings[i], new_code)
new_code[pos] [1] " integer acc(14,12), dir(14,12)"
[2] " integer outlet(2000,2),id,basin(14,12)"
[3] " integer basin2(14,12),basin3(14,12)"
[4] "\treal dis(14,12),dem(14,12)"
[5] " write(1,'(A14, A15)') xllcorner,'-2901.5200449'"
[6] " write(1,'(A14, A15)') yllcorner,'-3322.8548110'"
[7] " write(1,'(A14, A15)') xllcorner,'-2901.5200449'"
[8] " write(1,'(A14, A15)') yllcorner,'-3322.8548110'"
[9] " write(1,'(A14, A15)') xllcorner,'-2901.5200449' "
[10] " write(1,'(A14, A15)') yllcorner,'-3322.8548110' "
[11] " write(2,'(A14, A15)') xllcorner,'-2901.5200449'"
[12] " write(2,'(A14, A15)') yllcorner,'-3322.8548110'"
[13] "\t integer basin2(14,12)"
[14] " integer acc(14,12), dir(14,12)"
[15] " integer basin(14,12)"
[16] " integer acc(14,12), dir(14,12) "
[17] " integer acc(14,12), dir(14,12),mainstem(14,12) "
[18] "\tinteger dem(14,12),dis(14,12),ooi,ooj"
[19] "\tdo i=1,14"
[20] "\t\tdo j=1,12"
[21] " write(1,'(A, i)') 'ncols',12"
[22] " write(1,'(A, i)') 'nrows',14"
[23] " write(1,'(A, A20)') 'xllcorner','-2901.5200449'"
[24] " write(1,'(A, A20)') 'yllcorner','-3322.8548110'"
[25] " write(1,'(A, i)') 'cellsize',500"
[26] "\tdo i=1,14"
[27] " \tdo j=1,12"
[28] " \tinteger acc(14,12), dir(14,12)"
[29] " \tinteger acc(14,12), dir(14,12), basin(14,12)" Agora vamos salvar o código fortran modificado no diretório da bacia acrescentando o sufixo com nome da bacia.
## escrevendo novo código fortran
## novo nome para o arquivo
new_fortran_file_name <- paste0(gsub(pattern = "\\.f",
replacement = paste0("_", basin_name),
x = basename(fortran_file)
),
".f")
new_fortran_file_name[1] "priver_tatsch_pauliceia.f"new_fortran_file <- paste0(basin_dir, "/",new_fortran_file_name)
new_fortran_file[1] "../softs/priver/pauliceia/priver_tatsch_pauliceia.f"## escreve código modificado para arquivo
writeLines(text = new_code, con = new_fortran_file)Vamos entrar no diretório da bacia ../softs/priver/pauliceia para compilar e executar o código modificado priver_tatsch_pauliceia.f. Após isso, redefinimos o diretório de trabalho como o diretório de scripts /home/hidro2roilan/DBHM/data_process/geo_data/scripts.
setwd(basin_dir)
dir() [1] "acc.asc" "basin0.asc"
[3] "basin1.asc" "basin3_pol.dbf"
[5] "basin3_pol.prj" "basin3_pol.shp"
[7] "basin3_pol.shx" "dir.asc"
[9] "dis_alt.txt" "dis.asc"
[11] "mainstem.asc" "mainstem.dbf"
[13] "mainstem.prj" "mainstem.shp"
[15] "mainstem.shx" "pp_asciigrid_files"
[17] "priver_tatsch_pauliceia" "priver_tatsch_pauliceia.f"
[19] "wsdem.asc" ## nome do executável
exec <- gsub("\\.f", "", basename(new_fortran_file))
exec[1] "priver_tatsch_pauliceia" ## comando apra executar a rotina
cmd <- paste0("ifort ", basename(new_fortran_file), " -traceback -o ", exec, paste0("; ./", exec))
cmd[1] "ifort priver_tatsch_pauliceia.f -traceback -o priver_tatsch_pauliceia; ./priver_tatsch_pauliceia" ## compilando e executando
system(command = cmd, intern = T, wait = TRUE)[1] " wsdem lido"
[2] " nc,nr,x0,y0,size 12 14 -986359726 -984633939 500"
[3] " outlet acc = 44 outi= 10 outj= 5"
[4] " level 1 OK!"
[5] " 0 0 0 error!" ## check
# system("nedit *.asc &")
## retornando ao diretório de scripts
setwd(script_dir)4. Verificação dos rasters de saída
Vamos verificar as os arquivos ascii de saída da rotina de discretização da bacia em um objeto da classe stack do pacote raster, que permite empilhar diversos raster em uma única variável.
ascii_basin_files <- list.files(path = basin_dir,
pattern = "basin*.\\.asc",
full.names = T)
## stack com os 4 rasters de discretização da bacia
s <- stack(ascii_basin_files)
## aplicando mascara ao stack
s <- mask(s, acc)
plot(s)
Para melhor visualização dos níveis de discretização da bacia hidrográfica vamos utilizar cores aleatórias para cada sub-bacia.
## Gráfico para visuzlização das sub-bacias em diferentes níveis
## número total de cores disponíveis
ncores <- length(colors())
ncores[1] 657## parâmetros para tela gráfica
op <- par(mfrow = c(nlayers(s)/2, 2), mar = c(1, 1, 3, 4), las = 1)
## lopping em cada layer do stack composto pelos 4 arquivos ascii
for(i in 1:nlayers(s)) {
# i <- 1
## subset: seleciona um subconjunto do stack
ss <- subset(s, i)
## número total de sub-bacias
nbs <- length(sort(unique(getValues(ss))))
## sorteia aleatóriamente nbs cores das ncores
plot_cols <- colors()[sample(1:ncores, nbs)]
plot(ss, col = plot_cols, axes = F)
}# end for
## configura parâmetros para os valores default
par(op)Vamos contar o número de células por sub-bacia para os diferentes níveis de discretização.
## gráfico de barras
## definindo novos parâmetros para tela gráfica
op <-par(mfrow = c(nlayers(s)/2, 2), mar=c(4,4,2,0), oma=c(1.5,2,1,1))
## lopping em cada layer do stack composto pelos 4 arquivos ascii
l_grids_per_sub <- lapply(1:(nlayers(s)),
##
function(il){
# il = 3
## contagem de grades por
grids_per_b <- table(mask(subset(s, il), acc)[])
#grids_per_b[order(as.integer(names(grids_per_b)))]
barplot(sort(grids_per_b), las = 1)
eaxis(2,labels = F)
box()
mtext(text = paste0("Pontos de grade por subbacia",
"\n nível = ", il-1,
", máx = ", max(grids_per_b),
", min = ", min(grids_per_b)),
line = 0,
cex = 1)
grids_per_b
}) 
par(op)Entre os parâmetros de configuração do DBHM estão o número máximo de células dentro das sub-bacias e número máximo de sub-bacias. Vamos determinar esses valores e gerar um arquivo shapefile com os polígonos das sub-bacias de nível 3.
## número máximo de células por sub-bacia para os diferentes níveis
## parâmetro de entrada do DBHM
names(l_grids_per_sub) <- names(s)
max_grids_per_sub <- lapply(l_grids_per_sub, function(tab)tab[which.max(tab)])
max_grids_per_sub$basin0
1
45
$basin1
2
17 ## total de sub-bacias máximo nível de discrtização, em geral nível 3
## discretização de nível 3
b3 <- raster(s, layer = nlayers(s))
## intervalo de variação dos códigos de bacias nível 3
subb3_codes <- sort(unique(b3))
## mínimo e máximo código do nível 3
range(subb3_codes, na.rm = T)[1] 1 9 ## total de sub-bacias
max_sub <- length(subb3_codes)
max_sub[1] 9 ## código oi ID da bacia com máximo de células dentro
id_max_sub3 <- as.integer(names(max_grids_per_sub[[nlayers(s)]]))
id_max_sub3[1] 2 plot(b3 == id_max_sub3,
main = paste0("sub-bacia com maior número de células (", id_max_sub3, ")"))
## gerando polígonos de b3
b3_pol <- gdal_polygonizeR(x = b3, pypath = "/usr/bin/gdal_polygonize.py")
proj4string(b3_pol) <- laea
plot(b3_pol, add = T) 
## salvando como shapefile
b3_pol_file <- paste0(basin_dir, "/basin3_pol.shp")
shapefile(b3_pol, b3_pol_file, overwrite = TRUE)O DBHM pode ser usado com qualquer um dos níveis de discretização, dependendo da capacidade computacional disponível. No nosso caso optaremos pela discretização de nível 3, que para a bacia pauliceia:
o número máximo de sub-bacias é de 9;
o número de células dentro das sub-bacias é ;
5. Pós processamento dos rasters de saída
Os raster de saída da rotina priver (basin*.asc) serão salvos conforme o formato ascii do DBHM em um sub-diretório chamado ../softs/priver/pauliceia/pp_asciigrid_files.
## criando diretório para guardar arquivos ascii consistentes com formato dbhm
pp_dir <- paste0(basin_dir,"/pp_asciigrid_files")
dir.create(pp_dir)Warning in dir.create(pp_dir):
'../softs/priver/pauliceia/pp_asciigrid_files' already exists## por precaução usamos a extent do acc e a atribuímos à do stack s
extent(s) <- extent(acc)
## definindo valor para NA
NAvalue(s) <- -9999
## aplicando máscara ao stack
s <- mask(s, acc)
# writeRaster(x = s, filename = "t.asc", bylayer = TRUE, format = "ascii", suffix = names(s))
## looping para gerar ascii para cada layer do stack (s)
l_write_files <- lapply(1:nlayers(s),
function(i) {
# i <- 1
cat(i, "\n")
r <- raster(x = s, layer = i)
## output file name
ofn <- paste0(pp_dir, "/",names(s)[i], ".asc")
## escreve raster em formato ascii para o DBHM
raster2ascii_dbhm(r = r, arqname = ofn)
if(file.exists(ofn)) result <- ofn else result <- "arquivo não gerado"
return(result)
} # end fun
) ## end lapply1
2 l_write_files[[1]]
[1] "../softs/priver/pauliceia/pp_asciigrid_files/basin0.asc"
[[2]]
[1] "../softs/priver/pauliceia/pp_asciigrid_files/basin1.asc"6. Rio principal
Vamos importar o arquivo raster de saída mainstem.asc da rotina priver que localiza as células do rio principal da bacia.
mainstem_file <- paste0(basin_dir, "/mainstem.asc")
## raster mainstem com máscara aplicada
mainstem <- mask(raster(mainstem_file), acc)
mainstemclass : 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 : mainstem
values : 0, 1 (min, max)plot(mainstem, main = "Rio principal")
## salvando arquivo mainstem com máscara para bacia
raster2ascii_dbhm(r = mainstem, arqname = mainstem_file)As coordenadas do rio principal podem ser úteis futuramente para gráficos de caracterização da bacia. Então vamos salvá-las em um arquivo shapefile.
## células do rio principal
cells_mainstem <- Which(mainstem == 1, cells = T)
## coordenadas do rio principal
coords_mainstem <- data.frame(xyFromCell(mainstem, cells_mainstem),
acc = extract(acc, y = cells_mainstem))
head(coords_mainstem) x y acc
1 848.48 2427.1452 0
2 -151.52 1927.1452 5
3 348.48 1927.1452 2
4 -151.52 1427.1452 8
5 -151.52 927.1452 15
6 -151.52 427.1452 16 ## coordenadas em ordem da área acumulada
coords_mainstem_o <- coords_mainstem[order(coords_mainstem$acc),]
## criando spatialLines a partir das coordenadas
splines_mainstem <- SpatialLines(list(Lines(list(Line(coords_mainstem_o[,1:2])), ID = "1")))
proj4string(splines_mainstem) <- laea
## plot
plot(b3, col = plot_cols)
plot(b3_pol, add = T)
plot(splines_mainstem, add = T, lwd =2, col = "blue")
## salvando como shapefile
splines_mainstem_file <- gsub("\\.asc", "\\.shp", mainstem_file)
shapefile(x = splines_mainstem, splines_mainstem_file, overwrite = TRUE)7. 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.
8. Lista de verificação de arquivos
- basin[0-3].asc discretização da bacia para diferentes níveis, com -9999 (NA) atribuído para células fora da área da bacia;
unlist(l_write_files)[1] "../softs/priver/pauliceia/pp_asciigrid_files/basin0.asc"
[2] "../softs/priver/pauliceia/pp_asciigrid_files/basin1.asc"b3class : 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 : basin1
values : 1, 9 (min, max)- mainstem.asc raster com valor 1 para células que identificam rio principal com -9999 (NA) atribuído para células fora da área da bacia e 0 para células dentro da bacia mas que não são do rio principal;
mainstemclass : 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 : mainstem
values : 0, 1 (min, max)mainstem_file[1] "../softs/priver/pauliceia/mainstem.asc"9. Referências
Pfafstetter, Otto., 1989. Classification of hydrographic basins: coding methodology, unpublished manuscript, Departamento Nacional de Obras de Saneamento, August 18, 1989, Rio de Janeiro.
Verdin, K. L., Verdin, J. P., 1999. A topological system for delineation and codification of the Earth’s river basins. J. Hydrol. 218, 1-12.
Tang, Q., Oki, T., Hu, H., 2006. A distributed biosphere hydrological model (DBHM) for large river basin. Ann. J. Hydraul. Eng. JSCE 50, 37-42.
Tang, Q., 2006. A Distributed Biosphere-Hydrological Model for Continental Scale River Basins. The University of Tokyo, Tokyo, Japan, Ph.D. thesis.