DBHM-GIS: PARTE 1

---

Universidade Federal de Santa Maria - UFSM

Departamento de Física

Prof. Jônatan Tatsch

Aluno Roilan Hernandez Valdes

---

1. Introdução

1.1. DBHM e documentação

O Modelo Hidrológico Distribuído da Biosfera (DBHM) é espacialmente distribuído que integra processos hidrológicos e processos de transferência solo-vegetação-atmosfera (SVAT) na escala bacia hidrográfica. O DBHM foi desenvolvido pelos Dr. Qiuhong Tang, Taikan Oki and Shinjiro Kanae na Universidade de Tokyo. A referência para citação e agradecimentos é:

• Tang, Q., 2006. A Distributed Biosphere-Hydrological Model for Continental Scale River Basins. The University of Tokyo, Tokyo, Japan, Ph.D. thesis.

Diversos artigos científicos importantes relacionadas ao modelo são:

• Tang, Q., T. Oki, and S. Kanae, 2006. A distributed biosphere hydrological model (DBHM) for large river basin. Ann. J. Hydraul. Eng. JSCE, 50, 37–42.

• Tang, Q., Oki, T., Kanae, S., Hu, H., 2008. Hydrological Cycles Change in the Yellow River Basin During the Last Half of the 20th Century., J. Climate, 21(8), 1790-1806. DOI: 10.1175/2007JCLI1854.1.

• Tang, Q., Oki, T., Kanae, S., Hu, H., 2007. The influence of precipitation variability and partial irrigation within grid cells on a hydrological simulation. J. Hydromet., 8(3), 499-512. DOI: 10.1175/JHM589.1.

• Tang, Q., Oki, T., 2007. Daily NDVI relationship to cloud cover. J. Appl. Meteorol., 46(3), 377-387. DOI: 10.1175/JAM2468.1.

• Tang, Q., Oki, T., Kanae, S., Hu, H., 2007. A spatial analysis of hydro-climatic and vegetation condition trends in the Yellow River Basin. Hydrol. Process., 22, 451-458. DOI: 10.1002/hyp.6624.

O código da estimativa de radiação solar incidente foi desenvolvido por:

• Yang, K., and T. Koike (2005), A general model to estimate hourly and daily solar radiation for hydrological studies, Water Resour. Res., 41,W10403, doi:10.1029/2005WR003976.

A referência para o modelo de superfície terrestre empregado (versão revisada do Modelo Simples da Biosfera, SiB2) é:

• Sellers, P. J., and Coauthors, 1996. A revised land surface parameterization (SiB2) for atmospheric GCMs. Part I: Model formulation. J. Climate, 9, 676–705.

O seguinte endereço disponibiliza para download o código do DBHM, incluindo um tutorial e informações gerais.

1.2. DBHM-GIS

O DBHM requer diversas informações de entrada para sua aplicação ao domínio espacial da bacia hidrográfica de interesse. Entre os dados de sensoriamento remoto inclui-se o modelo digital de elevação do terreno (MDET).

Neste tutorial apresenta-se uma descrição detalhada sobre como gerar os dados de entrada estáticos do DBHM, em particular os que caracterizam a geomorfologia da bacia hidrográfica de interesse.

Na etapa anterior (arquivo dbhm_gis_parte0.html) foi realizado o download do MDET do SRTM e o recorte do MDET para uma região em torno da bacia hidrográfica de interesse.

Esta etapa consiste na definição do domínio espacial da bacia, com um número de colunas e linhas adequado para que no procedimento de upscaling de direções de fluxo da alta resolução (100 m) para a baixa resolução (1000 m) não ocorra inconsistências entre os domínios espaciais. Diversas operações de geoprocessamento serão necessárias para produzir os dados de entrada estáticos do DBHM, tais como o MDET preenchido (sem depressões), direções de fluxo, áreas de drenagem acumuladas, distância ao exutório, declividade, divisão da bacia em sub-bacias e etc. Todo esse processamento será composto por diversas partes e o conjunto formado é denominado DBHM-GIS.

1.3. Dados geomorfológicos de entrada do DBHM

As características geomorfológicas da bacia são determinadas por meio de técnicas de geoprocessamento a partir do MDET. Os dados são armazenados em arquivos texto, formato ARC/INFO ASCII GRID. Esses arquivos serão salvos no diretório GEO_INPUT_DATA (Figura 1).

Os arquivos desse diretório que serão gerados nesse documento referem-se aos arquivos geomorfológicos da bacia, os quais incluem-se:

• wsdem.asc MDET da bacia com altitude em m e com o valor -9999 (NA) atribuído para células fora da área da bacia;

• wsdemirr.asc MDET da bacia com altitude em m, incluindo células de áreas irrigadas que podem estar fora da bacia hidrográfica, o valor -9999 (NA) é atribuído para células fora da área da bacia e áreas não irrigadas;

• dir.asc códigos de direções de fluxo (rede de drenagem) seguindo a convenção do ArcGIS (Figura 2):

• acc.asc área de drenagem acumulada (número de células a montante da célula que escoam para ela);

• frac.asc fração de área dentro da bacia (%);

• fracirr.asc fração de área dentro da bacia (%) incluindo áreas irrigadas;

• gridarea.asc área de cada célula (km2);

• slope.asc declividade do terreno (%);

• slpsib.asc declividade do terreno (radianos) para o SiB2;

• demhdr.asc cabeçalho com as informações do domínio espacial da bacia hidrográfica;

• dis.asc distância a foz da bacia (km);

• rivlen.asc comprimento dos trechos de rio que conectam células vizinhas;

2. Pré-requisitos

2.1. Disponibilidade de dados de chuva e vazão

Antes de qualquer passo a frente é imprescindível verificar se há dados hidrometeorológicos disponíveis para bacia selecionada. O fato de existir estações de medida localizadas no domínio da bacia hidrográfica não significa que existam dados disponíveis para o seu período de interesse.

Então, faça o download dos dados fluviométricos (vazão) através do site hidroweb especificando o código da estação de interesse. O código pode ser obtido do arquivo /home/user/DBHM/data_process/geo_data/data_base/basin/flu6(8).kmz para o caso da bacia do Rio Mogi-Guaçú (ou Vacacaí), mogu (ou vacacai).

Após o download abra os arquivos das estações fluviométricas e no “olhômetro” verifique se o período de dados disponíveis coincide seu período de interesse. Em uma etapa futura, essa verificação dos dados fluviométricos será mais detalhada.

Para encontrar os códigos das estações pluviométricas faça a conversão do arquivo shapefile (pontos das estações pluviométricas da ANA) para KML ou KMZ, conforme feito no script /home/user/DBHM/data_process/geo_data/scripts/1_setRiverBasinDomain.R para as estações pluviométricas.

2.2. Estrutura de diretórios e organização dos dados

Para execução desse tutorial assume-se que o Sistema Operacional é Linux, distribuição Ubuntu. São necessários os seguintes dados organizados segundo a estrutura de diretórios pré-definida. Note a inclusão do sub-diretório basin (que pode ser nomeado como mogu ou vacacai). O arquivo essencial para a realização dessa etapa é mdet_basin.tif gerado na etapa anterior (dbhm_gis_parte0.html).

Os arquivos necessários para realização dessa etapa e os que serão gerados nela são identificados por diferentes símbolos na Figura 1, conforme indicado na legenda.

2.3. Instalação de softwares necessários

Alguns procedimentos que serão realizados nessa etapa de preparação dos dados estáticos (invariáveis no tempo) de entrada para o DBHM exigem a aplicação de executáveis gerados por rotinas em Fortran no formato do sistema operacional (SO) Windows (arquivos .exe). Essa rotinas foram desenvolvidas por Paz (2008) e trabalham com um formato de arquivos raster no formato padrão do software Idrisi. Tal software só é disponível no SO Windows. Por essa razão instalaremos o software wine que permite a execução de alguns softwares desenvolvidos exclusivamente para Windows dentro do Linux.

