[raster] MapAlgebraExpr can't handle conditions

[nicklas]

When running the example for the one raster version of ​MapAlgebraExpr:

UPDATE dummy_rast SET map_rast2 = ST_MapAlgebraExpr(rast,'2BUI','CASE WHEN [rast] BETWEEN 100 and 250 THEN 1 
WHEN [rast] = 252 THEN 2 
WHEN [rast] BETWEEN 253 and 254 THEN 3 ELSE 0 END', '0') WHERE rid = 2;

I don't get the values 1, 2, 3 and 0.

When I run, as in the example.

SELECT DISTINCT ST_Value(rast,1,i,j) As origval, ST_Value(map_rast2, 1, i, j) As mapval
FROM dummy_rast CROSS JOIN generate_series(1, 5) AS i CROSS JOIN generate_series(1,5) AS j
WHERE rid = 2;

I get:

origval;mapval
249;
250;
251;
252;
253;
254;

It is beta1

[nicklas]

It is the same with 2 raster version. At once a condition is involved it gives nothing (or maybe I had 1 as result in one case)

[Bborie Park]

Strange that this isn't working. Concise example:

DROP TYPE IF EXISTS cellset CASCADE;
CREATE TYPE cellset AS (
	x integer,
	y integer,
	value double precision
);

CREATE OR REPLACE FUNCTION dump_cells(rast raster, band int DEFAULT 1)
	RETURNS SETOF cellset
	AS $$
	DECLARE
		r cellset%rowtype;
		rows int;
		cols int;
	BEGIN
		SELECT width, height INTO rows, cols FROM ST_Metadata(rast);

		FOR r IN SELECT x, y, ST_Value(rast, band, x, y)
			FROM generate_series(1, rows) AS x
			CROSS JOIN generate_series(1, cols) AS y
			LOOP

			RETURN NEXT r;

		END LOOP;
	END;
	$$
	LANGUAGE 'plpgsql';

WITH foo AS (
	SELECT
		ST_AddBand(
			ST_MakeEmptyRaster(5, 5, 0, 0, 1, -1, 0, 0, 0)
			, '8BUI'::text, 1, 0
		)
	AS rast
)
SELECT dump_cells(rast)
FROM (
	SELECT
		ST_MapAlgebraExpr(rast, '8BUI', 'CASE WHEN [rast.x] % 2 = 0 AND [rast.y] % 2 = 0 THEN NULL ELSE 128. END') AS rast
	FROM foo
) bar

[Bborie Park]

As for two raster MapAlgebraExpr, it works as expected...

DROP TYPE IF EXISTS cellset CASCADE;
CREATE TYPE cellset AS (
	x integer,
	y integer,
	value double precision
);

CREATE OR REPLACE FUNCTION dump_cells(rast raster, band int DEFAULT 1)
	RETURNS SETOF cellset
	AS $$
	DECLARE
		r cellset%rowtype;
		rows int;
		cols int;
	BEGIN
		SELECT width, height INTO rows, cols FROM ST_Metadata(rast);

		FOR r IN SELECT x, y, ST_Value(rast, band, x, y)
			FROM generate_series(1, rows) AS x
			CROSS JOIN generate_series(1, cols) AS y
			LOOP

			RETURN NEXT r;

		END LOOP;
	END;
	$$
	LANGUAGE 'plpgsql';

WITH foo AS (
	SELECT
		ST_AddBand(
			ST_MakeEmptyRaster(5, 5, 0, 0, 1, -1, 0, 0, 0)
			, '8BUI'::text, 1, 0
		)
	AS rast
)
SELECT dump_cells(rast)
FROM (
	SELECT
		ST_MapAlgebraExpr(rast, rast, 'CASE WHEN [rast1.x] % 2 = 0 AND [rast1.y] % 2 = 0 AND [rast2.x] % 2 = 0 AND [rast2.y] % 2 = 0 THEN NULL ELSE 128. END') AS rast
	FROM foo
) bar

[Bborie Park]

By the looks of it, the solution for the 1-raster function is to cast the result from the prepared statement into a double precision.

Fixed in r9371. If there is still present after r9371, please provide a concise and portable example of the failure.

[nicklas]

Yes, I can confirm that 2-rasterworks as expected.

The reason I said it didn't was that I tested to make an expression like this:

(([rast1]-[rast2])>50)::int

and expected to get 0 and 1 for true or false as I do in a normal SELECT query. Bu that didn't work.

I also tried to populate pixeltype '1BB' with the above without the cast to integer but that didn't work either.

How to populate boolean pixeltype?

[Bborie Park]

The following works for me using your expression...

DROP TABLE IF EXISTS test_mapalgebra;
CREATE TABLE test_mapalgebra(
	rid serial PRIMARY KEY,
	rast raster
);
CREATE OR REPLACE FUNCTION make_test_raster()
	RETURNS void
	AS $$
	DECLARE
		width int := 5;
		height int := 5;
		x int;
		y int;
		rast raster;
	BEGIN
		rast := ST_MakeEmptyRaster(width, height, 0, 0, 1, -1, 0, 0);
		rast := ST_AddBand(rast, 1, '8BUI', 0, 0);

		FOR x IN 1..width LOOP
			FOR y IN 1..height LOOP
				rast := ST_SetValue(rast, 1, x, y, random() * 256);
			END LOOP;
		END LOOP;

		INSERT INTO test_mapalgebra (rast) VALUES (rast);

		RETURN;
	END;
	$$ LANGUAGE 'plpgsql';
SELECT make_test_raster();
SELECT make_test_raster();
DROP FUNCTION make_test_raster();

WITH foo AS (
	SELECT
		ST_MapAlgebraExpr(a.rast, b.rast, '(([rast1]-[rast2]) > 50)::int', '1BB') AS rast
	FROM test_mapalgebra a
	CROSS JOIN test_mapalgebra b
	WHERE a.rid = 1
		AND b.rid = 2
)
SELECT
	(ST_BandMetadata(rast)).*,
	(ST_SummaryStats(rast)).*
FROM foo;
DROP TABLE IF EXISTS test_mapalgebra;

[Bborie Park]

Just tested the example in my last comment using ​http://postgis.refractions.net/download/postgis-2.0.0beta2SVN.tar.gz and it works.

[Bborie Park]

Beta1 works as well for the example above.

[nicklas]

I don't think it is working. The output is wrong, at least from my box.

I get true in all cases or false in all cases.

I also get the same result if I send the result to '8BUI' pixeltype so the problem seems to be in the expression not the cast to boolean.

What casts is actually happening? Why do you cast to integer in your example when the result from the expression is boolean and the result pixeltype is boolean. The error-message if I don't do that cast is that it can't be casted to double precision.

Can it be all those casts that is making things bogus after all?

[nicklas]

I forgot to say that you should look in the attached sql-file too where I try to boil down the problem.

[pracine]

dustymug: You're funny with your dump_cells function. Didb't you give a try to ST_PixelAsPolygons

​http://postgis.refractions.net/documentation/manual-svn/RT_ST_PixelAsPolygons.html

[pracine]

False interpretation:

"So both pixels returns true as I understand it."

No, one pixel gets 1 (true) and the other get nodata. As yo can see the count is 1, not two.

"But at once I get under 8, instead of giving one false and one true I get both true."

As soon as your condition gets under 8, one pixel gets true = 1 and the other gets nodata value, hence the correct result.

I think we can close this one :-) One less!

