Georreferenciamento de dados vetoriais usando o R e o PostGIS

Quem já trabalhou com o ArcGIS provavelmente já deve ter utilizado a extensão Spatial Adjustement que permite o ajuste de bases vetoriais distintas a partir de pontos de controle (conhecido como georreferenciamento de vetores). Neste tutorial irei apresentar uma solução semelhante utilizando a função ST_Affine do PostGIS e o software estatístico R.

Na Figura 1 temos duas camadas vetoriais, a de cor vermelha teve problemas no seu georreferenciamento e será ajustada a partir da base de cor verde. Para isso, primeiro é necessário  importar os shapefiles para um banco de dados geográfico e depois visualizar as tabelas criadas no QGIS:

Figura 1 - Bases cadastrais levantadas por empresas diferentes

Devem ser criadadas duas camadas vetoriais do tipo ponto no BDG onde foram importados os shapefiles. Observa-se na Figura 2 que a localização dos pontos vermelhos correspondem aos pontos verdes (origem -> destino).

Figura 2 - Pontos de controle criados

Em  seguida, em uma sessão no R via terminal deve-se digitar os seguintes comandos:

#carregando a biblioteca e o driver do PostgreSQL
library(RPostgreSQL)
drv <- dbDriver("PostgreSQL")
#conexão com o BD
con <- dbConnect(drv, dbname="cadastro", 
                 user="marcello", 
                 pass="<minha_senha>")
#query 1 - coordenadas x,y dos pontos de destino
pts_ok <- data.matrix(dbGetQuery(con, 
                "SELECT ST_X(geom) AS x,
                 ST_Y(geom) AS y 
                 FROM pontos_ok"))
#query 2 - coordenadas dos pontos de origem
pts_err <- data.matrix(dbGetQuery(con, 
                "SELECT ST_X(geom) AS x,  
                 ST_Y(geom) AS y 
                 FROM pontos_errados"))
#regressão linear: modele pts_ok 
#como função estatística de pts_err 
parameters <- lm(pts_ok ~ pts_err)
#visualizando a matriz transposta dos coeficientes
t(coef(parameters))

Como resultado, teremos os coeficientes da equação que lineariza os pontos com as coordenadas erradas (pts_err) em função dos pontos de referência (pts_ok).

 (Intercept)  pts_errx  pts_erry
x   1312943.9 1.0138267 -0.143060
y   -163870.4 0.1389882  1.013313

Tais coeficientes são os parâmetros de entrada da função ST_Affine:

Formato de saída no R:
   | xoff a b |
   | yoff d e |

Formato de entrada no PostGIS:
 ST_Affine(geom, a, b, d, e, xoff, yoff)

Por último, são executadas as seguintes instruções no BDG:

--criação de uma cópia da tabela errada:
CREATE TABLE lotes_corrigidos AS
SELECT * FROM lotes_errados;

--em seguida o georreferenciamento dos lotes:
UPDATE lotes_corrigidos
SET geom =  
ST_Affine(geom, 
  1.0138267, -0.143060, 
  0.1389882, 1.013313, 
  1312943.9, -163870.4);

Finalizando, a Figura 3 mostra o resultado da base corrigida, o resultado foi satisfatório considerando a quantidade de pontos utilizada:

Figura 3 - Resultado do georreferenciamento (lotes_corrigidos)

PostGIS: Triggers e Procedures (Parte 2)

Complementando o post anterior, imagine o seguinte cenário: um mapa que mostre de forma esquemática, adutoras que serão projetadas para abastecer alguns municípios. Neste caso o objetivo é ter uma ideia dos comprimentos das adutoras (representados pela distância entre as sedes municipais), quantidade de municípios abastecidos e o nome dos mesmos.

CREATE TABLE adutoras (
 gid serial PRIMARY KEY,
 nome_adutora varchar(50),
 -- campo alimentado pela trigger
 comp_trecho_km numeric(10,3),
 -- campo alimentado pela trigger
 qtde_municipios_abast int,
 -- campo alimentado pela trigger
 municipios_abastecidos varchar(400)
 );

-- adição do campo geom
SELECT AddGeometryColumn (
 'public',
 'adutoras',
 'geom',
 '4291', --Lat/Long - Datum SAD69
 'LINESTRING',
 2
 );

-- Procedure
CREATE OR REPLACE FUNCTION pipeline_geometry_field()
RETURNS trigger AS
$$
BEGIN
SELECT string_agg(municipios.nome,', ')
INTO NEW.municipios_abastecidos
FROM municipios
WHERE ST_Intersects(NEW.geom, municipios.geom);

SELECT COUNT(municipios.*) 
INTO NEW.qtde_municipios_abast
FROM municipios
WHERE ST_Intersects(NEW.geom, municipios.geom);

NEW.comp_trecho_km:=ST_Length(Geography((NEW.geom))/1000);

RETURN NEW;
END;
$$
LANGUAGE 'plpgsql';

-- Trigger
CREATE TRIGGER fill_pipeline_geometry_field
BEFORE INSERT OR UPDATE ON adutoras
FOR EACH ROW EXECUTE PROCEDURE pipeline_geometry_field();

O vídeo a seguir mostra a implementação do código acima, a partir da edição da tabela adutora no QGIS, observa-se o preenchimento automático dos campos após o salvamento da feição no Banco de Dados:

Triggers and Procedures in PostGIS from Marcello Benigno on Vimeo.

PostGIS: Triggers e Procedures (Parte 1)

Quem já estudou normalização em Banco de Dados sabe que campos calculados devem ser eliminados por ferirem a terceira forma normal, no entanto, quando trabalhamos com SIG muitas vezes precisamos de campos como: área, perímetro, centróide, dentre outras medidas, calculadas a partir das geometrias das feições.

Porém, tais campos devem possuir um mecanismo que possibilite o seu cálculo de forma automatizada, tanto para novas feições, quanto para mudanças na forma original da geometria, através da edição da base. Uma forma de implementar tal comportamento em uma tabela espacial com campos calculados, é através de triggers e procedures.

Aqui temos um exemplo uma stored procedure e de uma trigger que calcula os valores de área (area_m2) e perímetro (perimetro_m) da tabela lote, a partir da geometria, armazenada no campo geom:

--Criação dos campos que serão calculados (caso não existam ainda)
ALTER TABLE lotes ADD area_m2 numeric(10,3);
ALTER TABLE lotes ADD perimetro_m numeric(10,3);

--Atualização dos valores dos campos calculados
UPDATE lotes SET area_m2 = ST_Area(geom);
UPDATE lotes SET perimetro_m = ST_Perimeter(geom);

-----------------------------
--Criação da Stored Procedure
-----------------------------
CREATE OR REPLACE FUNCTION calculate_geometry_fields()
RETURNS trigger AS
$$
BEGIN
NEW.area_m2:=ST_Area(NEW.geom);
NEW.perimetro_m:=ST_Perimeter(NEW.geom);
RETURN NEW;
END;
$$
LANGUAGE 'plpgsql';

--------------------
--Criação da Trigger
--------------------
CREATE TRIGGER fill_geometry_fields
BEFORE INSERT OR UPDATE ON lotes
FOR EACH ROW EXECUTE PROCEDURE calculate_geometry_fields();

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.

PostGIS: extração de medidas em tabelas com coordenadas geográficas

É possível obter valores de área, comprimentro, distância e perímetro em tabelas ondes os dados estão em lat/long, através da função Geography(). Seguem alguns exemplos: