DBHM-GIS: PARTE 0

0.1 Download do Modelo digital de elevação do terreno para a bacia de interesse
0.2. Pré-requisitos
- Software
- Arquivos de dados hidrográficos em formato SIG
0.3. Região em torno da bacia
0.4. Download do modelo digital de elevação do terreno (MDET)
0.6 Estrutura de diretórios e dados gerados
0.7. Lista de arquivos gerados

Universidade Federal de Santa Maria - UFSM

Departamento de Física

Prof. Jônatan Tatsch

Aluno Roilan Hernandez Valdes

0.1 Download do Modelo digital de elevação do terreno para a bacia de interesse

O modelo Hidrológico Distribuído da Biosfera (DBHM) requer diversas informações de entrada para sua aplicação ao domínio espacial da bacia hidrográfica de interesse. O DBHM utiliza dados de sensoriamento remoto, dados de estações meteorológicas e mapas de solo e vegetação. A manipulação de arquivos em formato SIG que armazenam dados do tipo vetorial e raster é realizada por sistemas de informação geográfica (SIG). Esse documento relata como fazer o download do modelo digital de elevação do terreno (MDET) para bacias hidrográficas. A bacia do Rio Mogi-Guaçú (Rio Vacacaí) localizada entre o nordeste do estado de SP e Sul de MG (Centro do RS) é usada como exemplo de aplicação.

0.2. Pré-requisitos

Para organização dessa etapa de processamento de dados recomenda-se criar um diretório com a estrutura my_path/DBHM/data_process/geo_data. Em que my_path deve ser um diretório que você tenha acesso irrestrito, como por exemplo /home/user/. O diretório obtido da descompactação do arquivo simCTL.tar.gz (com a simulação teste para a bacia do Rio Mogi-Guaçú com o compilafor Intel) pode ser armazenado no diretório my_path/DBHM/. Essa estrutura de diretórios (Figura 1.2) será usada como referência daqui para frente. Portanto recomenda-se mover todas simulação já realizadas para este diretório.

Software

Para manipulação e visualização de dados serão usados os softwares:

R

Para instalação do R consulte o arquivo html my_path/DBHM/data_process/geo_data/docs/R_install.html. Certifique-se de que os pacotes indicados no final desse arquivo foram instalados. Instale também o pacote rgl para visualização 3D da elevação do terreno.
```
$ sudo apt-get install r-cran-rgl
```

Google-earth

Para instalá-lo, no linux Ubuntu digite no terminal (ctrl+Alt+t):

$ sudo apt-get install googleearth-package
$ make-googleearth-package 
$ sudo dpkg -i googleearth_6.0.3.2197+1.1.0-1_amd64.deb
$ sudo apt-get -f install

Arquivos de dados hidrográficos em formato SIG

O diretório my_path/DBHM/data_process/geo_data/data_base/ANA contém os dados hidrográficos das regiões hidrográficas brasileiras do Rio Paraná (código 6) e Atlântico-Trecho Sudeste (código 8). Esses dados foram obtidos no site hidroweb) da Agência Nacional das Águas (ANA). A bacia do Rio Mogi-Guaçú localiza-se na região hidrográfica do Paraná (8).

Os arquivos mais relevantes são:

shapefile das regiões hidrográficas brasileiras
shapefile de ottobacias da ANA
shapefile de estações fluviométricas
shapefile de estações pluviométricas

Esses arquivos servirão como informações auxiliares na localização aproximada da domínio espacial da bacia. Esse domínio refere-se a um retângulo cujos os vértices envolvem a bacia hidrográfica de interesse. As coordenadas (longitude e latitude) dos vértices do retângulo definem o domínio espacial da bacia hidrográfica, que nesse momento devem ser definidas com uma boa folga em torno da bacia. Uma distância de afastamento em torno de 1/2 grau (~50 km) dos pontos extremos (do contorno) da bacia é suficiente. O contorno que delimita a bacia será inicialmente aproximado pelas informações da ANA (as Ottobacias).

0.3. Região em torno da bacia

Abra o Rstudio, crie um novo script denominado dbhm_gis_parte0.R e salve-o em my_path/DBHM/data_process/geo_data/scripts. Copie os comandos abaixo e cole nesse script a medida que for executando-os no RStudio.