[nicklas]

Great Pierre

almost there

First 1-raster version is still not working as I understand dustymugs comments above.

Second:

You are right I misunderstood things, but still, the casts doesn't make sense.

If I cast to Integer and put it in '8BSI' like this:

WITH foo AS (
	SELECT
		ST_MapAlgebraExpr(a.rast, b.rast, '(([rast1]-[rast2]) > 5)::int', '8BSI') AS rast
	FROM test_mapalgebra a
	CROSS JOIN test_mapalgebra b
	WHERE a.rid = 1
		AND b.rid = 2
)
SELECT
	(ST_BandMetadata(rast)).*,
	(ST_SummaryStats(rast)).*
FROM foo;

Then I would expect the false results to show as 0, right?

But I still get the same result (nodata as I have learned now :-) )

I also wonder if something is wrong because of performance. I am using this expressions to find clear cuttings (I will give a PostGIS talk next week and will show the functionality)

To do this calculation in PostGIS it takes about 45 seconds, but in QGIS with the same images as geotiff it takes less than 2 seconds.

But it can be because PostGIS can't gain fully from the SATA3 ssd on this box while maybe the file based gdal calculation through QGIS can?

Or is it some casts that is giving a false bottleneck?

[pracine]

Replying to nicklas:

If I cast to Integer and put it in '8BSI' like this: Then I would expect the false results to show as 0, right?

