[postgis-users] Contours of rasters without or with tiles

[Mathieu Basille]

Dear PostGIS users,

Here is a simple case: I'd like to extract the contour polygon of a given
raster. The dedicated function should be ST_Polygon. If I run it on the
raster imported without tiles, I get one polygon per island/continent,
which is exactly what I want. However, if the raster uses tiles, I get one
polygon per island/continent for each tile (which makes sense). I could
merge them using ST_Union, but I was unable to do that completely: some
polygons do get merged, but not all. This can be fairly problematic in my
case, since I'm using buffers of the resulting polygon, which means I need
one polygon per island/continent.

Example case, with the simple raster attached:

Importation without tiles:

raster2pgsql -s 4326 -I -C -M raster.tif test.raster | psql -h localhost -d
database -U pguser

Importation with tiles:

raster2pgsql -s 4326 -t 5x5 -I -C -M raster.tif test.raster_tile | psql -h
localhost -d database -U pguser

Using ST_Polygon:

CREATE TABLE test.pol AS
SELECT ST_Polygon(rast) FROM test.raster;

CREATE TABLE test.pol_t AS
SELECT ST_Polygon(rast) FROM test.raster_tile;

Using ST_Union on the result:

CREATE TABLE test.pol_tu AS
SELECT ST_Union(ST_Polygon(rast)) FROM test.raster_tile;

See the picture attached to see the result. 'pol' and 'pol_t' give the
expected results, but not 'pol_tu'. Is this an expected behaviour? Is there
a workaround to this problem?

Thanks in advance for any hint!
Mathieu.


# SELECT PostGIS_Full_Version();

postgis_full_version

-----------------------------------------------------------------------
POSTGIS="2.1.0SVN r10597" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23
September 2009" GDAL="GDAL 2.0dev, released 2011/12/29" LIBXML="2.8.0"
LIBJSON="UNKNOWN" TOPOLOGY RASTER


--

~$ whoami
Mathieu Basille, PhD

~$ locate --details
University of Florida \\
Fort Lauderdale Research and Education Center
(+1) 954-577-6314
http://ase-research.org/basille

~$ fortune
� Le tout est de tout dire, et je manque de mots
Et je manque de temps, et je manque d'audace. �
-- Paul �luard

[Francois Hugues]

Hello,

Did you try with the bounding boxes of your tiles using st_envelope ? I operated numerous tiles and unioned them with this function.

Your picture does not show us how many lines your last query returns. Here, it could be a multipolygon.

Hugues.



-----Message d'origine-----
De : postgis-us...@lists.osgeo.org [mailto:postgis-us...@lists.osgeo.org] De la part de Mathieu Basille
Envoyé : vendredi 23 novembre 2012 21:53
À : PostGIS-users
Objet : [postgis-users] Contours of rasters without or with tiles

-- Paul Éluard
_______________________________________________
postgis-users mailing list
postgi...@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[Mathieu Basille]

Dear François,

Thank you for your answer. Answers below, in your mesage:

Le 26/11/2012 12:02, Francois Hugues a écrit :
> Hello,
>
> Did you try with the bounding boxes of your tiles using st_envelope ? I operated numerous tiles and unioned them with this function.

I'm not quite sure to understand your suggestion. I don't see how bounding
boxes could be useful for my problem: What I would like to get is the
polygon surrounding not null pixels of the raster, not bounding boxes of
the tiles. Could you elaborate here please?

> Your picture does not show us how many lines your last query returns. Here, it could be a multipolygon.

Here are the results for the three commands:

SELECT ST_Polygon(rast) FROM test.raster;
1 line

SELECT ST_Polygon(rast) FROM test.raster_tile;
9 lines

SELECT ST_Union(ST_Polygon(rast)) FROM test.raster_tile;
1 line