Vamos começar importando os dados hidrográficos em formato shapefile da ANA serão importados no R. O arquivo com a delimitação dos estados brasileiros foi criado anteriormente na etapa de instalação do R e do pacote raster (veja o arquivo my_path/DBHM/data_process/geo_data/docs/R_install.html para gerá-lo).

Para importar os dados no R são utilizados alguns pacotes que devem ser instalados. Use o comando install.packages(pkgs = "nome_pacote", dependencies = TRUE) para instalar um pacote.

Caso tenha dúvidas ou erros na execução de alguma função, consulte o help da função utilizada no comando através do comando, por exemplo, ?gsub.

Configurações básicas no R

Vamos começar carregando os pacotes necessários.

## carregando pacotes necessários
require(maptools)
require(raster)
require(rgdal)
require(sp)
require(ggplot2)
require(rasterVis)
require(rgdal)

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 que podem 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. Nessa parte nós só usaremos o script getData.R.

## configurando opção para mostrar barra de progresso durante operações com raster
rasterOptions(progress = "text",tmpdir = "/tmp")
## defina seu nome de usuario que definira o caminho até o diretório scripts
user_name <- system("echo $USER", intern =  TRUE )
#user_name <- "pqgfapergs1"
## 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"

## carregando scripts necessários
## carregando função getData corrigida
source("getData2.R")

Definição de diretórios e nome de arquivos

Vamos definir a localização do diretório com a base de dados (home/user/DBHM/data_process/geo_data/data_base) e definir o nome da bacia de interesse.

## local para salvar dados
data_base_path <- "../data_base"
## nome da bacia
#basin_name <- "mogu"
#basin_name <- "vacacai"
basin_name <- "pauliceia"
data_base_path_basin <- paste0(data_base_path, "/", basin_name)
data_base_path_basin

[1] "../data_base/pauliceia"

Em seguida definimos o código da região hidrográfica em que está inserida a bacia de interesse: 8 para Atlântico leste e 6 para Paraná. Atere a variável basin_id conforme a sua bacia. Essa variável será usada para definir o caminho do shapefile da ANA da região hdrográfica.

basin_id <- 8
if(basin_name %in% c("mogu", "pauliceia")) basin_id <- 6
## caminho para os arquivos da bacia de interesse
basin_path <- gsub(pattern = "X",
                   replacement = basin_id,
                   x = "../data_base/ANA/BaciaX")
basin_path

[1] "../data_base/ANA/Bacia6"

## lista com nome dos arquivos shapefiles da baciaX
shapefiles_list <- list.files(basin_path,
                              pattern = "shp$", 
                              recursive = T, 
                              full.names = T)
shapefiles_list

 [1] "../data_base/ANA/Bacia6/Bacia 6.shp"                           
 [2] "../data_base/ANA/Bacia6/Est_Fluviometricas_ANA_2010.shp"       
 [3] "../data_base/ANA/Bacia6/Est_Fluviometricas_outras_ent_2010.shp"
 [4] "../data_base/ANA/Bacia6/Est_Pluviometricas_ANA_2010.shp"       
 [5] "../data_base/ANA/Bacia6/Est_Pluviometricas_outras_ent_2010.shp"
 [6] "../data_base/ANA/Bacia6/Est_Qualidade_ANA_2010.shp"            
 [7] "../data_base/ANA/Bacia6/Est_Qualidade_outras_ent_2010.shp"     
 [8] "../data_base/ANA/Bacia6/Est_Sedimentos_ANA_2010.shp"           
 [9] "../data_base/ANA/Bacia6/Est_Sedimentos_outras_ent_2010.shp"    
[10] "../data_base/ANA/Bacia6/Est_Telemetricas_ANA_2010.shp"         
[11] "../data_base/ANA/Bacia6/Est_Telemetricas_outras_ent_2010.shp"  
[12] "../data_base/ANA/Bacia6/Hidrografia 1000000.shp"               
[13] "../data_base/ANA/Bacia6/Hidrografia 250000.shp"                
[14] "../data_base/ANA/Bacia6/Municipios.shp"                        
[15] "../data_base/ANA/Bacia6/Rodovias.shp"                          
[16] "../data_base/ANA/Bacia6/Sedes munipais.shp"                    
[17] "../data_base/ANA/Bacia6/Sub_bacias.shp"

