DBHM-GIS: PARTE 2
15-11-2014
Universidade Federal de Santa Maria - UFSM
Departamento de Física
Prof. Jônatan Tatsch
Aluno Roilan Hernandez Valdes
1 Preenchimento de depressões do terreno e geração de direções de fluxo
1.1 Rotinas em Fortran do IPH (executáveis)
A seguir serão executadas rotinas em Fortran desenvolvidas por Paz (2008) que utilizam arquivos no formato formato padrão do software Idrisi. Nesse formato, cada imagem raster é composta por dois arquivos de mesmo nome, um arquivo com extensão .rst que contém a informação propriamente dita da imagem (ou seja, os valores referentes a cada pixel), e um arquivo com extensão .rdc, que constitui a documentação da imagem (contém informações sobre o sistema de referência e projeção, número de linhas e de colunas, resolução, vértices, valores máximos e mínimos da imagem, etc). O arquivo .rdc pode ser aberto e editado tanto no próprio Idrisi quanto em um editor de texto qualquer, ou mesmo via alguma rotina específica (trata-se de um arquivo em formato ascii). Mas atenção: dependendo decomo o arquivo .rdc foi gerado e de qual informação se quer editar, é possível ou não utilizar algum processador de texto ou o próprio software Idrisi.
Os executáveis Fortran necessários nas etapas de preparação de dados do MDET estão disponíveis na página web do Prof. Dr. Adriano Paz. Uma cópia do manual MGB-GIS (Paz, 2008) está disponpível na página da disciplina no Moodle. As partes desse texto que descrevem a utilização dos executáveis Fortran são baseadas no manual descrito em Paz (2008). Algumas etapas foram adaptadas para execução no R.
1.1.1 Descrição
Esta etapa consiste em gerar as direções de fluxo de alta resolução (100 m nos casos da bacia mogu e vacacai), ou seja, gerar uma imagem raster onde cada elemento (pixel) contém um código que indica para qual pixel vizinho ocorre a drenagem. Atribui-se uma única direção para cada pixel, que pode assumir um dos oito valores possíveis definidos por convenção de acordo com o software utilizado (Figura 1). A convenção de direção de fluxo do pacote raster do R é a mesma utilizada no ArcGis (Figura 1).
As direções de fluxo são determinadas com a rotina MNTAlta4A.exe. O algoritmo de cálculo segue o algoritmo D8 (eight-deterministic neighbours, Jenson e Domingue (1988)) com algumas variações. O algoritmo D8 consiste basicamente em atribuir a direção de fluxo de um pixel para aquele pixel vizinho conforme a maior declividade (diferença de elevação/distância entre os pixels). A rotina MNTAlta4A.exe emprega um algoritmo que leva em conta ainda um fator de aleatoriedade na busca por direções de fluxo em áreas planas, tal qual adotado por Fairfield e Leymarie (1991), o que têm evitado satisfatoriamente a tendência do algoritmo D8 em criar drenagens excessivamente paralelas.
1.1.2 Aplicação
Crie um diretório nomeado
../softs/Paz2006_routines/MNTAlta/basin, substituindo o basin p.ex.: mogu, vacacai.$ mkdir -p ../softs/Paz2006_routines/MNTAlta/basinCopie os arquivos
../data_base/basin/mdet_basin_laea_int_crop.[rst,rdc]gerados na etapa anterior para o diretório criado../softs/Paz2006_routines/MNTAlta/basin, renomeando-os paraMNT.[rst,rdc].$ cp ../data_base/basin/mdet_basin_laea_int_crop.* ../softs/Paz2006_routines/MNTAlta/basin/Copie o executável
../softs/Paz2006_routines/MNTAlta/MNTAlta.exepara../softs/Paz2006_routines/MNTAlta/basine o arquivo combina.dir também.$ cp ../softs/Paz2006_routines/MNTAlta/MNTAlta4A.exe ../softs/Paz2006_routines/MNTAlta/basin/ $ cp ../softs/Paz2006_routines/MNTAlta/combina.dir ../softs/Paz2006_routines/MNTAlta/basin/Renomeie os arquivos
../softs/Paz2006_routines/MNTAlta/basin/mdet_basin_laea_int_crop.[rst,rdc]para../softs/Paz2006_routines/MNTAlta/basin/MNT.[rst,rdc].$ mv ../softs/Paz2006_routines/MNTAlta/basin/mdet_basin_laea_int_crop.rst ../softs/Paz2006_routines/MNTAlta/basin/MNT.rst $ mv ../softs/Paz2006_routines/MNTAlta/basin/mdet_basin_laea_int_crop.rdc ../softs/Paz2006_routines/MNTAlta/basin/MNT.rdcCertifique-se de o arquivo combina.dir está no diretório
../softs/Paz2006_routines/MNTAlta/basin.Execute a rotina clicando com o botão direito em
MNTAlta4A.exe> open with wine windows …leia as informações na nova tela aberta e tecle
, novamente , digite 1 (para ler MNT.rst), <enter>, <enter>, digite 0 (para não ler mascara.rst) e aguarde alguns minutos …
quando terminar a execução dessa rotina: faça uma cópia do arquivo MNT.rdc e o renomeie para MNTfill.rdc; clique com botão direito do mouse no MNTfill.rdc > open with idrisi32 …
na janela aberta do Idrisi clique em Tools > calculate resolution e depois Tools > Calculate Min/Max; então clique em File > save, feche o arquivo.
clique duas vezes no arquivo MNTfill.rst para visualizá-lo no Idrisi
repita os 3 últimos passos para o arquivo Dir.rst.
1.1.3 Visualização no R
Configurações básicas
Para visualizarmos as saídas da rotina MNTAlta.exe no R precisamos primeiro carregar os pacotes necessários e definir algumas variáveis com a lozalização dos diretórios e dados, similar ao realizado nas etpas 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 <- "hidro2roilan"
#user_name <- "pqgfapergs1"
## 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")
## scipt para gerar plot de drenagem
source("flow_lines.R")
## scipt para gerar plot de drenagem
source("gdal_polygonizeR.R")Para preparação de dados é fundamental definirmos 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.
## defina o nome de sua bacia
#basin_name <- "mogu"
# basin_name <- "vacacai"
basin_name <- "pauliceia"
## mdet da bacia de alta resolução (100 m) gerado na etapa anterior
mdet_alta_file <- gsub("basin", basin_name,"../data_base/basin/mdet_basin_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
laea <- projection(raster(gsub("basin", basin_name,"../data_base/basin/mdet_basin_laea_int.tif")))
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) <- laeaVamos importar também os polígonos de bacias da ANA para visualização no R. Os polígonos das bacias da ANA são definidos para diferentes níveis de discretização, quanto maior o nível, maior o detalhamento das sub-bacias. Para as bacias mogu e vacacai o nível 4 ou 5 é suficiente. Defina o nível de discretização (X) de sub-bacias de seu interesse conforme a escala espacial da bacia de interesse.
## nível de discretização
X <- 5
otto_file <- gsub("X", X, "../data_base/ANA/Ottobacias_NivelX/Ottobacias_NivelX.shp")
otto <- shapefile(otto_file, encoding = "latin1")
## definindo projeção
otto <- spTransform(x = otto, CRSobj = CRS(proj4string(mdet_alta)))
## recortando para o domínio do mdet
otto <- crop(otto, mdet_alta)
## plot sobrepondo as bacias no MDET
plot(mdet_alta, main = "MDET na projeção geográfica")
plot(otto, main = paste0("ottobacias - nível: ",X), add = TRUE)
Visualização dos rasters de saída
Vamos verificar se as direções de fluxo estão com valores plausíveis. Para isso usamos os comandos abaixo para para visualizar os rasters de saída da rotina MNTAlta.exe vamos importá-los para o R usando a função raster do pacote de mesmo nome.
Vamos importar o raster de direção de fluxo gerada com rotinas IPH e definir sua projeção.
fdir_alta_file <- gsub("basin", basin_name,"../softs/Paz2006_routines/MNTAlta/basin/Dir.rst")
fdir_alta <- raster(fdir_alta_file)
## define sist proj.
projection(fdir_alta) <- projection(mdet_alta)
fdir_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. : +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/MNTAlta/pauliceia/Dir.rst
names : Dir
values : 1, 128 (min, max) plot(fdir_alta, main = "Direção de fluxo do MDET de alta resolução (90 m)")
Repetir procedimento para as direções de fluxo.
## importa raster do mdet preenchido
mdet_alta_fill_file <- gsub(pattern = "basin",
replacement = basin_name,
x = "../softs/Paz2006_routines/MNTAlta/basin/MNTfill.rst")
mdet_alta_fill <- raster(mdet_alta_fill_file)
projection(mdet_alta)[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" ## define sist proj.
projection(mdet_alta_fill) <- projection(mdet_alta)
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. : +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/MNTAlta/pauliceia/MNTfill.rst
names : MNTfill
values : 572, 751 (min, max) plot(mdet_alta_fill)
Vamos verificar se as direções de fluxo estão com valores plausíveis e então reescrever o arquivo original para salvar as alterações.
## extraindo valores de direção de fluxo de alta res
fdir_alta_v <- getValues(fdir_alta)
sum(is.na(fdir_alta_v))[1] 0u_dir <- unique(fdir_alta_v)
u_dir[1] 32 128 2 64 16 8 4 1## frequência de ocorrência das direções de fluxo
barplot(prop.table(table(fdir_alta_v, useNA = "always"))*100,
main = "Frequência de ocorrência das direções de fluxo")
## distribuição espacial de NAs
plot(is.na(fdir_alta), main = "Distribuição espacial de dados inválidos (NA)")
## definindo constante para valores faltantes
fdir_alta[is.na(fdir_alta)] <- -9999
## verificando substituição
unique(fdir_alta)[1] 1 2 4 8 16 32 64 128## precisamos atualizar o arquivo com novos valores
writeRaster(x = fdir_alta,
filename = fdir_alta_file,
NAflag = -9999,
format='IDRISI',
overwrite = T,
datatype = "INT2S")class : 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/Dir.rdc
names : layer
values : 1, 128 (min, max)## Comparação entre os mdets
## diferença entre o mdet preenchido e o mdet bruto
mdet_diff <- mdet_alta_fill - mdet_alta
plot(mdet_diff, main = "Diferença entre os MDETs preenchido e o bruto")
# zoom(mdet_diff)2 Determinação das áreas de drenagem acumuladas
2.1 Rotinas TauDEM
A derivação da área de drenagem acumulada será realizado com a ferramenta aread8 do TauDEM (Terrain Analysis Using Digital Elevation Models). A convenção de direção de fluxo adotada pelas rotinas TauDEM é diferente daquela usada nas rotinas do IPH (Figura 2).
Por isso precisamos converter o raster de direções de fluxo ../softs/Paz2006_routines/MNTAlta/basin/Dir.rst para essa convenção usada pelo TauDEM. Mas primeiro vamos criar um sub-diretório com nome da bacia de interesse no diretório do TauDEM.
## definindo local do diretório taudem
taudem_dir <- "../softs/TauDEM-master"
## criando sub-diretório da bacia no diretório do taudem
dir_taudem_basin <- paste0(taudem_dir,"/basin")
dir_taudem_basin <- gsub("basin", basin_name, dir_taudem_basin)
dir_taudem_basin[1] "../softs/TauDEM-master/pauliceia"## cria diretório da bacia se ele não existir
if(!basin_name %in% dir(dirname(dir_taudem_basin))) {
dir.create(path = dir_taudem_basin, showWarnings = T)
}Vamos então reclassificar as direções de fluxo do raster fdir_alta para a convenção de direção de fluxo do TauDEM. O raster resultante será salvo em format .tif e .rst com nome ../softs/TauDEM-master/basin/dir_mntalta_taudem.*. Para verificar a reclassificação, vamos comparar a frequência de ocorrência das direções de fluxo nas duas convenções.
## releitura do raster de dir fluxo no formato idrisi (fda: flux. direc. alta)
fdir_alta <- raster(fdir_alta_file)
fdir_alta_v <- getValues(fdir_alta)
## distribuição das direções de fluxo obtidas com o MNTAlta
round(prop.table(table(fdir_alta_v))*100, 2)fdir_alta_v
1 2 4 8 16 32 64 128
5.81 15.14 10.86 14.07 12.74 21.19 11.48 8.71 ## reclassificando da convenção do iph para a do taudem
fdir_alta_taudem <- recode_flowdir(fd = fdir_alta,
from = "iph",
to = "taudem")## extraindo valores do raster fdir_alta_taudem
fdir_alta_taudem_v <- getValues(fdir_alta_taudem)
## conferindo distribuição das direções de fluxo na nova convenção
## a proporção deve ser mantida entre as direções prop. IPH = 2, TD = 1 e etc.
round(prop.table(table(fdir_alta_v))*100, 2)fdir_alta_v
1 2 4 8 16 32 64 128
5.81 15.14 10.86 14.07 12.74 21.19 11.48 8.71 round(prop.table(table(fdir_alta_taudem_v))*100, 2)fdir_alta_taudem_v
1 2 3 4 5 6 7 8
15.14 5.81 8.71 11.48 21.19 12.74 14.07 10.86 A seguir, vamos salvar os rasters de direção reclassificados em format .tif e .rst.
## nome do arquivo raster de dir. fluxo gerado com o MNTAlta em formato tid para o taudem
fdir_alta_taudem_file <- paste0(dir_taudem_basin,"/dir_mntalta_taudem.tif")
fdir_alta_taudem_file <- gsub("basin", basin_name, fdir_alta_taudem_file)
fdir_alta_taudem_file[1] "../softs/TauDEM-master/pauliceia/dir_mntalta_taudem.tif" ## fdir_alta_taudem em dormato SpatialGridDataFrame (pacote sp) para escrita em formato requerido pelo taudem
fdir_alta_taudem_sgdf <- as(fdir_alta_taudem, "SpatialGridDataFrame")
## salvar raster de dir. fluxo na convenção do taudem em formato tif (requerido pelo TauDEM)
writeGDAL(dataset = fdir_alta_taudem_sgdf,
fname = fdir_alta_taudem_file,
drivername = "GTiff",
type="Int16",
mvFlag = -32768)
## salvar raster de dir. fluxo na convenção do taudem em formato .rst (se quiser converter para tif com idrisi)
writeRaster(x = fdir_alta_taudem,
filename = gsub("\\.tif","\\.rst",fdir_alta_taudem_file),
NAflag = -9999,
format='IDRISI',
overwrite = T,
datatype = "INT2S")class : 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/TauDEM-master/pauliceia/dir_mntalta_taudem.rdc
names : layer
values : 1, 8 (min, max)A rotina aread8 do TauDEM será aplicada usando o raster de direções de fluxo ../softs/Paz2006_routines/MNTAlta/mogu/dir_mntalta_taudem.tif (gerado com a rotina MNTAlta.exe) para obter o raster de área de drenagem acumulada (arquivo que será nomeado como ../softs/TauDEM-master/basin/area_taudem.tif).
area_taudem_alta_file <- paste0(dir_taudem_basin,"/area_taudem_alta.tif")
## Area acumulada com TauDEM
cmd <- paste0("mpiexec -n 1 aread8 -p ", fdir_alta_taudem_file," -ad8 ", area_taudem_alta_file)
## atualizando comando de acordo com o sub-diretório da bacia
cmd <- gsub("basin", basin_name, cmd)
cmd[1] "mpiexec -n 1 aread8 -p ../softs/TauDEM-master/pauliceia/dir_mntalta_taudem.tif -ad8 ../softs/TauDEM-master/pauliceia/area_taudem_alta.tif"system(cmd)2.2 Conversão da saída do TauDEM para o formato raster do Idrisi
A rotina AreaAcu.exe poderia ser usada para geração da área de drenagem acumulada. Entretanto ela tem a desvantagem do longo tempo de execução, por isso como alternativa usamos a rotina do TauDEM que tem um tempo de execeução muito menor para a derivação da área de drenagem acumulada. A saída de área acumulada do TauDEM é em termos de células acumuladas e para uso na rotina de upscaling de direções de fluxo (dirfluxo4.exe) precisamos converter o arquivo ../softs/TauDEM-master/area_taudem.tif para área acumulada em km2, deixando-o assim consistente com o arquivo ../softs/Paz2006_routines/AreaAcu/basin/AreaAcu.rst que poderia ser gerado pela rotina AreaAcu.exe do IPH.
Além disso, precisamos converter o raster do formato .tif para rst.
## importando area acumulada TauDEM gerado do mdet preenchido pelo TauDEM
area_taudem_alta <- raster(area_taudem_alta_file)
## atualiza estatisticas do raster
plot(area_taudem_alta, col = tim.colors(32), main = "Área de drenagem acumulada (TauDEM)")
## área de cada délula em km2
area_cel_km2 <- prod(res(area_taudem_alta)/1000) # 0.1 * 0.1
area_cel_km2[1] 0.01 ## multiplicar pela area de cada célula (km2) para compatibilidade com area do iph
area_taudem_alta_km2 <- area_taudem_alta * area_cel_km2
area_taudem_alta_km2class : 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 : area_taudem_alta
values : 0.01, 25.02 (min, max) plot(area_taudem_alta_km2, col = tim.colors(32), main = "Área de drenagem acumulada (TauDEM) em Km²")
## plot em escala log
# plot(log10(area_taudem_km2), col = tim.colors(32))
## gerar nome do arquivo raster de área acum. em km2 no formato do Idrisi
area_taudem_alta_km2_file <- gsub(pattern = "area_taudem_alta.tif",
replacement = "area_taudem_alta_km2.rst",
x = area_taudem_alta_file)
writeRaster(x = area_taudem_alta_km2,
filename = area_taudem_alta_km2_file,
NAflag = -9999,
format='IDRISI',
overwrite = T,
datatype = "FLT8S")class : 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/TauDEM-master/pauliceia/area_taudem_alta_km2.rdc
names : layer
values : 0.01, 25.02 (min, max)Vamos visualizar a rede de drenagem gerada, mas para isso o valor da variável limiar deve ser definida para que a rede de drenagem seja gerada. Um valor de 5 vezes a resolução espacial do raster é assumido como default visualização.
## raster de drenagem, assumindo limiar de 100 km2 para formação de rios, conforme anteriormente
fator <- 1
limiar <- fator*5*sqrt(prod(res(fdir_alta)[1]/1000)) # km2
limiar[1] 1.581139 limiar <- 0.5 # km2 (para bacia paulicéia)
dren <- fdir_alta
dren[area_taudem_alta_km2 < limiar] <- NA
plot(dren, col = "blue", main = paste0("Rede de drenagem (área >", limiar, "km²"))
2.3 Geração de polígono da bacia de alta resolução
O passo principal para determinar o contorno da bacia de alta resolução é localizar o exutório da bacia. Siga os passos a seguir para encontrá-lo e salvá-lo em um arquivo shapefile.
## plot área de drenagem
plot(area_taudem_alta_km2, col = tim.colors(32), main = "Área de drenagem (TauDEM) em Km²")
## linhas de drenagem
flow_lines(dren, length = 0.001, col = "yellow", lwd = 3)
## se a bacia de estudo da bacia for microescala é possível que o nível de discretização da ANA
## não apresente bacias para visualização
# plot(otto, border = 2, add = TRUE)
shapefile_exut_file <- paste0(dir_taudem_basin,"/exutorio_alta.shp")
sp_xy_exut <- shapefile(x = shapefile_exut_file)
projection(sp_xy_exut) <- projection(mdet_alta)
plot(sp_xy_exut, add = T,lwd = 6, col = 2)
sp_xy_exutclass : SpatialPointsDataFrame
features : 1
extent : -951.0219, -951.0219, -1075.847, -1075.847 (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
variables : 1
names : ID
min values : 1
max values : 1 # shapefile_exut_jdtatsch <- paste0("../data_base/data_pauliceia","/exutorio_alta.shp")
# sp_xy_exut_tatsch <- shapefile(x = shapefile_exut_jdtatsch)
# plot(sp_xy_exut_tatsch, add = T,lwd = 6, col = 6, lty = 5)O polígono com contorno das bacias da ANA é apenas para referência aproximada da bacia, a localização do exutório deve ser realizada considerando a rede drenagem e área acumulada. O pixel mais à jusante da bacia e que possui maior área de drenagem, antes de confluência com outro rio, corresponde ao exutório da bacia. Por isso vamos dar um zoom na área de drenagem e encontrar o ponto com tais características.
## selecione um retângulo em torno do exutório da bacia
zoom(area_taudem_alta_km2, col = tim.colors(32))
## clique no ponto exutório para extrair as coordenadas xy
## talvez seja necessário dar mais zoom (repita o zoom quanto necessário for necessário!)
xy_exut <- locator(n=1) A região na qual deve ser feito zoom para localização do exutório deve ter aspecto semelhante ao mostrado nas Figura 5.
Na Figura 6 o ponto exutório da bacia (Vacacaí) é indicado por uma cruz.
## converte coordenadas para SpatialPoints
sp_xy_exut <- SpatialPoints(xy_exut)
## projeta para laea
proj4string(sp_xy_exut) <- laea
plot(sp_xy_exut, add = TRUE)
## nome do arquivo shapefile como localizacao do exutório
shapefile_exut_file <- paste0(dir_taudem_basin,"/exutorio_alta.shp")
## salva shapefile
shapefile(x = sp_xy_exut, filename = shapefile_exut_file, overwrite = TRUE)
## mostra informações do shapefile do ponto exutório
sp_xy_exutCom a localização do ponto exutório podemos obter a célula do raster de alta de resolução (area_taudem_km2) e então a área de drenagem da bacia.
## determina a célula do exutório no raster de alta resolução
cell_exut_alta <- cellFromXY(object = area_taudem_alta_km2,
xy = sp_xy_exut)
cell_exut_alta[1] 2840## área acumulada da bacia naquela célula (km2)
area_bacia <- raster::extract(x = area_taudem_alta_km2,
y = cell_exut_alta)
area_bacia
11.94 Vamos gerar a área de drenagem acumulada fornecendo o arquivo shapefile do ponto exutório da bacia de alta resolução (100 m), fazendo com que a rotina gere um raster de área somente para as células que drenam para aquele ponto.
area_taudem_alta_exutorio_file <- paste0(gsub("\\.tif","",area_taudem_alta_file), "_exutorio.tif")
## Area acumulada com TauDEM
cmd <- paste0("mpiexec -n 4 aread8 -p ",
fdir_alta_taudem_file,
" -ad8 ",
area_taudem_alta_exutorio_file,
" -o ",
shapefile_exut_file)
## atualizando comando de acordo com o sub-diretório da bacia
system(cmd)Importanto área de drenagem de alta resolução determinada somente para a região a montante do ponto exutório para produzirmos uma mascara da bacia de alta resolução.
## raster de área de drenagem para a bacia formada pelo exutório
area_taudem_alta_exutorio <- raster(x = area_taudem_alta_exutorio_file)
plot(area_taudem_alta_exutorio,
col = tim.colors(32),
main = "Área de drenagem formada pelo exutório")
## mascara da bacia de alta resolução
## constrói raster com propiedades do raster area_taudem_exutorio
mascara_bacia_alta <- raster(x = area_taudem_alta_exutorio)
## atribui valores das células: 1 dentro, NA fora da bacia
mascara_bacia_alta[!is.na(area_taudem_alta_exutorio)] <- 1Convertendo de raster para polígono, para então gerar um shapefile do polígono da bacia de alta resolução. Primeiro verifique se está instalado o Python. No R digite Sys.which("gdal_polygonize.py"). Se retornar uma sctring vazia, ou o script gdal_polygonize.py não existe no seu sistema, ou o caminho para o diretório dele está faltandona dus variável PATH. Nesse último caso, podemos informar o caminho através do argumento pypath aceito pela função usada abaixo para especificar o caminho, ou modificar sua variável PATH.
# Sys.which("gdal_polygonize.py")
## aplica função polygonize (script gdal_polygonizeR.R) para gerar polígono da bacia de alta resolução
mascara_bacia_alta_pol <- gdal_polygonizeR(x = mascara_bacia_alta,
pypath = "/usr/bin/gdal_polygonize.py")
## define projeção
proj4string(mascara_bacia_alta_pol) <- laea
## informações básicas do SpatialPolygonsDataFrame
class(mascara_bacia_alta_pol)[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"summary(mascara_bacia_alta_pol)Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -1401.520 2298.480
y -1622.855 2877.145
Is projected: TRUE
proj4string :
[+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 attributes:
Min. 1st Qu. Median Mean 3rd Qu. Max.
1 1 1 1 1 1 ## gráfico
plot(mascara_bacia_alta_pol,
col = "gray",
border = "black",
lwd = 4)
title(paste0("Polígono da bacia gerado do MDET \n de alta resolução ",
sqrt(prod(res(mascara_bacia_alta))),
" m"))
plot(otto, border = 2, add = TRUE)
Vamos salvar o shapefile do contorno da bacia derivado a partir do MDET de alta de resolução.
## nome do arquivo shapefile
mascara_bacia_alta_pol_file <- gsub(pattern = basename(area_taudem_alta_exutorio_file),
replacement = "mascara_bacia_alta_pol.shp",
x = area_taudem_alta_exutorio_file)
## salva polígono em shapefile
shapefile(x = mascara_bacia_alta_pol,
filename = mascara_bacia_alta_pol_file,
overwrite = TRUE)3 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.
4 Lista de arquivos gerados
- MNTfill.rst MDET modificado pelo algoritmo abrangendo o domínio espacial da bacia, na projeção cartográfica LAEA (formato Idrisi/integer/binary)
mdet_alta_fill_file[1] "../softs/Paz2006_routines/MNTAlta/pauliceia/MNTfill.rst"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. : +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/MNTAlta/pauliceia/MNTfill.rst
names : MNTfill
values : 572, 751 (min, max)- Dir.rst arquivo raster (formato Idrisi/integer/binary) de direções de fluxo de alta resolução (100 m), abrangendo o domínio espacial da bacia onde cada pixel contém um valor referente a sua direção seguindo codificação do IPH (Figura 2);
fdir_alta_file[1] "../softs/Paz2006_routines/MNTAlta/pauliceia/Dir.rst"fdir_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/MNTAlta/pauliceia/Dir.rst
names : Dir
values : 1, 128 (min, max)- area_taudem_alta_km2.rst arquivo raster no formato tif com área de drenagem acumulada em Km² para todo domínio espacial da bacia (sem máscara aplicada);
area_taudem_alta_km2_file[1] "../softs/TauDEM-master/pauliceia/area_taudem_alta_km2.rst"area_taudem_alta_km2class : 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 : area_taudem_alta
values : 0.01, 25.02 (min, max)5. Referências
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.