Indeed, the last one seems to be a multipolygon. I admit that I'm
completely sure about the consequences of that, and how I could solve it.
But I see what you mean. From what I understand, ST_Union is unable to
correctly 'merge' all the polygons of the main land, but only three of
them, but correctly 'union' them. My assumption that ST_Union was wrong was
then maybe misleading.

But the bottom line is that I get different results from tiled and not
tiled rasters, which should not happen IMHO (or I'm just misunderstanding
the rationale of using tiles). I just used ST_Union as an attempt to solve
the behaviour of ST_Polygon, which gives different results on tiled and not
tiled rasters. In my example, I think ST_Polygon should return one line in
both cases with the same output (here a multipolygon with three polygons
for the 2 islands and the main land).
My question ("Is this an expected behaviour?") should thus rather apply to
ST_Polygon.

Mathieu.

[Francois Hugues]

Hello,

There's no mess ! You used my firstname correctly and I thank you for that.

To (try to) answer your questions :
1. I would have use st_envelope to see what was the result and if the problem would comes from no data values. St_envelope will also give you one polygon by tile and the result will be totally different if you use it with a tiled or non tiled raster.
2. As st_envelope, I think the correct behaviour of st_polygon is to proceed tile by tile if you tiled your raster (a tile by row = a polygon by row).

Hugues.


--
Hugues FRANÇOIS
Ingénieur recherche

Irstea Grenoble
2 rue de la papeterie
38400 Saint-Martin d'Hères

tel : +33 (0)4.76.76.27.44
port : +33 (0)6.77.66.21.31
fax : +33 (0)4.76.51.38.03


-----Message d'origine-----
De : Mathieu Basille [mailto:bas...@ase-research.org]
Envoyé : mardi 27 novembre 2012 00:14
À : Francois Hugues
Objet : Re: [postgis-users] Contours of rasters without or with tiles