## nome com caminho completo para os arquivos
basin_file <- grep(x = shapefiles_list, 
                   pattern = paste("Bacia", basin_id), 
                   value = T)
basin_file

[1] "../data_base/ANA/Bacia6/Bacia 6.shp"

Agora vamos carregar o shapefile com localização dos postos fluviométricos da ANA, definir o nome do arquivo com a hidrografia da ANA e dos arquivos de Ottobacias da ana para os níveis 1,2 e 4.

## nome do shapefile das estações fluviométricas
flu_file    <- grep(x = shapefiles_list, 
                     pattern = "Est_Fluviometricas_ANA_2010", 
                     value = T)
flu_file

[1] "../data_base/ANA/Bacia6/Est_Fluviometricas_ANA_2010.shp"

## nome do shapefile com a rede drenagem da bacia de interesse
dren_file   <- grep(x = shapefiles_list, 
                     pattern = "Hidrografia 1000000", 
                     value = T)
dren_file

[1] "../data_base/ANA/Bacia6/Hidrografia 1000000.shp"

## nome dos arquivos de ottobacias
otto1_file  <- "../data_base/ANA/Ottobacias_Nivel1/Ottobacias_Nivel1.shp"
otto2_file  <- "../data_base/ANA/Ottobacias_Nivel2/Ottobacias_Nivel2.shp"
otto4_file  <- "../data_base/ANA/Ottobacias_Nivel4/Ottobacias_Nivel4.shp"
## lista de todos onbjetos criados nessa sessão do R
ls()

 [1] "basin_file"           "basin_id"             "basin_name"          
 [4] "basin_path"           "data_base_path"       "data_base_path_basin"
 [7] "dren_file"            "flu_file"             "getData2"            
[10] "my_path"              "otto1_file"           "otto2_file"          
[13] "otto4_file"           "script_dir"           "shapefiles_list"     
[16] "user_name"

## lista de objetos (ou variáveis) criadas na sessão do R com a string "_file" no nome
obj_list <- ls(pattern = "_file")
obj_list

[1] "basin_file" "dren_file"  "flu_file"   "otto1_file" "otto2_file"
[6] "otto4_file"

## usando lopping para verificar se cada arquivo existe no caminho definido
file.exists(sapply(obj_list, get))

[1] TRUE TRUE TRUE TRUE TRUE TRUE

Leitura dos shapefiles e conversão para formato KML

Com os caminhos dos arquivos definidos vamos agora importar os diferentes dados, definir o sistema de projeção de cada um e gerar arquivos KML das Ottobacias e Hidrografia da ANA para visualização no Google Earth. Para salvar os dados em formato KML eles devem ser projetados para a projeção geográfica (lonlat).

## nome do arquivo com delimitação dos estados brasileiros
br_states_file <- "../data_base/br_states.rds"
## importa delimitação dos estados brasileiros
br_states <- readRDS(file = br_states_file)
## estacoes flu (vazao)
flu <- shapefile(flu_file, encoding="latin1")
  ## projetando flu para lonlat para salvar como kml
  flu <- spTransform(x = flu, CRSobj = CRS(proj4string(br_states)))
## bacia brasileira
basin <- shapefile(basin_file, encoding="latin1")
  ## projetando basin para lonlat para salvar como kml
  basin <- spTransform(x = basin, CRSobj = CRS(proj4string(br_states)))
## rede de drenagem
dren <- shapefile(dren_file,  encoding="latin1")
 ## selecionando somente as variáveis de interesse no slot de dados 
 ## do objeto dren (SpatialLinesDataFrame)  
 dren@data <- subset(dren@data, sel = c("NORIO", "NORIOCOMP", "COBACIA", "CORIO"))                                     
  ## projetando dren para lonlat para salvar como kml
  dren <- spTransform(x = dren, CRSobj = CRS(proj4string(br_states)))
  ## gerando arquivo KML para visualização no google-earth
  ## salva no dir data_base_path
  KML(x = dren, file = paste0(data_base_path, "/","dren.kml"), overwrite = TRUE)