• wine

  $ sudo apt-get install wine

  pressione a tecla “super” do teclado (ou tecla com símbolo windows) e clique em Configure Wine > ok.

  • Idrisi

    • Caso seu linux esteja no idioma Português é necessário alterá-lo para o idioma Inglês para que o Idrisi funcione adequadamente

      • Vá em congigurações do sistema > Suporte de Linguagem > Instale o suporte de linguagem (clique install) > configure o idioma para inglês como primeira opção > aplicar a todo sistema;
      • Na aba configurações regionais > defina inglês como primeira opção > aplicar a todo sistema
      • Reinicie sua máquina para que as alterações sejam aplicadas

  Acesse o diretório /home/user/DBHM/GIS/softs/IdrisiKilimanjaro clique com botão direito em setup.exe > Open with Wine > next > I Agree > next > (não preencha) next > next > next > next > next > finish

  Para instalar uma atualização clique com botão direito em /home/user/DBHM/GIS/softs/IdrisiKilimanjaro/atualizacao/setup_sapatch1401.exe e > next > next > next > finish.

  Copie o arquivo de licença para o diretório de instalação do Idrisi:

  $ cp /home/user/DBHM/GIS/softs/IdrisiKilimanjaro/license_not_comercial/lservrc /home/user/.wine/drive_c/IDRISI\ Kilimanjaro/

  Instale a segunda atualização: clique com botão direito em /home/user/DBHM/GIS/softs/IdrisiKilimanjaro/atualizacao/setup_sapatch1402.exe e > next > next > next > finish.

  Abra o dash (segure a tecla “super”) e digite Idrisi, clique no icone para abrir o Idrisi.

• Instalação de pacotes do linux necessários para uso de alguns pacotes no R:

  $ sudo apt-get install libwxgtk-media3.0-0
  $ sudo apt-get install libwxgtk-media3.0-dev
  $ sudo apt-get install libwxgtk3.0-dev
  $ sudo apt-get install gdal-bin 
  $ sudo apt-get install proj-bin
  $ sudo apt-get install tk-dev

• Instalação da biblioteca python-gdal

  $ sudo apt-get install python-gdal python3-gdal

• Instalação das rotinas TauDEM (versão 5.1.2)

  O download é opcional pois o arquivo zip encontra-se já disponível no diretório $HOME/DBHM/data_process/geo_data/softs/TauDEM-master.zip

  $ cd $HOME/DBHM/data_process/geo_data/softs/ 
  $ wget https://github.com/bmfekete/TauDEM/archive/master.zip
  $ mv master.zip master.zip TauDEM-master.zip
  $ unzip TauDem-master.zip; cd TauDem-master

  Vamos instalar o mpi (mpich2) eo cake necessários para compilação do TauDEM:

  $ sudo apt-get install mpich2 cmake
  $ cd src

  Abra o arquivo linearpart.h e adicione depois da linha

  #include "mpi.h"

  Adicione uma nova linha com

  #include <stdint.h>

  Salve as mudanças e feche o arquivo. Abra o arquivo tiffIO.h, encontre a linha #include "stdint.h" e substitua as aspas (“”) por <>, então você obterá:

  #include <stdint.h>

  Salve as mudanças e feche o arquivo. Crie um diretório chamado build e entre nele

  $ mkdir build
  $ cd build
  $ cp -r ../../shape ../
  $ CXX=mpicxx cmake -DCMAKE_INSTALL_PREFIX=/usr/local ..
  $ make

  Finalmente para instalar o TauDEM em /usr/local/bin, rode:

  $ sudo make install

  Teste o TauDEM

  $ cd $HOME/DBHM/data_process/geo_data/data_base/basin
  $ mpiexec pitremove -z srtm_26_18.tif -fel srtm_26_18.tif_fel.tif

• Pacotes no R:

