Mostrando postagens com marcador GRASS. Mostrar todas as postagens
Mostrando postagens com marcador GRASS. Mostrar todas as postagens

segunda-feira, junho 16, 2014

Perfil Longitudinal do Rio Principal Utilizando o GRASS e o R

Por definição, o perfil longitudinal de um rio apresenta a relação existente entre a variação altimétrica e o comprimento do mesmo, desde a nascente até a foz ou ponto de confluência.

Neste tutorial daremos continuidade ao post "Caracterização Morfométrica de Bacias Hidrográficas utilizando o GRASS" onde a partir do rio principal, será gerado o gráfico do perfil longitudinal do mesmo, utilizando o GRASS e o R.

Com o GRASS aberto no terminal, digite os comandos a seguir:

Como resultado, teremos o gráfico abaixo:


Referências:

terça-feira, dezembro 03, 2013

Caracterização Morfométrica de Bacias Hidrográficas utilizando o GRASS

A caracterização morfométrica tem como objetivo esclarecer várias questões associadas ao relevo e a dinâmica do escoamento superficial de uma bacia hidrográfica. Tal caracterização é feita com base nas principais medidas físicas e serve de instrumento para diversas aplicações relacionadas aos recursos hídricos.

Utilizando a extensão r.basin do GRASS, é possível extrair os atributos hidrológicos que compõem a caracterização morfométrica, para tanto, é necessário que o usuário informe a coordenada do exutório da bacia hidrográfica. Neste exemplo irei utilizar a imagem SRTM SB-24-Z-C como área de estudo.

1 - Convertendo o sistema de projeção da imagem SRTM para UTM


Utilizando o terminal, entre na pasta que contém a imagem e digite o seguinte comando:
$ gdalwarp -s_srs "EPSG:4326" -t_srs "EPSG:31984" SB-24-Z-C.tif dem.tif

Desta forma o raster será reprojetado para a projeção UTM Fuso 24 Sul (Datum SIRGAS 2000).

2 - Configurando a Location e o Mapset do GRASS


Abra o GRASS na janela do terminal (digitando "grass") e siga os próximos passos:

  • Escolha Location Wizard:

  • Defina o nome da Location e clique em next:

  • Na próxima tela, escolha a opção indicada abaixo:

  • Clique em "Browse" e procure o raster dem.tif:

  • Será exibida a tela de sumário, clique em "Finish"

  • No próximo passo, clique em "Yes":

  • Defina o nome do mapset e clique em "Start GRASS":



3 - Modificando as cores do Raster:


No terminal onde a sessão do GRASS está aberto, digite:
$ r.colors map=dem color=elevation

Para visualizar o resultado no  "Map Display", adicione esta camada através da opção "Add rast map layer" (Control + Shift + r):


4 - Preparando o MDE para a análise morfométrica:


Com o comando r.fill.dir são removidos os "sinks" (células cercadas por outras com maiores elevações), que produzem a descontinuidade do escoamento superficial descendente para uma célula vizinha.
$ r.fill.dir input=dem elevation=dem.filled direction=filled.dir


O próximo passo consiste em criar duas superfícies para análises hidrológicas através do comando r.watershed, são elas: a superfície de fluxo acumulado (accumulation) e a rede de drenagem matricial (streams).
$ r.watershed -fa elevation=dem.filled accumulation=accumulation \
stream=streams convergence=5 threshold=62

OBS: O threshold define o tamanho mínimo a ser considerado como uma área de contribuição para a bacia hidrográfica.
Para áreas de 0,5 km2 (500.000 m2) foi utilizado um valor igual a 62  (500.000 / 90x90), onde 90x90 é a resolução espacial de uma imagem SRTM).

5 - Instalando as extensões necessárias:


Instale uma a uma, as extensões a seguir:
$ g.extension extension=r.stream.angle
$ g.extension extension=r.stream.basins
$ g.extension extension=r.stream.del
$ g.extension extension=r.stream.distance
$ g.extension extension=r.stream.extract
$ g.extension extension=r.stream.order
$ g.extension extension=r.stream.pos
$ g.extension extension=r.stream.preview
$ g.extension extension=r.stream.stats
$ g.extension extension=r.ipso
$ g.extension extension=r.wf