## ottobacias
otto1 <- shapefile(otto1_file, encoding="latin1")
 otto2 <- shapefile(otto2_file, encoding="latin1")
  otto4 <- shapefile(otto4_file, encoding="latin1")
   otto4 <- spTransform(x = otto4, CRSobj = CRS(proj4string(br_states)))
   ## gerando arquivo KML para visualização no google-earth
   ## salva no dir data_base_path
   KML(x = otto4, file = paste0(data_base_path, "/","otto4.kml"), overwrite = TRUE)

          used (Mb) gc trigger  (Mb) max used  (Mb)
Ncells 1589796 85.0    3205452 171.2  3205452 171.2
Vcells 4367523 33.4    9127424  69.7  9084285  69.4

Visualização dos dados hidrográficos até a área de interesse

Para chegarmos a nossa região de interesse vamos gerar gráficos com a delimitação dos estados brasileiros e sobrepor os polígonos de Ottobacias da ANA.

## estados
#x11()
plot(br_states, border = "white", col = "gray", lwd = 2, axes = T)
## ottobacias nível de divisão 4
plot(otto4, add = T, border = "blue")
## ottobacias nível de divisão 2
plot(otto2, add = T, border = "red", lwd = 3)
## ottobacias nível de divisão 1
plot(otto1, add = T, border = "black", lwd= 4)

A função zoom pode ser usada para aproximar-se da região da bacia de interesse. Ao executá-la, aguarde até que o cursor do mouse forme uma cruz. Desenhe um retângulo, começando com um clique no vértice do canto esquerdo superior do retângulo que cobrirá a região alvo e depois outro clique no vértice do canto direito inferior. Após clicar nos dois pontos surgirá um tela gráfica em branco. Rode novamente os comandos para sobrepor os estados e bacias para visualização. Repita esse processo o quanto for necessário para chegar ao domínio da sua bacia de interesse.

## zoom na regia de interesse
zoom(otto4)
  ## estados
  # plot(br_states, add = T, border = "white", col = "gray", lwd = 2)
   ##  grande bacia
   # plot(basin, add = T, axes = T, lwd = 2, lty = 2)
    ## ottobacias nível de divisão 4
    plot(otto4, add = T, border = "blue")
     ## drenagem
     plot(dren, add = T, col = 5)
      #3 postos fluviométricos
      plot(flu, add = T, pch = 20,  cex = 1.2, col = 2)
    ## repita os comandos acima até chegar no domínio desejado
    zoom(otto4)
     plot(otto4, add = T, border = "blue")
      ## drenagem
      plot(dren, add = T, col = 5)
       ## postos fluviométricos
       plot(flu, add = T, pch = 20,  cex = 1.2, col = 2)

0.4. Download do modelo digital de elevação do terreno (MDET)

Quando estiver com a região de interesse na tela gráfica use os comandos abaixopara definir o domínio espacial e as coordenadas centrais desse domínio que servirão para fazer o download do MDET.

if (basin_name == "pauliceia") {
  #*Rhv>> Definindo as coordenadas centrais do dominio da bacia.  
lon_med <- mean(x = c(-47.64479, -47.60735) )
lat_med <- mean(x = c(-21.6588, -21.61294) )
  
} else {

## definindo domínio aproximado
 e <- drawExtent()
 e
  ## criando raster a partir da extent
  r <- raster(ext = e)
  r
   ## preenchendo raster com numero sequencial de celulas
   r <- setValues(r, 1:ncell(r))
   plot(r)
    ## obtendo as coordendas de cada célula do raster r
    coords <- xyFromCell(r, cell = 1:ncell(r))
    ## long central da região
    lon_med <- mean(coords[, "x"])
    ## lat central da região
    lat_med <- mean(coords[, "y"])
}
lon_med

[1] -47.62607

lat_med

[1] -21.63587

O download do MDET requer os valores de um ponto no domínio, nesse caso usaremos as coordendas centrais do domínio.