## pacotes para manipulação de dados espaciais
install.packages("sp", dep = TRUE)
install.packages("maptools", dep = TRUE)
install.packages("rgdal", dep = TRUE)
install.packages("RSAGA", dep = TRUE)
install.packages("adehabitat", dep = TRUE)
## pacotes para algumas funções gráficas
install.packages("fields", dep = TRUE)
install.packages("sfsmisc", dep = TRUE)
install.packages("lattice", dep = TRUE)
install.packages("rasterVis", dep = TRUE)
## pacote para manipulação de strings
install.packages("stringr", dep = TRUE)
install.packages("doParallel", dep = TRUE)

3. Reprojeção do MDET da bacia para a Sistema de Projeção espacial do DBHM

3.1. Projeções Cartográficas

Projeções cartográficas são tentativas de representar a superfície terrestre (ou uma porção dela) sobre uma superfície plana. Algumas distorções na direção, escala, distância e área do mapa são inevitáveis neste processo. Dependendo do tipo de projeção as distorções são reduzidas em algumas destas propriedades ao custo do aumento dos erros em outras. Algumas projeções buscam distribuir o erro da transformação de forma moderada entre todas essas propriedades. A maioria das projeções cartográficas pode ser pelo menos parcialmente visualizadas geometricamente como sendo projetadas em três superfícies: o cilindro, o cone e o plano. A projeção cartográfica cilíndrica pode ser visualizada geometricamente como um cilindro envolvendo a Terra: após os pontos do mapa são projetados, o cilindro é então cortado longitudinalmente e desenrolado como uma superfície plana. A Projeção cônica pode ser visualizada como se um cone fosse colocado sobre a Terra com o pico do cone ao longo do eixo da terra e a superfície do cone em contato com a Terra ao longo de uma faixa de latitude constante. A projeção azimutal ou plana pode ser visualizada como um mapa gerado a partir de um plano tangente a superfície terrestre. Se o plano é tangente ao pólo Norte ou Sul, a projeção é chamada polar.

3.2. Projeção Cartográfica Azimutal de área equivalente de Lambert

A projeção cartográfica usada no DBHM é a projeção cartográfica Azimutal de Área Equivalente de Lambert (LAEA). Ela preserva a área constante em toda superfície do mapa e representa realisticamente as direções a partir de seu centro de projeção, mas as deformações na escala aumentam com o distanciamento do ponto central da projeção. A escala diminui com a distância do centro ao longo dos raios e aumenta com a distância do centro em uma direção perpendicular ao raios. A base de dados hidrográficos globais HYDRO1k usa esta projeção. Portanto para usar dados em projeções diferentes (p.ex.: a geográfica também conhecida por longlat ou latlon) no DBHM é necessário projetá-los na projeção LAEA. A projeção LAEA tem as seguintes equações de transformação para o sistema de referência geográfica:

Equações

Os parâmetros da projeção LAEA variam por continente. Os parâmetros para o continente sul americano são 60º W para longitude de origem e 15ºS para latitude de origem e o raio da esfera de influência é de 6370997 m. Para as bacias hidrográficas definiremos esses dois parâmetros como as cooredenadas do centróide da bacia.

3.3. Reamostragem de dados espaciais

Reamostragem é o processo de determinar novos valores para os pixels em uma grade de saída que resulta da aplicação de alguma transformação geométrica para a grade de entrada. A grade de entrada pode estar em sistema de coordenadas diferente, numa resolução diferente, ou pode ser rotaciononado em relação a grade de saída. A reamostragem entre a Projeção Interrompida Homolosina de Goode e a LAEA está inclusa no DBHM. Entretanto para usar o MDET do sensor SRTM com resolução de 90 m é essencial fazer a reamostragem dos dados do sistema de referência geográfica e o elipsóide de referência WSGS 1984 para a projeção LAEA.

3.4 Configurações básicas

Vamos começar carregando os pacotes necessários.

## 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 

