domingo, dezembro 30, 2012

Acessando Dados do PostGIS no R

Banco de Dados Geográficos aliados a softwares estatísticos possibilitam as mais complexas análises, isso pode ser feito através da biblioteca rgdal no R. Dessa forma é possível fazer a leitura das tabelas espaciais do PostGIS neste ambiente, veja o exemplo:

1 - Para abrir uma sessão do R, abra o terminal e digite R (Figura 1):

Figura 1 - Sessão do R no Terminal

2 - Dentro da sessão, instale os pacotes abaixo: 
install.packages("sp")
install.packages("gstat")
install.packages("rgdal")

3 - Agora execute os comandos abaixo:

# carregando a biblioteca rgdal
library(rgdal)
# criando o objeto streams, a partir da leitura da tabela streams
# no banco de dados ghydroweb
streams <-readOGR("PG:dbname=ghydroweb", "streams")
# vendo os dados, através da plotagem do gráfico 
plot(streams, axes="true", col="blue")
# adicionando uma grade ao gráfico
grid()

4 - O resultado pode ser conferido na Figura 2:

Figura 2 - Gráfico dos dados da tabela streams


Fiquem de olho nas próximas publicações sobre o R! Um abraço e até o ano que vem :)




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



sábado, novembro 24, 2012

Geocodificação e Geocodificação Reversa no QGIS


A Geocodificação é o processo de conversão de endereços (ex: Rua Augusta, 2862 - Jardim América, São Paulo) em coordenadas geográficas (X,Y). Já a Geocodificação Reversa permite que você encontre um endereço a partir de uma coordenada informada.

Através do plugin GeoCoding no Quantum GIS é possível realizar os dois processos, da seguinte forma:

1- Abra o instalador de plugins do QGIS:




2- Digite a palavra geocoding e em seguida clique em Install plugin:



3- Se a instalação for bem sucedida, clique em OK na próxima janela:

 

4- É necessário definir o sistema de referência antes de utilizar o plugin, caso não tenha feito, digite Control+Shift+P para abrir a próxima janela. Neste exemplo, foi utilizado o EPSG:29195:




5- No menu Plugins, escolha a opção Geocode -> Geocode:



6- Digite o endereço que você quer geocodificar (ex: Av. Rio Grande do Sul, 40, Bairro dos Estados, João Pessoa):

 7- Se o endereço for encontrado, uma nova camada será criada (GeoCoding Plugin Results), contendo o ponto geocodificado:



8- Para realizar a geocodificação reversa, escolha a opção Geocode -> Reverse GeoCode no menu Plugins:




9 - Clique no local em que se deseja obter o endereço:




10- Se o endereço for encontrado, surgirá uma janela com o resultado, clique em OK:



11- Será criada uma feição pontual, contendo o endereço encontrado:




Espero que seja útil, até a próxima! o/

sábado, novembro 10, 2012

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 

terça-feira, fevereiro 28, 2012

Curso online de Mapserver 6



Olá Amigos,

Estarei ministrando um treinamento online sobre Mapserver 6 a partir do dia 24/03/2012 através da Acesso Livre Consultoria.

Maiores detalhes sobre o curso neste link.

Garanta já a sua vaga realizando a pré-incrição através deste formulário: http://goo.gl/LqrEm

Um abraço e até lá :)

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