E por último, a extensão r.basin:
$ g.extension extension=r.basin


6 - Realizando a caracterização morfométrica:


Deve-se escolher como coordenada do exutório um ponto que toque a rede de drenagem matricial, neste caso a coordenada escolhida foi a seguinte: (619.241,988946 E || 9.207.668,119428 N)
$ r.basin map=dem.filled prefix=out \ 
easting=619241.988946 northing=9207668.119428 threshold=62

7- Resultados:


O comando r.basin gera os seguintes resultados:

  • Layers: limite da bacia hidrográfica, rio principal, exutório e rede de drenagem da bacia:
  • Gráfico: Curva hipsométrica:

  • Planilha: é gerado um arquivo csv contendo o resultado da análise:


  • Saída no terminal:


sexta-feira, setembro 27, 2013

Krigagem Ordinária no R

O código abaixo exemplifica a utilização da Krigagem Ordinária utilizando o R dentro de uma sessão do GRASS. O intuito desta postagem é mostrar a codificação básica necessária para se chegar a um primeiro resultado e não será abordado aqui como realizar o ajuste do variograma.



Foram utilizados 300 pontos cotados para gerar a superfície mostrada na Figura 1.

Figura 1 - Resultado da Krigagem

quinta-feira, maio 30, 2013

Script para criação de mapas de declividade no GRASS

O script a seguir cria um mapa de declividade a partir de um Modelo Digital de Terreno, de acordo com classificação proposta por Lepsch (1991), descrita na Tabela 1.

Classe de Declividade  Porcentagem
Plano 0 - 2
Suave Ondulado  2 - 5
Moderadamente Ondulado 5 - 10
Ondulado 10 - 15
Forte Ondulado 15 - 45
Montanhoso 45 - 70
Escarpado > 70
Tabela 1 - Classes de declividade (LEPSH, 1991)

Para utilizar o script é necessário fazer o seu download (salve-o na pasta home, por exemplo) e lhe dar permissões de execução através do comando abaixo:
$ sudo chmod +x declividade.sh

Em seguida execute o comando dentro de uma sessão do GRASS:



Resultado:


REFERÊNCIAS:


LEPSCH, I. F. Manual para levantamento utilitário do meio físico e classificação de terras no sistema de capacidade de uso. Campinas : Sociedade Brasileira de Ciência do Solo, 1991.

segunda-feira, fevereiro 11, 2013

Como importar vários arquivos para o GRASS via Terminal

Este laço deve ser aplicado dento do diretório que contém os dados que serão importados para o GRASS:
for i in *.<extensao>; do r.in.gdal -o in=$i out=`basename $i`; done
Onde <extensao> é a extensão dos arquivos raster do diretório atual.

quarta-feira, dezembro 19, 2012

GRASS GIS: Cálculo do comprimento do rio principal de uma Bacia Hidrográfica

O comprimento do rio principal de uma bacia é um parâmetro frequentemente utilizado em análises hidrológicas, neste tutorial demonstrarei como realizar a extração desse parâmetro utilizando o software GRASS GIS e imagens de radar SRTM. 

Nossa área de estudo compreenderá a cena da imagem SB-24-Z-C (Figura 1).

Figura 1 - Modelo Digital de Elevação da área de estudo no GRASS

#(1) Preenchimento das falhas do MDE:

r.fill.dir input=dem elevation=dem.fill direction=dem.dir     

#(2) Criação das superfícies para análises hidrológicas:

r.watershed -fma elevation=dem.fill accumulation=accum drainage=drainage convergence=5 threshold=123

#(3) Extração da rede de drenagem (ainda no formato raster)

r.stream.extract elevation=dem.fill threshold=123 d8cut=infinity mexp=0 stream_rast=streams direction=drain_dir stream_vect=streams

#(4) Delimitação da Bacia Hidrográfica para o exutório de coordenada
# x=617345.784270 y=9205935.470789 (Figura 2):

r.water.outlet drainage=drainage basin=b1 easting=617345.784270 northing=9205935.470789

#(5) Utilizando álgebra de mapas com o comando r.mapcalc,
# são recortadas as superfícies para a bacia gerada (Figura 3):