A seguir definimos a opção de barra de progresso do pacote raster como tipo texto (“=====”), o nome do usuário e o diretório de trabalho (WD). O WD é definido no diretório /home/user/DBHM/data_process/geo_data/scripts. Nesse diretório encontram-se scripts em R gerados para facilitar a realização da preparação de dados de entrada estáticos do DBHM. Em geral, esses scripts contém funções que requerem a definição de seus argumentos. Esses argumentos podem ter valores default. Os scripts ficam armazenados no diretório script_dir, definido logo a seguir. Esses scripts devem ser carregados para que as funções neles definidas fiquem disponíveis na sessão do R em uso.

## 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("adjust_extent.R")
## script com função para recorte de raster baseado linhas e colunas
source("cropByRowCol.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. No diretório da base dados (data_base) estão os arquivos do MDET da bacia e os polígonos das OttoBacias Brasileiras (para diferentes níveis de divisão) disponíveis no website da ANA. Esses polígonos são utilizados para auxiliar na localização do domínio da bacia.

## defina o nome de sua bacia
basin_name <- "pauliceia"
# basin_name <- "mogu"
#basin_name <- "vacacai"
## diretório com base de dados
basin_db_dir <- gsub("basin", basin_name, "../data_base/basin")
basin_db_dir

[1] "../data_base/pauliceia"

dir(basin_db_dir)

[1] "mdet_pauliceia.asc" "mdet_pauliceia.tif"

Vamos importar os arquivos do MDET da bacia (vacacaí ou mogu) e o polígono da bacias da ANA para visualização no R.

## caminho completo até o arquivo tif da bacia
file_name <- gsub(pattern = "basin",
                  replacement = basin_name, 
                  x = paste0(basin_db_dir,"/mdet_basin.tif"))
file_name

[1] "../data_base/pauliceia/mdet_pauliceia.tif"

## modelo digital de elevação do terreno do SRTM
mdet <- raster(file_name)
mdet

class       : RasterLayer 
dimensions  : 85, 75, 6375  (nrow, ncol, ncell)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : -47.66, -47.5975, -21.67417, -21.60333  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/data_base/pauliceia/mdet_pauliceia.tif 
names       : mdet_pauliceia 
values      : 566, 752  (min, max)

plot(mdet, main = "MDET na projeção geográfica")

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)))
## recortando para o domínio do mdet
otto <- crop(otto, mdet)
## plot sobrepondo as bacias no MDET
plot(mdet, main = "MDET na projeção geográfica")
plot(otto, main = paste0("ottobacias - nível: ",X), add = TRUE)

As vezes podem ocorrer falhas de dados, valores inválidos ou dados inconsistentes que por algum motivo foram removidos. Vamos verificar esse é o caso para o domínio da bacia de estudo.

## há pixels inválidas ou faltantes?
tot_miss <- sum(is.na(getValues(mdet)))
tot_miss

[1] 0

3.5. Projeção do MDET para projeção cartográfica do DBHM

Um atributo espacial básico para localização da bacia é o centróide do seu domínio espacial.

## usando coordenadas próximas ao centroide da bacia
coords <- xyFromCell(object = mdet, cell = 1:ncell(mdet))
## long central da região
lon_med <- mean(coords[, "x"])
## lat central da região
lat_med <- mean(coords[, "y"])
lon_med; lat_med

[1] -47.62875

[1] -21.63875

## centroide do dominio
# centroid_ll <- round(c(lon_med, lat_med))
centroid_ll <- round(c(lon_med, lat_med),2)
 centroid_ll_sp <- SpatialPoints(data.frame(lon = centroid_ll[1], lat = centroid_ll[2]))
  projection(centroid_ll_sp) <- projection(mdet)
## removendo objetos que não serão mais utilizados (limpando memória) 
rm(coords)
centroid_ll_sp

class       : SpatialPoints 
features    : 1 
extent      : -47.63, -47.63, -21.64, -21.64  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

No chunck (pedaço de cógigo ou script) abaixo faremos a projeção do MDET do sistema longlat para laea (Lambert Azimutal Equal Area - Sistema de Coordenadas usada no DBHM). Se o arquivo projetado em LAEA já estiver disponível (i.e. gerado em uma rodada anterior) ele será carregado ao invés de executarmos a mudança de projeção, pois tal procedimentos pode levar alguns minutos dependendo do tamanho do domínio espacial.