[sorry for the mess, I mixed up your first and last name, and forgot to CC you... I think it's time to get back home]


Dear Hugues,

Thank you for your answer. Answers below, in your mesage:

Le 26/11/2012 12:02, Francois Hugues a écrit :

> Hello,
>
> Did you try with the bounding boxes of your tiles using st_envelope ? I operated numerous tiles and unioned them with this function.

I'm not quite sure to understand your suggestion. I don't see how bounding
boxes could be useful for my problem: What I would like to get is the
polygon surrounding not null pixels of the raster, not bounding boxes of
the tiles. Could you elaborate here please?

> Your picture does not show us how many lines your last query returns. Here, it could be a multipolygon.

Here are the results for the three commands:

SELECT ST_Polygon(rast) FROM test.raster;
1 line

SELECT ST_Polygon(rast) FROM test.raster_tile;
9 lines

SELECT ST_Union(ST_Polygon(rast)) FROM test.raster_tile;
1 line

Indeed, the last one seems to be a multipolygon. I admit that I'm
completely sure about the consequences of that, and how I could solve it.
But I see what you mean. From what I understand, ST_Union is unable to
correctly 'merge' all the polygons of the main land, but only three of
them, but correctly 'union' them. My assumption that ST_Union was wrong was
then maybe misleading.

But the bottom line is that I get different results from tiled and not
tiled rasters, which should not happen IMHO (or I'm just misunderstanding
the rationale of using tiles). I just used ST_Union as an attempt to solve
the behaviour of ST_Polygon, which gives different results on tiled and not
tiled rasters. In my example, I think ST_Polygon should return one line in
both cases with the same output (here a multipolygon with three polygons
for the 2 islands and the main land).
My question ("Is this an expected behaviour?") should thus rather apply to
ST_Polygon.

Mathieu.

> Hugues.
>
>
>
> -----Message d'origine-----
> De : postgis-us...@lists.osgeo.org [mailto:postgis-us...@lists.osgeo.org] De la part de Mathieu Basille
> Envoyé : vendredi 23 novembre 2012 21:53
> À : PostGIS-users
> Objet : [postgis-users] Contours of rasters without or with tiles
>
> Dear PostGIS users,
>

> Here is a simple case: I'd like to extract the contour polygon of a givenraster. The dedicated function should be ST_Polygon. If I run it on the raster imported without tiles, I get one polygon per island/continent, which is exactly what I want. However, if the raster uses tiles, I get one polygon per island/continent for each tile (which makes sense). I could merge them using ST_Union, but I was unable to do that completely: some polygons do get merged, but not all. This can be fairly problematic in my case, since I'm using buffers of the resulting polygon, which means I need one polygon per island/continent.

>
> Example case, with the simple raster attached:
>
> Importation without tiles:
>
> raster2pgsql -s 4326 -I -C -M raster.tif test.raster | psql -h localhost -d database -U pguser
>
> Importation with tiles:
>
> raster2pgsql -s 4326 -t 5x5 -I -C -M raster.tif test.raster_tile | psql -h localhost -d database -U pguser
>
> Using ST_Polygon:
>
> CREATE TABLE test.pol AS
> SELECT ST_Polygon(rast) FROM test.raster;
>
> CREATE TABLE test.pol_t AS
> SELECT ST_Polygon(rast) FROM test.raster_tile;
>
> Using ST_Union on the result:
>
> CREATE TABLE test.pol_tu AS
> SELECT ST_Union(ST_Polygon(rast)) FROM test.raster_tile;
>

> See the picture attached to see the result. 'pol' and 'pol_t' give the expected results, but not 'pol_tu'. Is this an expected behaviour? Is there aworkaround to this problem?

[Mathieu Basille]

Dear Hugues,

I see now what you mean. ST_Envelope gives the behaviour you described
(i.e. one polygon per tile, which is 9 lines). If I try to ST_Union them, I
get the exact same problem I had before with ST_Polygon: not all
polygon-tiles are unioned, as can be seen with a ST_Dump:

SELECT ST_Dump(ST_Union(ST_Envelope(rast))) FROM test.raster_tile;
st_dump
--------------------------------------------------
({1},0103000020E6100000010000000900000024068...
({2},0103000020E6100000010000000B0000009F1A2...
(2 lines)

With this line of code on the tile version of the raster, one would expect
a single line as output (as is the case on the non-tiled raster).

So my questions are really now:

1) Is this a desirable feature that vector functions such as ST_Envelope or
ST_Polygon behave differently on tiled and non-tiled rasters? From my point
of view, it shouldn't be the case, as the decision to tile a raster should
be transparent from a user point of view (it's just a trick to get better
performance for most GIS procedures). The same code should run fine on both
types of rasters in my opinion.

2) How can we circumvent this "feature"? Using ST_Union should do it, but
it's obviously not working as expected (at least by me).

I feel like I'm getting quite far from my initial request, so maybe it
should be better to start a new thread... Any opinion?

Best,
Mathieu.

Le 27/11/2012 01:43, Francois Hugues a écrit :
> Hello,
>

> There's no mess ! You used my firstname correctly and I thank you for that.
>
> To (try to) answer your questions :
> 1. I would have use st_envelope to see what was the result and if the problem would comes from no data values. St_envelope will also give you one polygon by tile and the result will be totally different if you use it with a tiled or non tiled raster.
> 2. As st_envelope, I think the correct behaviour of st_polygon is to proceed tile by tile if you tiled your raster (a tile by row = a polygon by row).
>
> Hugues.
>
>

--

[Mathieu Basille]

OK, I found some sense out of this... Of course, ST_Union has to be called
directly on the raster *before* ST_Envelope or ST_Polygon... Something like:

SELECT ST_Dump(ST_Envelope(ST_Union(rast))) FROM test.raster_tile;
st_dump
----------------------------------------------------
({},0103000020E610000001000000050000009F1A2FDD24...
(1 line)

However, I wonder whether ST_Envelope/Polygon should not do it
automagically... It can be rather heavy on large raster files, though.

Mathieu.

[Pierre Racine]

Mathieu,

> I see now what you mean. ST_Envelope gives the behaviour you described
> (i.e. one polygon per tile, which is 9 lines). If I try to ST_Union them, I
> get the exact same problem I had before with ST_Polygon: not all
> polygon-tiles are unioned, as can be seen with a ST_Dump:
>
> SELECT ST_Dump(ST_Union(ST_Envelope(rast))) FROM test.raster_tile;
> st_dump
> --------------------------------------------------
> ({1},0103000020E6100000010000000900000024068...
> ({2},0103000020E6100000010000000B0000009F1A2...
> (2 lines)
>
> With this line of code on the tile version of the raster, one would expect
> a single line as output (as is the case on the non-tiled raster).

This is a problem with ST_Union. The resulting polygon SHOULD be a single polygon. I suggest you create a minimal query not involving rasters reproducing the problem and that you submit it back to the list as a geometry only problem. I suggest you get the ST_AsText() version of your envelope and create a simple new self contained SQL query reproducing the problem. This way you will get attention from the right people.

While you get an answer... There are two ways to fix the problem:

1) Snap the geometries returned by ST_Envelop() (or ST_Polygon) to a very fine grid with ST_Snap(). ST_Union should work better.

2) Buffer each geometry returned by ST_Envelop() (or ST_Polygon) a little bit with ST_Buffer() before passing them to ST_Union and then remove the buffer from the result of ST_Union. For this to work properly you have to tell ST_Buffer 'endcap=square join=mitre' as the "buffer_style_parameters".

> 1) Is this a desirable feature that vector functions such as ST_Envelope or
> ST_Polygon behave differently on tiled and non-tiled rasters? From my point
> of view, it shouldn't be the case, as the decision to tile a raster should
> be transparent from a user point of view (it's just a trick to get better
> performance for most GIS procedures). The same code should run fine on both
> types of rasters in my opinion.

This depend on you definition of a raster... Apparently you decided that a whole table IS one raster. In PostGIS, one row is definitely a raster but a whole table can and can not be interpreted as something like a raster coverage and a raster coverage can still be quite different from a simple raster. It can be not rectangular, having missing tiles, overlaps, storing unrelated raster, etc... PostGIS API works on rasters. YOU work on a raster coverage (depending on the way you decided to store you rasters or your raster tiles and hence to interpret the meaning of a table containing rasters). It is then YOUR job to ST_Union the result of a query on a coverage.

So just keep in mind: one RASTER column in a row is a raster. A table containing a RASTER column is what you decide it to be. I agree however that the fact that the loader can "tile" a raster to many rows may let you think that a table is still a raster... I prefer to call that a raster coverage.

Pierre

[Mathieu Basille]

Dear Pierre,

Thanks for the detailed information. This makes perfect sense now. But as
you rightly guessed, I was misled by the loader tiling possibility. With
the loader, the tiling appears to be a mere feature of the raster, just as
the SRID for instance. But in the end, it is not just a feature of the
raster, as it gives a completely different object! Maybe this should be
emphasized in the documentation, to make things clear: tiling a raster can
dramatically improve performance, but this results in a different object (a
raster coverage of many rasters), hence possibly different code -- I insist
on the code, since this is what the user really sees and uses. Given your
explanations, it also makes sense to have different versions of the same
raster layer for different tasks, e.g. one non tiled (i.e. a raster
coverage with only one raster/row) to get the polygon, and one tiled (i.e.
a raster coverage with multiple rasters/tiles/rows) for intersections...

I'll try to repost a reproducible example of the ST_Union problem as soon
as I get some time for it!

Thanks again,
Mathieu.

--

~$ whoami
Mathieu Basille, PhD

~$ locate --details
University of Florida \\
Fort Lauderdale Research and Education Center
(+1) 954-577-6314
http://ase-research.org/basille

~$ fortune
« Le tout est de tout dire, et je manque de mots
Et je manque de temps, et je manque d'audace. »
-- Paul Éluard