Terça-feira, Fevereiro 28, 2012

Curso on line de Mapserver 6



Olá Amigos,

Estarei ministrando um treinamento on line 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 :)