Yes, but still 0 is nodata.

But I still get the same result (nodata as I have learned now :-) )

Seems normal to me.

To do this calculation in PostGIS it takes about 45 seconds, but in QGIS with the same images as geotiff it takes less than 2 seconds.

I guess a real comparative test with QGIS would be to open the image and make the calculation but yes, we still have a lot of work to do on our map algebra to compare with QGIS or Arc.

But it can be because PostGIS can't gain fully from the SATA3 ssd on this box while maybe the file based gdal calculation through QGIS can?

I'm afraid it won't make much difference...

Or is it some casts that is giving a false bottleneck?

We still have to find the bottleneck but I guess it is in the great flexibility of our expression parser. With flexibility comes slowiness. We will have to think about the avenues we have to make it faster. I think there are many.

[Bborie Park]

Replying to pracine:

dustymug: You're funny with your dump_cells function. Didb't you give a try to ST_PixelAsPolygons

​http://postgis.refractions.net/documentation/manual-svn/RT_ST_PixelAsPolygons.html

Heeheehee. I didn't realize that ST_PixelAsPolygons did what I want... need to read the docs more often!

[Bborie Park]

The expression "(([rast1]-[rast2]) > 5)::int" must be cast from boolean to integer first because a cast from boolean to double precision is not permitted. An error was thrown by PostgreSQL 9.1 when I attempted that.

[Bborie Park]

Last month, it was found that some casts (integer vs double precision) caused noticeable bottlenecks when executing prepared statements. We haven't had time to ask the PostgreSQL crowd why using integer instead of the usual double precision would cause a significant slowdown.

I don't expect PostGIS Raster to significantly beat any dedicated desktop raster processing application ever. If we're lucky, we'll be able to match the speed though realistically we'll probably be happy incremental performance boosts over time.

[pracine]

To close or not to close?

[Bborie Park]

+1 for closing. There was a genuine issue with 1-raster ST_MapAlgebraExpr where the returning value from the expression had to be a double precision. Anything else being returned from the expression resulted in NULL.

The output I'm seeing from my examples above look correct and are behaving as expected. I'll try nicklas' attached example and see if the output matches what the expression indicates.

[pracine]

We'll wait for ok from Nicklas for closing.

[nicklas]

Yes I am convinced

Thanks a lot :-)

I missed the fix before the examples.

A side note.

From the expression we get a logical true or false, then that is cast to integer because PostgreSQL otherwise won't let us cast it to double precision. So we can put it in a 1 bit pixel ...

I guess this is a result of the flexibility Pierre was talking about, but it seems like a long way for a boolean.

Thanks again Nicklas

[nicklas]

side note part 2:

Does that mean that each pixel gets stored as double precision during the process?

[Bborie Park]

In the C code, all pixels are expressed and worked on as double precision values. When it comes time to store a value, the value is clamped/truncated to band's pixel type. I suppose that when the raster project started that a union variable could have been used for working with a pixel but that adds its own bag of nightmares.

And yes, it is does look like a long way to go for a simple boolean. Sadly, the choices are simplicity or flexibility and we picked flexibility.