GDAL 1.6, FME e os datums brasileiros

Retornando ao problema de configuração de datums brasileiros, estou postando os dois arquivos que considero fundamentais: o de configurar a gdal e o de configurar o fme.
Espero que seja útil.
Arquivo do FME , Arquivo da GDAL

GDAL e Sistemas de Coordenadas

Como vimos anteriormente, uma das grandes aplicações da GDAL é na conversão de coordenadas. Esta semana me peguei querendo uma transformação rápida de coordenadas e utilizei uma sequência de comandos Python no terminal interativo para resolver esse problema.

Uma das vantagens do Python é que no terminal interativo pode-se acompanhar os resultados conforme eles aparecem.




Neste exemplo cada sistema de coordenadas é definido por um objeto da classe "osr.SpatialReference" (Referência Espacial). Para definir os sistemas de coordenadas novamente utilizamos a referência EPSG de cada um destes no método "osr.SpatialReference.ImportFromEPSG()".

Por fim, é criado um objeto "osr.CoordinateTransformation" (Transformação de Coordenadas) partindo de srE para srS. Este objeto vai ser o responsável por transformar todos pares (ou triplas) de coordenadas que eu passar para ele através do método "osr.CoordinateTransformation.TransformPoint(x,y,z)" Repetindo esse comando várias vezes pode-se transformar várias coordenadas. É uma calculadora rápida e precisa de coordenadas. E serve como exemplo de como definir sistemas de referência na GDAL em Python.

Espero que seja útil.
E até o próximo post.

Reprojetando todo o servidor Raster

function teste {
}
import sys, os
try:
import gdal #gdal <1.6
except:
from osgeo import gdal #gdal >1.6
import os, sys
gdalwarp=r"gdalwarp.exe" #Para funcionar, coloque o gdalwarp dentro do path usando o script set_fw.bat
lista = []

#-->procurando os fotolitos dentro da estrutura de diretorios - inicio
for root, dirs, file in os.walk(os.getcwd()):
for i in file:
if (i[-4:]==".cit" ) or (i[-4:]==".tif"):
if ("Georreferenciados" in root):
lista.append(os.path.join(os.getcwd(),root,i))
#-->procurando os fotolitos dentro da estrutura de diretorios - fim

def corrigeGeo(arquivo): #transforma de centimetros para metro
ds=gdal.Open(arquivo)
georef=ds.GetGeoTransform()
a=[]
for i in range(0,len(georef)):
a.append(georef[i]*0.01)
georef=ds.SetGeoTransform(a)
del ds

print "Iniciando conversoes"
for i in lista:
print i
corrigeGeo(i) #transforma de centimetros para metro
#converte os arquivos listados do fuso 22 Corrego Alegre para fuso 22 WGS e salva em imagens de 1 bit
os.system(gdalwarp+" "+i+" "+i.split(".")[0]+'_WGS.tif'+' -s_srs EPSG:22522 -t_srs EPSG:32722 -co "NBITS=1"')
print "Terminado"

GDAL, Corrego Alegre e definições de DATUM

No Brasil vivemos uma transição de datum geodésico. Devemos converter os dados de SAD 69 e Córrego Alegre (em diversas materializações) para o SIRGAS 2000. Aqui no Rio Grande do Sul todo o mapeamento é em Córrego Alegre na escala 1:50.000. Então surge um problema: Utilizando as resoluções Resolução 1/2005 e Resolução 23/1989, a transformação de Córrego Alegre para SIRGAS 2000 pode ser feita somando as coordenadas de transformação geocêntrica utilizando a transformação de Molodensky com os parâmetros X=-206,048m Y=168,279m Z=-3,823m.

Porém, no site do IBGE foi dada uma nova forma de transformar. O programa PROGRID utiliza grades de valores e efetua transformações que difereriram em alguns metros das coordenadas calculadas com os parâmetros previstos nas resoluções. Perguntei a um professor meu da graduação e ele me orientou que o PROGRID foi uma solução que buscava obter melhor coerência entre as transformações de coordenadas.

Apesar disso, nenhum software de GIS atualmente adotou essa transformação e não achei nenhuma referência no IBGE falando sobre alguma resolução que oficialize esse método de transformação. Como estou trabalhando com cartas 1:100.000, os 5m de erro absoluto não são relevantes durante a conversão de dados analógicos para digitais.

Depois de contar a história triste que motivou esse estudo, vamos pra parte prática: Tínhamos em torno de 50 arquivos TIF para converter para esta semana e mais uns 500 no servidor aguardando oportunidade. Decidimos automatizar a conversão.

Como prefiro, optamos por utilizar tecnologia aberta, escolhendo a GDAL. Para quem não conhece, esta é a biblioteca principal de programas C/C++ e python para dados geográficos matriciais e vetoriais. Ela vem com diversos programas auxiliares, entre eles, um programa de reprojeção de arquivos matriciais ("gdalwarp"). Utilizei o pacote pré-compilado FWTools em Windows XP durante os testes.

Para definir os parâmetros corretos, deve-se alterar arquivos na pasta "data" da instalação da gdal (no meu caso "C:\Program Files\FWTools2.4.2\data"). O arquivo "data\coordinate_operation.csv" guarda o nome das transformações, o identificador delas e os identificadores dos sistemas de coordenadas (referentes à EPSG). Procurando por "rrego" neste arquivo achei a linha:

1132,Corrego Alegre to WGS 84 (1),transformation,4225,4326,...

Que representa exatamente a transformação pré-definida na GDAL.

Dentro do arquivo "data\coordinate_operation_parameter_value.csv" procurei o identificador 1132 e achei as linhas:

1132,9603,8605,-206,,9001
1132,9603,8606,172,,9001
1132,9603,8607,-6,,9001

Que eu troquei para:

1132,9603,8605,-206.048,,9001
1132,9603,8606,168.279,,9001
1132,9603,8607,-3.823,,9001

E desta forma, minha instalação da GDAL passou a seguir os parâmetros previstos nas resoluções.

Agora meu gdalwarp já fazia as transformações corretamente, mas queríamos transformar muitos arquivos. No próximo post comento o código que eu e César Soares escrevemos para fazer isso em toda a estrutura.

Novamente, agradeço a todos que se interessaram pelas informações aqui postadas pois é esse o objetivo.
Até a próxima!