## projeção do mdet
projection(x = mdet)

[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

## projecao do DBHM
## para América do Sul +lat_0=-15 +lon_0=-60
# laea <- "+proj=laea +lat_0=-15 +lon_0=-60 +x_0=0 +y_0=0 +ellps=sphere"
## definindo projecao baseado no centróide do domínio da bacia
laea <- "+proj=laea +lat_0=LAT +lon_0=LON +x_0=0 +y_0=0 +ellps=sphere"
laea <- gsub("LON", as.character(centroid_ll[1]), laea)
laea <- gsub("LAT", as.character(centroid_ll[2]), laea)
laea

[1] "+proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +ellps=sphere"

## projetando bacias da ana para mesma projeção do DBHM
otto_laea <- spTransform(x = otto,  CRSobj = CRS(laea))
## projetando a extent do mdet (proj longlat) para a projeção usada no DBHM: 
##  Lambert Azimutal Equal Area (laea)
mdet_laea <- projectExtent(object = mdet, crs = laea)
mdet_laea

class       : RasterLayer 
dimensions  : 85, 75, 6375  (nrow, ncol, ncell)
resolution  : 86.15333, 92.66642  (x, y)
extent      : -3101.52, 3359.98, -3799.501, 4077.145  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=-21.64 +lon_0=-47.63 +x_0=0 +y_0=0 +ellps=sphere 

## definir resolucao para reamostragem
res(mdet_laea) <- 100
mdet_laea

class       : 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 +ellps=sphere 

## nome do mdet na proj laea, valores inteiros
mdet_laea_filename_i <- gsub(pattern = "basin", 
                             replacement = basin_name, 
                             x = paste0(basin_db_dir,"/mdet_basin_laea_int.tif")
                             )
## projetando e salvando em arquivo tif 
## Se tiver sido criado o arquivo mdet_laea_filename, numa rodada anterior por exemplo
# if(file.exists(mdet_laea_filename_i)){
#   ## le arquivo tif criado
#   mdet_laea_i <- raster(mdet_laea_filename_i)
# ## se não
# } else {
  ## faz a projeção do mdet para laea
  ## reprojetando de long lat para laea e salvando arquivo em formato tif
  cat("pode levar algns minutos ... dependendo do domínio e da resolução espacial", "\n")

pode levar algns minutos ... dependendo do domínio e da resolução espacial 

    ## projetando raster com vizinho mais próximo (resulta em valores integer)
    mdet_laea_i <- projectRaster(from = mdet,
                                 to = mdet_laea,
                                 method = "ngb",
                                 filename = mdet_laea_filename_i,
                                 overwrite = TRUE)  
# } # end if

Vamos visualizar do MDET na projeção LAEA e a frequência de ocorrência de altitudes na região da bacia.

## gráfico do mdet
plot(mdet_laea_i, main = "MDET na Proj. Cartográfica\n Lambert Azimutal de área equivalente (LAEA)")

## extraindo valores de altitude do raster mdet
vals_mdet_laea_i <- getValues(mdet_laea_i)
## histograma da altitude do domínio espacial
hist(vals_mdet_laea_i[!is.na(vals_mdet_laea_i)], 
     xlab = "altitude (m)", 
     ylab = "Freq. ocorrência", 
     main = "Histograma da elevação do terreno (m)")

## Função densidade acumulada da altitude do domínio espacial
plot(ecdf(vals_mdet_laea_i[!is.na(vals_mdet_laea_i)]), 
     lwd = 4, 
     xlab = "altitude (m)", 
     ylab = "Frequencia Acumulada",
     main = "Função distribuição acumulada (ecdf)")

Para concluir essa parte inicial vamos salvar o MDET no formato do software Idrisi. São gerados dois arquivos: um com a extensão .rst (contendo os dados em formato binário) e outro .rdc com as informações do domínio espacial.

mdet_laea_idrisi_file <- gsub("tif","rst", mdet_laea_filename_i)
writeRaster(x = mdet_laea_i, 
            filename = mdet_laea_idrisi_file, 
            NAflag = -9999,
            format ='IDRISI', 
            overwrite = T,
            datatype = "INT2S")

  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |===========                                                      |  17%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |================================                                 |  50%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |======================================================           |  83%
  |                                                                       
  |=================================================================| 100%

class       : 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. : NA 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/data_base/pauliceia/mdet_pauliceia_laea_int.rdc 
names       : mdet_pauliceia_laea_int 
values      : 566, 751  (min, max)

Esse arquivo será utilizado na etapa seguinte, como entrada para um executável de uma rotina Fortran desenvolvida pelo Instituto de Pesquisas Hidráulicas (IPH) como parte da preparação de dados de entrada do modelo hidrológico MGB-IPH (Paz 2008, Paz et al. 2006).

4. Definição do domínio da bacia

A definição do domínio espacial da bacia hidrográfica é etapa crucial da preparação de dados, pois o número de colunas e linhas (e consequentemente as coordenadas dos vértices do domínio espacial) devem ser definidas levando-se em conta a resolução na qual será aplicado o DBHM. Em geral essa resolução é na ordem de 1 km. então para que seja possível aplicar uma metodologia de upscaling as linhas e colunas devem ser múltiplas de 10. Assim um upscaling de 10 vezes como no caso de 100 m para 1000 m o domínio espacial será consistente nas duas resoluções. Então nós vamos truncar o número de colunas e linhas do MDET (mdet_laea_i) para números múltiplos de 10. Essa mudança deve ser feita nas extent do raster, recortando-o para uma nova extent, tal que o total de linhas e colunas resultantes seja múltiplo de 10. A função adjust_extent (definida no script ../scripts/adjust_extent.R) foi desenvolvida para isso. O script está comentado, consulte-o para detalhes do procedimento.

## ajuste da extent do raster par linhas e colunas multiplas de 10
mdet_laea_crop_file <- adjust_extent(ifile = mdet_laea_idrisi_file)

------------------------------------------------- 
../data_base/pauliceia/mdet_pauliceia_laea_int.rst 
class       : 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. : NA 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/data_base/pauliceia/mdet_pauliceia_laea_int.rst 
names       : mdet_pauliceia_laea_int 
values      : 566, 751  (min, max)

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 : in memory
names       : mdet_pauliceia_laea_int 
values      : 572, 751  (min, max)

mdet_laea_crop_file

[1] "../data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst"

## importando arquivo raster recortado 
mdet_laea_crop <- raster(mdet_laea_crop_file)
mdet_laea_crop

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/data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst 
names       : mdet_pauliceia_laea_int_crop 
values      : 572, 751  (min, max)

O raster recortado foi salvo com o nome ../data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst. Este arquivo raster é o principal para realização da etapa seguinte.

5. Lista de arquivos gerados

• mdet_basin_laea_int.tif MDET da bacia na projeção cartográfica do DBHM (LAEA);

mdet_laea_i

class       : 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)

mdet_laea_filename_i

[1] "../data_base/pauliceia/mdet_pauliceia_laea_int.tif"

• mdet_basin_laea_int_crop.rst MDET do domínio definitivo da bacia na projeção cartográfica do DBHM (LAEA) e no formato do software Idrisi;

mdet_laea_crop_file

[1] "../data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst"

mdet_laea_crop

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/data_base/pauliceia/mdet_pauliceia_laea_int_crop.rst 
names       : mdet_pauliceia_laea_int_crop 
values      : 572, 751  (min, max)

7. Referências citadas

Paz, A.R. MGBgis - Modelo hidrológico distribuído MGB-IPH, ETAPA 1 de preparação das informações de entrada. Manual do Usuário, IPH-UFRGS, Porto Alegre (RS). 2008.

Paz, A. R., W. Collischonn, e A. L. L. Silveira (2006), Improvements in large scale drainage networks derived from digital elevation models, Water Resourc. Res., 42 (8), doi: 10.1029/2005WR004544.