## CHUNK OCULTO
tile_file <- list.files(data_base_path_basin, pattern = "^srtm.*tif$|mosaico", full.names = TRUE)
if(basin_name == "vacacai") {
  tile_file <- grep("mosaico",tile_file, value = T)
} else {
  if(basin_name == "mogu" ) tile_file <- grep("^srtm.*tif$",tile_file, value = T)
  if(basin_name == "pauliceia") {
      tile_file <- grep("^srtm.*tif$",tile_file, value = T)
#      tile_file <- "asterV2TilePauliceia.img"
#      tile_mdet <- raster(paste0(data_base_path_basin, "/",tile_file))
  }
}
if(file.exists(tile_file)) tile_mdet <- raster(tile_file)
tile_file

## download do modelo digital de elevação do terreno (SRTM)

cat("Fazendo download do Tile do SRTM ... pode ir tomar um café", "\n")

Fazendo download do Tile do SRTM ... pode ir tomar um café

tile_mdet <- getData2(name = "SRTM", 
                     lon = lon_med, 
                     lat = lat_med, 
                     download = T,
                     path = data_base_path)

../data_base/srtm_27_17.ZIP 
../data_base/srtm_27_17.tif

Vamos visualizar tile do MDET que cobre uma ampla região que provavelmente inclui a bacia. Os dados hidrográficos da ANA serão sobrepostos ao tile do MDET.

## gráfico para visualização do mdet
plot(tile_mdet)
 plot(otto4, add = T, border = 4, lwd = 2)
   plot(dren, add = T, col = 5)
    plot(flu, add = T, pch = 20,  cex = 1.2, col = 2)
e <- extent(-47.64479-0.015, -47.60735+0.01, -21.6588-0.015, -21.61294+0.01)
e

class       : Extent 
xmin        : -47.65979 
xmax        : -47.59735 
ymin        : -21.6738 
ymax        : -21.60294

plot(e, add = TRUE)

Vamos repetir o precedimento usado para determinar as coordenadas centrais da região da bacia, ou seja dar zoom até definirmos um domínio espacial entorno da bacia da interesse. Novamente esse domínio deve ser definido com uma boa folga em torno da bacia. Uma distância de afastamento em torno de 1/2 grau (~50 km) dos pontos extremos (do contorno) da bacia é suficiente. Repita o procedimento se necessário. Em seguida vamos recortar o tile do MDET de acordo com a extent definida.

## definindo domínio aproximado
e <- extent(-47.64479-0.015, -47.60735+0.01, -21.6588-0.015, -21.61294+0.01)
e
plot(e, add = TRUE)

## mdet recortado para o domínio definido anteriormente
mdet_c <- crop(tile_mdet, e)
dren_c <- crop(dren, mdet_c)

## gráfico para visualização do mdet_c
plot(mdet_c)
 #plot(otto4, add = T, border = 4, lwd = 2)
   plot(dren_c, add = T, col = 5)
    plot(flu, add = T, pch = 20,  cex = 1.2, col = 2)

Por fim salvamos o raster recortado para região em torno da bacia em formato tif e ascii.

## salvando mdet da bacia em formato ARC/INFO ASCII GRID e TIF
mdet_ascii_file <- paste0(data_base_path_basin, "/mdet_", basin_name,".asc")
  mdet_tif_file <- paste0(data_base_path_basin, "/mdet_", basin_name,".tif")

writeRaster(mdet_c, filename = mdet_ascii_file, overwrite = T)

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. : NA 
data source : /home/hidro2roilan/DBHM/data_process/geo_data/data_base/pauliceia/mdet_pauliceia.asc 
names       : mdet_pauliceia

writeRaster(mdet_c, filename = mdet_tif_file, overwrite = T)

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)

O arquivo salvo no formato ARCINFO ASCII GRID pode ser aberto em um editor de texto para visualização e/ou edição de seus dados. Ambos arquivos podem ser importados no R pela função raster (r<- raster('file.tif'); r<- raster('file.asc')) do pacote raster.

O Plot 3D do mdet_c permite a interação com o gráfico. Navegue pela bacia para conhecer a topografia da região.

plot3D(mdet_c)

0.6 Estrutura de diretórios e dados gerados

A figura abaixo mostra a estrutura de diretórios com destaque para os arquivos gerados nessa etapa que são mais importantes .

0.7. Lista de arquivos gerados

mdet_basin.tif MDET da região em torno da bacia na projeção geográfica (lonlat);

mdet_c

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 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : srtm_27_17 
values      : 566, 752  (min, max)

mdet_tif_file

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