Blog Geo.NET Geoprocessamento, SIG e Sensoriamento Remoto

No post anterior, falamos um pouquinho sobre como o WKT Raster funciona. Hoje gostaria de mostrar um pouquinho sobre as maquinações que ando desenvolvendo para o PostGIS e WKT Raster.

Em tese, se o raster está armazenado no banco e conseguimos determinar o valor de cada pixel, poderíamos usar a SQL para realizar álgebra de mapas? A resposta é sim, podemos. É claro, já existem softwares especializados em realizar álgebra de mapas, como o IDRISI, ENVI, SPRING, ArcGIS (com Spatial Analyst), o módulo Sextante do gvSIG, WhiteBox, entre outros.

A vantagem de fazermos isto no PostGIS é que podemos aproveitar de algumas funcionalidades do banco de dados para realizar atualizações automáticas (funciona, mas a perfomance não deve ser das melhores), controlar melhor o acesso aos rasters - centralizando tudo em uma única localização, backup, etc.

Existem diversas funções para análise espacial de rasters, incluindo as matemáticas, lógicas, entre outras. Hoje irei mostrar apenas algumas matemáticas, que são bem simples. O módulo WKT Raster já planeja construir em C (linguagem nativa do PostgreSQL) algumas funções para análise espacial, mas enquanto isto não acontece, vamos levando com estas aqui.

Importante: estas funções só operam com rasters de uma única banda. Todas as operações são realizadas na banda 1. Nenhuma operação é realizada na banda 2 ou 3 - se bem que o suporte para tal é simples de ser construído.

Cuidado ao utilizar estas funções. Caso você sobreescreva seu raster, não o conseguirá de volta - caso precise.

ST_Plus

Esta função adiciona um valor double precision ao nosso raster de interesse. Se o pixel tem valor 1 e nosso argumento é 1, o resultado daquele pixel será 2. Esta função tem como entrada um raster e um valor numérico. O retorno da função é um objeto raster completo.

CREATE OR REPLACE FUNCTION ST_Plus(p_grid raster, z double precision)
	RETURNS raster AS
$$
DECLARE
	w integer;
	h integer;
	z1 integer;

	grid raster;
BEGIN

	SELECT
		ST_Width(p_grid),
		ST_Height(p_grid)
	 into w,h;

	-- realizar cópia do raster original
	grid := p_grid;

	for i in 1..h loop
		for j in 1..w loop

			-- selecting cell value
			z1 := ST_Value(p_grid,i,j);
			grid := st_setvalue(grid,1,i,j,z1 + z);

		end loop;
	end loop;

	return grid;
END
$$
language 'plpgsql';

A função é bem simples. Ela copia o raster original em memória e adiciona o valor à cada pixel. Sem segredos. Como usá-la? Recomendo primeiramente, duplicar sua coluna do tipo raster, pois a função de adição de colunas raster é bastante longa e com parâmetros meio obscuros. Procure o arquivo SQL gerado pelo gdal2wktraster e ache a linha que adiciona a coluna. Modifique o nome da coluna e execute este comando novamente. Depois disso popule a coluna com os novos valores.

SELECT AddRasterColumn('public','teste_raster','grid_stplus',-1, ARRAY['16BUI'], false,
       true, null, 0.000833, -0.000833, 9, 9,
           ST_Envelope(
               ST_SetSRID('POLYGON((-49.500416 -18.000416,-49.500416 -19.000416,-48.000416 -18.000416,-48.000416 -19.000416,-49.500416 -18.000416))'::geometry, -1)));

UPDATE teste_raster SET grid_stplus = ST_Plus('nome_da_coluna_raster_original',10);

Faça um teste com ST_Value para identificar se tudo correu bem. Neste teste, somamos 10 unidades ao valor de cada pixel da coluna original raster. Por que estamos criando uma coluna extra para armazenar nosso novo raster? Se alterarmos o raster original, não conseguiremos ele de volta, portanto, cuidado ao utilizar as funções.

ST_Minus

De forma similar à ST_Plus, ST_Minus vai subtrair determinado valor de cada pixel de nosso raster e nos trará como resultado o raster modificado.

CREATE OR REPLACE FUNCTION ST_Minus(p_grid raster, z double precision)
	RETURNS raster AS
$$
DECLARE
	w integer;
	h integer;
	z1 integer;

	grid raster;
BEGIN

	SELECT
		ST_Width(p_grid),
		ST_Height(p_grid)
	 into w,h;

	-- realizar cópia do raster original
	grid := p_grid;

	for i in 1..h loop
		for j in 1..w loop

			-- selecting cell value
			z1 := ST_Value(p_grid,i,j);
			grid := st_setvalue(grid,1,i,j,z1 - z);

		end loop;
	end loop;

	return grid;
END
$$
language 'plpgsql';

Recomendo testar da mesma maneira que a função anterior. O uso é o mesmo - modificando o nome da função.

ST_Plus

Esta função realiza a multiplicação do valor do pixel. Simples e direta.

CREATE OR REPLACE FUNCTION ST_Times(p_grid raster, z double precision)
	RETURNS raster AS
$$
DECLARE
	w integer;
	h integer;
	z1 integer;

	grid raster;
BEGIN

	SELECT
		ST_Width(p_grid),
		ST_Height(p_grid)
	 into w,h;

	-- realizar cópia do raster original
	grid := p_grid;

	for i in 1..h loop
		for j in 1..w loop

			-- selecting cell value
			z1 := ST_Value(p_grid,i,j);
			grid := st_setvalue(grid,1,i,j,z1 * z);

		end loop;
	end loop;

	return grid;
END
$$
language 'plpgsql';

Use-a do mesmo jeito que as anteriores.

ST_Divide

CREATE OR REPLACE FUNCTION ST_Divide(p_grid raster, z double precision)
	RETURNS raster AS
$$
DECLARE
	w integer;
	h integer;
	z1 integer;

	grid raster;
BEGIN
	if (z = 0) then
		return p_grid;
	end if;

	SELECT
		ST_Width(p_grid),
		ST_Height(p_grid)
	 into w,h;

	-- realizar cópia do raster original
	grid := p_grid;

	for i in 1..h loop
		for j in 1..w loop

			-- selecting cell value
			z1 := ST_Value(p_grid,i,j);
			grid := st_setvalue(grid,1,i,j,z1 / z);

		end loop;
	end loop;

	return grid;
END
$$
language 'plpgsql';

ST_Power

Esta função eleva o valor da célula ao valor passado como parâmetro na função. Se temos uma célula de valor 2, e passamos 2 como parâmetro, teremos o 4 como resultado (2^2);

CREATE OR REPLACE FUNCTION ST_Power(p_grid raster, z double precision)
	RETURNS raster AS
$$
DECLARE
	w integer;
	h integer;
	z1 integer;

	grid raster;
BEGIN

	SELECT
		ST_Width(p_grid),
		ST_Height(p_grid)
	 into w,h;

	-- realizar cópia do raster original
	grid := p_grid;

	for i in 1..h loop
		for j in 1..w loop

			-- selecting cell value
			if z = 0 then
				grid := st_setvalue(grid,1,i,j,1);
			else
				z1 := st_value(p_grid,i,j);
				grid := st_setvalue(grid,1,i,j,z1 ^ z);
			end if;

		end loop;
	end loop;

	return grid;
END
$$
language 'plpgsql';

E assim podemos construir funções que alterem nossos rasters originais. São funções simples, praticamente demonstrativas. Na sequência demonstraremos algumas funções mais interessantes, como declividade, cálculo de fluxo de direção, identificação de depressões, eliminação de depressões, entre outras.

Um abraço,

George R. C. Silva