r.mapcalc'dem.clip=dem.fill*b1'
r.mapcalc 'accum.clip=accum*b1'
r.mapcalc 'drain.clip=drainage*b1'
r.mapcalc 'streams.clip=streams*b1'
r.mapcalc 'drain_dir.clip=drain_dir*b1'

#(6) Classificação da rede de drenagem de acordo com 
# strahler (Figura 4) e hack (Figura 5):

r.stream.order stream=streams.clip dir=drain_dir.clip strahler=r_strahler hack=r_hack

#(7) Determinação do rio principal

r.mapcalc 'r_mainchannel=if($r_hack==1,1,null())'

#(8) Refinamento do rio principal e conversão para vetor (Figura 6):

r.thin input=r_mainchannel output=r_mainchannel_thin

r.to.vec input=r_mainchannel_thin output=v_mainchannel feature=line

#(9) Adição da coluna length a camada v_mainchannel,
# em seguida, o cálculo do comprimento em km:

v.db.addcol map=v_mainchannel columns='length double precision'

v.to.db map=v_mainchannel type='line,boundary' layer=1 qlayer=1 option='length' units='kilometers' columns='length'

#(10) Visualizando o valor do comprimento (Figura 7):

echo "SELECT length FROM v_mainchannel" | db.select

Figura 2 - Bacia gerada com o comando r.water.outlet

Figura 3 - Superfícies recortadas em função do limite da bacia


Figura 4 - Classificação da rede de drenagem de acordo com Strahler


Figura 5 - Classificação da rede de drenagem de acordo com Hack


Figura 6 - Rio principal em destaque, já no formato vetorial

Figura 7 - Valor do comprimento do rio principal em quilômetros

Sugestões de Leitura:


  1. Manual do GRASS: http://grass.osgeo.org/grass65/manuals/index.html
  2. Add-Ons para rede de drenagem: http://grasswiki.osgeo.org/wiki/R.stream.*
  3. Delimitação de bacias no GRASS: http://grasswiki.osgeo.org/wiki/Creating_watersheds
  4. Importação de imagens SRTM: http://grasswiki.osgeo.org/wiki/SRTM



terça-feira, novembro 06, 2012

Determinação da Curva Hipsométrica de uma Bacia Hidrográfica com o GRASS GIS

Uma curva hipsométrica representa o estudo da variação da elevação dos vários terrenos da bacia hidrográfica com referência ao nível do mar (Villela e Mattos, 1975).

 Esta curva é traçada lançando-se em sistema cartesiano a cota versus o percentual da área de drenagem com cota superior.

Através do comando r.ipso, um add-on do GRASS é possível calcular esta curva, a partir de um modelo digital de elevação de uma bacia (figura 1), como mostra o exemplo a seguir:


Figura 1 - Modelo Digital de Elevação da Bacia


GRASS 6.5.svn (Litoral):~ > r.ipso -a map=srtm.clip \
image=/home/marcello/Desktop/curva_ipsometrica.png
 100%
Tot. cells 72802.0
===========================
Ipsometric | quantiles
===========================
171 | 0.025
160 | 0.05
143 | 0.1
111 | 0.25
80  | 0.5
50  | 0.75
58  | 0.7
22  | 0.9
6   | 0.975

Done!

O comando é muito simples, é necessário informar a flag -a (para o cálculo da curva hipsométrica) e mais dois parâmetros: map, nome do MDE da bacia e image, local onde será gerada a imagem.


O resultado do gráfico pode ser visto na figura 2.



Figura 2 - Curva Hipsométrica calculada com o comando r.ipso 

quarta-feira, fevereiro 15, 2012

GRASS: Snapping points to lines

Trabalhar com dados espaciais de diversas fontes nem sempre é uma tarefa fácil, ainda mais quando espera-se que os mesmos obedeçam regras topológicas para que as análises espaciais sejam consistentes.

Um problema topológico comumente enfrentado pode ser visto na Figura 1, onde determinada análise espacial requer que uma camada de pontos tenha que estar devidamente conectada a uma camada linear. Os pontos vermelhos da direita indicam as feições com problemas de conectividade.

Figura 1: Problema topológico - os pontos devem estar conectados aos arcos.

Este problema pode ser resolvido através de uma edição vetorial na camada de pontos, utilizando o recurso snapping, porém se o usuário tiver de editar centenas (ou até milhares) de pontos, este tipo de edição torna-se inviável.

Uma solução é utilizar os recursos do GRASS aliado a versatilidade dos comandos linux para edição/manipulação de arquivos textuais (cat, sed, tail, awk, head, etc.).

Neste exemplo, temos uma rede de drenagem e pontos de controle, que deveriam estar conectados a esta rede (Figura 2). Aparentemente o problema não ocorre (por conta da escala), porém, utilizando a ferramenta zoom é possível identificar que os pontos não tocam os arcos (Figura 3).

Figura 2 - Rede de drenagem e os pontos de controle não conectados.

Figura 3 - Detalhe do problema de conectividade

Segue a solução com os comandos do GRASS:

#adicionar duas colunas que receberão as coordenadas (to_x,to_y), da linha mais próxima do ponto
v.db.addcolumn pontos_controle col="to_x double precision, to_y double precision"

#calcular a menor distância entre os pontos e as linhas e povoar as colunas com as coordenadas do endpoint
v.distance from=pontos_controle to=rios  upload=to_x,to_y column="to_x,to_y"

A Figura 4 mostra a tabela de atributos da camada "pontos_controle" após os dois comandos acima, o par (x,y) representa a coordenada correta, obtida através do cálculo de uma distância euclidiana entre o ponto e coordenada mais próxima da linha.

Figura 4 - Campos povoados com as coordenadas
O comando v.out.ascii permite exportar uma camada vetorial para um arquivo txt:

v.out.ascii input=pontos_controle@analises layer=1 output=/home/marcello/Desktop/pontos_ok.txt columns=to_x,to_y format=point

A Figura 5 mostra o arquivo gerado, observe que as coordenadas corretas são as que estão após o terceiro pipe ( | ), na quarta e quinta coluna:

Figura 5 - Arquivo gerado com o comando v.out.ascii.

Por último, o comando v.in.ascii importa as coordenadas corretas:

#Obs: x=4 e y=5 indicam as colunas com as coordenadas que devem ser lidas na importação
v.in.ascii input=/home/marcello/Desktop/pontos_ok.txt output=pontos_ok x=4 y=5

O resultado pode ser observado na Figura 6, o ponto verde indica a posição correta do ponto  e o vermelho indica a coordenada original.
Figura 6 - Coordenada corrigida (ponto verde).

Caso fosse necessário, poderíamos retirar as colunas que não são do nosso interesse do arquivo ASCII, através da seguinte instrução:

cat pontos_ok.txt | awk -F "|" '{print $4,$5}' > coords_pontos_ok

Ficando o arquivo ASCII como mostra a Figura 7.

Figura 7 - Arquivo  ASCII apoś a remoção das colunas desnecessárias.

That's all folks! até o próximo tutorial.


terça-feira, janeiro 31, 2012

Reprojeção e Clip de Imagens "On the Fly" com o GRASS e GDAL

Lendo o livro "Open Source GIS: A GRASS GIS Approach" (pág. 49) encontrei um procedimento muito útil para resolver o seguinte problema: recortar uma imagem SRTM e reprojetá-la para os mesmos parâmetros definidos no meu projeto. Dessa forma fica bem mais fácil importar a imagem gerada com o comando r.in.gdal. Segue o passo a passo:

1 - Verificar as informações da imagem a ser reprojetada e recortada:
GRASS 7.0.svn (gramame):~ > gdalinfo Geoprocessamento/raster/srtm/SB-25-Y-C.tif
Driver: GTiff/GeoTIFF
Files: Geoprocessamento/raster/srtm/SB-25-Y-C.tif
Size is 1800, 1200
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-36.000415601243731,-7.000415916030761)
Pixel Size = (0.000833333353512,-0.000833333353512)
Metadata:
  TIFFTAG_SOFTWARE=IMAGINE TIFF Support
Copyright 1991 - 1999 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: etif.c $ $Revision: 1.10.1.9 $ $Date: 2003/11/13 19:25:49EST $
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -36.0004156,  -7.0004159) ( 36d 0' 1.50"W,  7d 0' 1.50"S)
Lower Left  ( -36.0004156,  -8.0004159) ( 36d 0' 1.50"W,  8d 0' 1.50"S)
Upper Right ( -34.5004156,  -7.0004159) ( 34d30' 1.50"W,  7d 0' 1.50"S)
Lower Right ( -34.5004156,  -8.0004159) ( 34d30' 1.50"W,  8d 0' 1.50"S)
Center      ( -35.2504156,  -7.5004159) ( 35d15' 1.50"W,  7d30' 1.50"S)
Band 1 Block=1800x2 Type=UInt16, ColorInterp=Gray

2 - Imprimir a região definida como padrão na location:
GRASS 7.0.svn (gramame):~ > g.region -p
projection: 1 (UTM)
zone:       -25
datum:      sam69
ellipsoid:  sam69
north:      9206640
south:      9180540
west:       259380
east:       300780
nsres:      90
ewres:      90
rows:       290
cols:       460
cells:      133400

3 - Usar o comando "eval", para setar como variáveis, a saída do comando anterior (veja mais sobre esse comando em: http://goo.gl/70sAi):
GRASS 7.0.svn (gramame):~ > eval `g.region -g`


4 - Verificar o resultado do comando anterior com o comando "echo", observe que é exibido o retângulo envolvente da região definida como default:
GRASS 7.0.svn (gramame):~ > echo "$w $s $e $n"
259380 9180540 300780 9206640


5 - Verificar a projeção definida na location:
GRASS 7.0.svn (gramame):~ > g.proj -wf
PROJCS["UTM Zone 25, Southern Hemisphere",GEOGCS["sam69",DATUM["South_American_Datum_1969",SPHEROID["South_American_Datum_1969",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]


6 - Usar o comando "gdalwarp" para recortar e reprojetar a imagem a partir da projeção definida no grass e do retângulo envolvente da region:
GRASS 7.0.svn (gramame):~ > gdalwarp -t_srs "`g.proj -wf`" -te $w $s $e $n \
Geoprocessamento/raster/srtm/SB-25-Y-C.tif  \ 
Geoprocessamento/raster/srtm/srtm_clip_utm25.tif
Creating output file that is 450P x 283L.
Processing input file Geoprocessamento/raster/srtm/SB-25-Y-C.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.


7 - Verificar a imagem reprojetada e cortada:
GRASS 7.0.svn (gramame):~ > gdalinfo Geoprocessamento/raster/srtm/srtm_clip_utm25.tif 
Driver: GTiff/GeoTIFF
Files: Geoprocessamento/raster/srtm/srtm_clip_utm25.tif
Size is 450, 283
Coordinate System is:
PROJCS["UTM Zone 25, Southern Hemisphere",
    GEOGCS["sam69",
        DATUM["South_American_Datum_1969",
            SPHEROID["South_American_Datum_1969",6378160,298.25,
                AUTHORITY["EPSG","7036"]],
            AUTHORITY["EPSG","6291"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-33],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (259380.000000000000000,9206640.000000000000000)
Pixel Size = (92.000000000000000,-92.226148409893995)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  259380.000, 9206640.000) ( 35d10'43.72"W,  7d10'19.95"S)
Lower Left  (  259380.000, 9180540.000) ( 35d10'47.83"W,  7d24'29.27"S)
Upper Right (  300780.000, 9206640.000) ( 34d48'14.67"W,  7d10'25.80"S)
Lower Right (  300780.000, 9180540.000) ( 34d48'18.07"W,  7d24'35.32"S)
Center      (  280080.000, 9193590.000) ( 34d59'31.09"W,  7d17'27.72"S)
Band 1 Block=450x9 Type=UInt16, ColorInterp=Gray


8 - Importar a imagem reprojetada e cortada:
GRASS 7.0.svn (gramame):~ > r.in.gdal input=Geoprocessamento/raster/srtm/srtm_clip_utm25.tif output=mde
Projection of input dataset and current location appear to match
 100%
r.in.gdal complete. Raster map created.

Um abraço e até o próximo post :)