Re: [postgis-users] Slow MapAlgebra performance

[Pierre Racine]

Hi Kirill,

Using the callback version of ST_MapAlgebra() should be much faster but you 
have to write your own callback function.

https://postgis.net/docs/RT_ST_MapAlgebra.html

Let us know the difference.

Pierre

De : postgis-users  De la part de ? 
?? ??
Envoyé : April 6, 2020 2:57 PM
À : postgis-users@lists.osgeo.org
Objet : [postgis-users] Slow MapAlgebra performance

[Externe UL*]
Dear All!



I calculate vegetation indices on rasters (Landsat 8 scenes).

The calculation of SAVI=(NIR-RED)/(NIR+RED+L)*(1+L) takes about 2 minutes to 
complete.

This is enormous time compared to the same operation in QGIS which takes 
several seconds.



Here is the query (it’s defined as postgres UDF and $1 parameter is L):



BEGIN

DROP TABLE IF EXISTS savi;

CREATE TABLE savi AS

SELECT bandt5.rid, ST_MapAlgebra(

bandt4.rast,

bandt5.rast,

'CASE WHEN ([rast2.val] + [rast1.val] = 0) THEN NULL ELSE (([rast2.val] - 
[rast1.val])*(1+$1) / ([rast2.val] + [rast1.val]+$1)::float) END', '32BF'

 )

FROM bandt5, bandt4 WHERE bandt5.rid=bandt4.rid;

RETURN 1;

END;



I used the data from USGS, Landsat 8 Level 2, bands are stored in separate 
tables.

DB indices are created for every band by raster2pgsql.

Tile size is 128x128.

The whole raster2pgsql command: for each tif band:

"raster2pgsql -s 4326 -N -32767 -t 128x128 -I -C -M -d bandX.tif public.bandtX 
> bandX.sql".

Afterwards bandX.sql files were ingested by psql.



Are there any ways to improve the performance?



Thank you!



---

Regards,

Kirill


*ATTENTION : L’émetteur de ce courriel est externe à l’Université Laval.
Évitez de cliquer sur un hyperlien, d’ouvrir une pièce jointe ou de transmettre 
des informations si vous ne connaissez pas l’expéditeur du courriel. En cas de 
doute, contactez l’équipe de soutien informatique de votre unité ou 
hameconn...@ulaval.ca.


___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Point geom, ST_Value, ST_Intersects and out of range raster coordinates

[Pierre Racine]

Solution is to filter out NULL returned when using ST_Value() or DISTINCT or 
GROUP BY multiple values returned when using ST_Intersect().

Using (only) one solution consistently returns consistent results.

Pierre

-Message d'origine-
De : postgis-users  De la part de Peter 
Devoy
Envoyé : 8 janvier 2020 09:29
À : PostGIS Users Discussion 
Objet : [postgis-users] Point geom, ST_Value, ST_Intersects and out of range 
raster coordinates

[Externe UL*]

Hi list

I am trying to sample a value in a raster using ST_Value and a point geometry.
With the same point geometry I tried to use ST_Intersects to sample from only 
the relevant row in the raster table but if the point is on the boundary of 
multiple raster tiles a row is returned for each.

However, at the cell level, ST_Value seems to play by different rules and raises
"NOTICE: Attempting to get pixel value with out of range raster coordinates..."
for all but one of the rows.  Resulting in NULL values for all but one row.

Please does anyone know if there is there a method or function I can use to get 
more "consistent" spatial logic without modifying the input point geometry?

E.g. either to filter out all rows where the sampled cell would be out of 
range; or, to have ST_Value return a value for each row found by ST_Intersects.

My postgis_version() = "2.5 USE_GEOS=1 USE_PROJ=1 USE_STATS=1".

Kind regards


Peter
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users
*ATTENTION : L’émetteur de ce courriel est externe à l’Université Laval.
Évitez de cliquer sur un hyperlien, d’ouvrir une pièce jointe ou de transmettre 
des informations si vous ne connaissez pas l’expéditeur du courriel. En cas de 
doute, contactez l’équipe de soutien informatique de votre unité ou 
hameconn...@ulaval.ca.


___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] raster2pgsql fail with JPEG encoded image

[Pierre Racine]

GDAL is complaining that your image is rotated and that it does not support it 
yet. What if you erase (or just rename) the .jgw, temporally to remove the 
rotation?

Pierre

-Message d'origine-
De : postgis-users  De la part de Jorge 
Gustavo Rocha
Envoyé : 27 septembre 2019 18:18
À : PostGIS Users Discussion 
Objet : [postgis-users] raster2pgsql fail with JPEG encoded image

Hi Postgisers,

I'm trying to find out why a particular JPEG image does not work with 
raster2pgsql.

I've described the problem as a question in gis.stackexchange:

https://gis.stackexchange.com/questions/337127/raster2pgsql-fail-with-jpeg-encoded-image

There are no errors when inserting the raster in the database.
gdalinfo works fine, both with the original image and the Postgis raster.

The only error I got is if I try to save the Postgis raster back to a file 
system image, like:

gdal_translate -of JPEG PG:"service=pg_geotuga dbname=geotuga schema=raster 
table=geologica mode=2" 05.jpg

When I do this, I get the error: "GDAL PostGIS Raster driver can not work with 
rotated rasters yet." (but the image is written, but I'm not able to render it).

I suspect that the encoding/compression of this specific image is not supported.

identify '05-D Braga.jpg'
05-D Braga.jpg JPEG 10614x7739 10614x7739+0+0 8-bit sRGB 11.1453MiB 0.000u 
0:00.000

This is just one of the many geological rasters of the country. I would like to 
upload them all to Postgis. I can recode the images, if necessary, but even 
changing them from JPEG to GTiff is not enough.

Preferably, I would like to identify and maybe fix the issue in raster2pgsql.

The GDAL version reported by Postgis is 2.3.2 (does it use static libs?). The 
GDAL version on the system is 3.0.1.

raster2pgsql -?
RELEASE: 2.5.2 GDAL_VERSION=23 (r17328)

gdalinfo --version
GDAL 3.0.1, released 2019/06/28

select postgis_full_version();
POSTGIS="2.5.2 r17328" [EXTENSION] PGSQL="110" GEOS="3.7.0-CAPI-1.11.0 
673b9939" PROJ="Rel. 5.2.0, September 15th, 2018" GDAL="GDAL 2.3.2, released 
2018/09/21" LIBXML="2.9.4" LIBJSON="0.12.1" LIBPROTOBUF="1.3.1"
RASTER

The database has raster support enabled and I can upload other rasters to the 
database without any issues.

Best regards,

Jorge Gustavo
--
Jorge Gustavo Rocha
Departamento de Informática
Universidade do Minho
4710-057 Braga
Gabinete 3.29 (Piso 3)
Tel: +351 253604480
Fax: +351 253604471
Móvel: +351 910333888
skype: nabocudnosor
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] raster regions

[Pierre Racine]

Right. To minimize the amount of intermediate computation you can add a 
discriminant to the first query:

WITH floodareas AS (
  SELECT ST_MapAlgebra(rast, ‘1BB’, ‘CASE WHEN [rast] < 10 THEN 1 ELSE NULL 
END’) rast
  FROM elevation e, districts d
  WHERE ST_Intersects(d.geom, e.rast)
)
SELECT d.*
FROM districts d, floodareas f
WHERE ST_Intersects(d.geom, f.rast)

Pierre
De : postgis-users  De la part de 
joao.corde...@inpe.br
Envoyé : 15 juillet 2019 14:40
À : postgis-users@lists.osgeo.org
Objet : Re: [postgis-users] raster regions


Thanks, Pierre,

It seems ok to model the query the way you showed, however, that implies the 
generation of an entire intermediate binary raster, just to identify the region 
inside the whole mapped area. My issue is on avoiding such intermediate data by 
incorporating map algebra expression into the same query expression, Another 
example, ""... d.geom touch (ndvi2018 < ndvi2019)" where ndvi are raster 
data.

Att.

Em 15-07-2019 14:43, Pierre Racine escreveu:
Hi,

That could translate into something like this in PostGIS:

WITH floodareas AS (
  SELECT ST_MapAlgebra(rast, ‘1BB’, ‘CASE WHEN [rast] < 10 THEN 1 ELSE NULL 
END’) rast
  FROM elevation
)
SELECT d.*
FROM districts d, floodareas f
WHERE ST_Intersects(g.geom, f.rast)

NULL becomes nodata values. You could also use ST_Reclass() instead of 
ST_Mapalgebra().

Pierre

De : postgis-users 
mailto:postgis-users-boun...@lists.osgeo.org>>
 De la part de Joao Cordeiro
Envoyé : 15 juillet 2019 10:19
À : 'PostGIS Users Discussion' 
mailto:postgis-users@lists.osgeo.org>>
Objet : [postgis-users] raster regions


Hi there,

I am not actually a POSTGIS user, or developer, although I am doing some 
research of integration of vector -raster in a dbms context. My main question s 
about identifying regions represented in a domain like the "hybrid raster" so 
that they can be referred to in a query. For example, I would like to write 
something like "select d from districts where d.geometry touches flooding", 
where "flooding" indicates an expression intended to decide on the pertinence 
of locations to a flooding ​​risk area described, for instance through map 
algebra.

Thanks for any hints or suggestions.

Att.


___

postgis-users mailing list

postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>

https://lists.osgeo.org/mailman/listinfo/postgis-users



___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] raster regions

[Pierre Racine]

Hi,

That could translate into something like this in PostGIS:

WITH floodareas AS (
  SELECT ST_MapAlgebra(rast, ‘1BB’, ‘CASE WHEN [rast] < 10 THEN 1 ELSE NULL 
END’) rast
  FROM elevation
)
SELECT d.*
FROM districts d, floodareas f
WHERE ST_Intersects(g.geom, f.rast)

NULL becomes nodata values. You could also use ST_Reclass() instead of 
ST_Mapalgebra().

Pierre

De : postgis-users  De la part de Joao 
Cordeiro
Envoyé : 15 juillet 2019 10:19
À : 'PostGIS Users Discussion' 
Objet : [postgis-users] raster regions


Hi there,

I am not actually a POSTGIS user, or developer, although I am doing some 
research of integration of vector -raster in a dbms context. My main question s 
about identifying regions represented in a domain like the "hybrid raster" so 
that they can be referred to in a query. For example, I would like to write 
something like "select d from districts where d.geometry touches flooding", 
where "flooding" indicates an expression intended to decide on the pertinence 
of locations to a flooding ​​risk area described, for instance through map 
algebra.

Thanks for any hints or suggestions.

Att.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Strange behavior with empty geometries GROUP BY

[Pierre Racine]

:)

Well if you have some geometric aggregators you have to expect that at some 
point some crazy people will want to GROUP BY or DISTINCT geometries. A simple 
example (not crazy at all) is if I want to search for duplicates geometries:

SELECT geom, count(*) cnt
FROM mytable
GROUP BY geom
HAVING count(*) > 1;

Now I know it doesn't work well enough since two different geometries having 
the same bounding box will look like duplicates. I should probably GROUP BY 
ST_AsEWKB(geom) instead:

SELECT ST_AsEWKB(geom)::geometry, count(*) cnt
FROM mytable
GROUP BY ST_AsEWKB(geom)
HAVING count(*) > 1;

I guess that a good replacement?

Pierre

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 11 septembre 2017 16:32
À : 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Cc : 'PostGIS Development Discussion' <postgis-de...@lists.osgeo.org>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Pierre,

Jump on the boat of all those folks who want it changed.  You are in crowded 
company.
Only reason I think  we haven't done it is it will require a lot of testing to 
see how it breaks peoples code, any performance bottlenecks, and people relying 
on that odd behavior for well over 15 years.

I personally don't care as why anyone would want to group by something as 
largish as a geometry always puzzled me
And no definition of geometric equality ever satisfied me to the point where I 
wouldn't be asked questions

How come ST_Equals says these are equal but   geom::text = geom::text says they 
are not.

Seems to me like trading one set of problems for another.

Thanks,
Regina
From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Monday, September 11, 2017 4:20 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Any reason why it's not something more precise then the equality of the 
bounding boxes? Is it redefinable?

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 11 septembre 2017 15:59
À : 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

BTree = operator

Which is bounding box equality

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Monday, September 11, 2017 3:54 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Strange behavior with empty geometries GROUP BY

What is the operator used by GROUP BY?

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 11 septembre 2017 15:42
À : 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Hmm interesting.  I was about to say they probably have the same bounding box 
but they don't quite according to box2d.  Perhaps our box2d returns float4 but 
our internal box is double precision?

I don't know.
--

WITH polys AS (
  SELECT 1 id, ST_GeomFromText('MULTIPOLYGON(((251000 5243597.61702128,
   251038 5243558,
   251056 5243527,
   251057.227272727 5243500,
   251000 5243500,
   251000 5243597.61702128)),
 ((251021 5243690,
   251000 5243681.87096774,
   251000 5243750,
   251250 5243750,
   251250 5243749.86206897,
   251196 5243694,
   251168 5243687,
   251021 5243690)))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((251000 5243681.87096774,
 251021 5243690,
 251168 5243687,
 251196 5243694,
 251250 5243749.86206897,
 251250 5243598.15,
 251167.008064516 5243500,
 251057.227272727 5243500,
25105

Re: [postgis-users] Strange behavior with empty geometries GROUP BY

[Pierre Racine]

Any reason why it's not something more precise then the equality of the 
bounding boxes? Is it redefinable?

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 11 septembre 2017 15:59
À : 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

BTree = operator

Which is bounding box equality

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Monday, September 11, 2017 3:54 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Strange behavior with empty geometries GROUP BY

What is the operator used by GROUP BY?

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 11 septembre 2017 15:42
À : 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Hmm interesting.  I was about to say they probably have the same bounding box 
but they don't quite according to box2d.  Perhaps our box2d returns float4 but 
our internal box is double precision?

I don't know.
--

WITH polys AS (
  SELECT 1 id, ST_GeomFromText('MULTIPOLYGON(((251000 5243597.61702128,
   251038 5243558,
   251056 5243527,
   251057.227272727 5243500,
   251000 5243500,
   251000 5243597.61702128)),
 ((251021 5243690,
   251000 5243681.87096774,
   251000 5243750,
   251250 5243750,
   251250 5243749.86206897,
   251196 5243694,
   251168 5243687,
   251021 5243690)))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((251000 5243681.87096774,
 251021 5243690,
 251168 5243687,
 251196 5243694,
 251250 5243749.86206897,
 251250 5243598.15,
 251167.008064516 5243500,
 251057.227272727 5243500,
251056 5243527,
 251038 5243558,
 251000 5243597.61702128,
 251000 5243681.87096774))') geom
)
SELECT BOX2D(geom)
FROM polys;


BOX(251000 5243500,251250 5243750)
BOX(251000 5243500,251250 5243749.86206897)




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Monday, September 11, 2017 3:08 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Strange behavior with empty geometries GROUP BY

I get a similar bug but now with full pledged geometries. Those two polygons, 
though very different, GROUP BY together. Which mysterious PostGIS operator is 
used by GROUP BY and make those two polygons equivalent?

  WITH polys AS (
  SELECT 1 id, ST_GeomFromText('MULTIPOLYGON(((251000 5243597.61702128,
   251038 5243558,
   251056 5243527,
   251057.227272727 5243500,
   251000 5243500,
   251000 5243597.61702128)),
 ((251021 5243690,
   251000 5243681.87096774,
   251000 5243750,
   251250 5243750,
   251250 5243749.86206897,
   251196 5243694,
   251168 5243687,
   251021 5243690)))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((251000 5243681.87096774,
 251021 5243690,
 251168 5243687,
 251196 5243694,
  

Re: [postgis-users] Strange behavior with empty geometries GROUP BY

[Pierre Racine]

What is the operator used by GROUP BY?

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 11 septembre 2017 15:42
À : 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Hmm interesting.  I was about to say they probably have the same bounding box 
but they don't quite according to box2d.  Perhaps our box2d returns float4 but 
our internal box is double precision?

I don't know.
--

WITH polys AS (
  SELECT 1 id, ST_GeomFromText('MULTIPOLYGON(((251000 5243597.61702128,
   251038 5243558,
   251056 5243527,
   251057.227272727 5243500,
   251000 5243500,
   251000 5243597.61702128)),
 ((251021 5243690,
   251000 5243681.87096774,
   251000 5243750,
   251250 5243750,
   251250 5243749.86206897,
   251196 5243694,
   251168 5243687,
   251021 5243690)))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((251000 5243681.87096774,
 251021 5243690,
 251168 5243687,
 251196 5243694,
 251250 5243749.86206897,
 251250 5243598.15,
 251167.008064516 5243500,
 251057.227272727 5243500,
251056 5243527,
 251038 5243558,
 251000 5243597.61702128,
 251000 5243681.87096774))') geom
)
SELECT BOX2D(geom)
FROM polys;


BOX(251000 5243500,251250 5243750)
BOX(251000 5243500,251250 5243749.86206897)




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Monday, September 11, 2017 3:08 PM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Strange behavior with empty geometries GROUP BY

I get a similar bug but now with full pledged geometries. Those two polygons, 
though very different, GROUP BY together. Which mysterious PostGIS operator is 
used by GROUP BY and make those two polygons equivalent?

  WITH polys AS (
  SELECT 1 id, ST_GeomFromText('MULTIPOLYGON(((251000 5243597.61702128,
   251038 5243558,
   251056 5243527,
   251057.227272727 5243500,
   251000 5243500,
   251000 5243597.61702128)),
 ((251021 5243690,
   251000 5243681.87096774,
   251000 5243750,
   251250 5243750,
   251250 5243749.86206897,
   251196 5243694,
   251168 5243687,
   251021 5243690)))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((251000 5243681.87096774,
 251021 5243690,
 251168 5243687,
 251196 5243694,
 251250 5243749.86206897,
 251250 5243598.15,
 251167.008064516 5243500,
 251057.227272727 5243500,
 251056 5243527,
 251038 5243558,
 251000 5243597.61702128,
 251000 5243681.87096774))') geom
) SELECT ST_AsText(geom)
  FROM polys
  GROUP BY geom;

Pierre

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Pierre Racine
Envoyé : 29 juin 2017 15:42
À : PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Objet : Re: [p

Re: [postgis-users] Strange behavior with empty geometries GROUP BY

[Pierre Racine]

I get a similar bug but now with full pledged geometries. Those two polygons, 
though very different, GROUP BY together. Which mysterious PostGIS operator is 
used by GROUP BY and make those two polygons equivalent?

  WITH polys AS (
  SELECT 1 id, ST_GeomFromText('MULTIPOLYGON(((251000 5243597.61702128,
   251038 5243558,
   251056 5243527,
   251057.227272727 5243500,
   251000 5243500,
   251000 5243597.61702128)),
 ((251021 5243690,
   251000 5243681.87096774,
   251000 5243750,
   251250 5243750,
   251250 5243749.86206897,
   251196 5243694,
   251168 5243687,
   251021 5243690)))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((251000 5243681.87096774,
 251021 5243690,
 251168 5243687,
 251196 5243694,
 251250 5243749.86206897,
 251250 5243598.15,
 251167.008064516 5243500,
 251057.227272727 5243500,
 251056 5243527,
 251038 5243558,
 251000 5243597.61702128,
 251000 5243681.87096774))') geom
) SELECT ST_AsText(geom)
  FROM polys
  GROUP BY geom;

Pierre

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Pierre Racine
Envoyé : 29 juin 2017 15:42
À : PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Much better...

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 23 juin 2017 23:26
À : 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Pierre,

Thanks for pointing out the issue.  Sandro has fixed the problem

https://trac.osgeo.org/postgis/ticket/3777#comment:17


Can you try with your larger workload.  If you happen to have a windows box, 
64-bit  pre-compiled binaries here:

http://postgis.net/windows_downloads/

(I have 32-bit building turned off at moment cause it errors out on the gui 
tests and haven't had a chance to troubleshoot).


If you have a patched version:

SELECT postgis_full_version();

Your version should have  r15450 (or above) for PostGIS 2.4 and r15451 (or 
above)for PostGIS 2.3.

POSTGIS="2.3.3dev r15451" GEOS="3.6.1-CAPI-1.10.1 r4317" PROJ="Rel. 4.9.1, 04 
March 2015" GDAL="GDAL 2.1.3, released 2017/20/01" LIBXML="2.7.8" 
LIBJSON="0.12" RASTER

Thanks,
Regina

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Tuesday, June 20, 2017 5:42 PM
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: [postgis-users] Strange behavior with empty geometries GROUP BY

Hi,

When I try to count the occurrence of a number of empty geometries like this:

SELECT ST_AsText(geom) txt, count(*)
FROM (
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 1)'::geometry geom
UNION ALL
SELECT 'LINESTRING(0 0,0 1)'::geometry geom
UNION ALL
SELECT 'GEOMETRYCOLLECTION EMPTY'::geometry geom
UNION ALL
SELECT 'POINT EMPTY'::geometry geom
) foo
GROUP BY geom;

The 'GEOMETRYCOLLECTION EMPTY' and the 'POINT EMPTY' do not aggregate and I get:

txt 
count
POINT(0 0) 2
POINT EMPTY   1
POINT(0 0) 1
LINESTRING(0 0,0 1)  1
GEOMETRYCOLLECTION EMPTY1
POINT(0 1) 1

If I remove any geometry other than the two EMPTY ones like this (I removed the 
first one 'POINT(0 0)'), the two EMPTY geometries DO aggregate:

SELECT ST_AsText(geom) txt, count(*)
FROM (
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT

Re: [postgis-users] ST_Polygonize give different results on same set of linestrings

[Pierre Racine]

I simplified the example a bit, removing ST_ExteriorRing():

WITH polys AS (
  SELECT ST_GeomFromText('LINESTRING(1 1, 3 2, 3 0, 1 1)') geom
  UNION ALL
  SELECT ST_GeomFromText('LINESTRING(2 1, 4 2, 4 0, 2 1)') geom
), extrings AS (
SELECT 1 test, geom
FROM polys
UNION ALL
SELECT 2 test, ST_Union(geom) geom
FROM polys
UNION ALL
SELECT 3 test, ST_Collect(geom) geom
FROM polys
)
SELECT test, ST_AsText(ST_Polygonize(geom))
FROM extrings
GROUP BY test;

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Pierre Racine
Envoyé : 7 septembre 2017 09:46
À : PostGIS Users Discussion <postgis-users@lists.osgeo.org>
Objet : [postgis-users] ST_Polygonize give different results on same set of 
linestrings

Hi,

Is it normal that ST_Polygonize() gives two different results depending on if 
the linestrings are ST_Union() before calling it or not?

In the following example, test 1 (no union) and 3 (ST_Collect) give the same 
result, constructing two overlapping polygons. Test 2 (ST_Union) on the other 
side constructs 3 disjoint polygons removing overlaps. All tests provide, 
however, the same basic set of linestrings to ST_Polygonize.

WITH polys AS (
  SELECT ST_GeomFromText('POLYGON((1 1, 3 2, 3 0, 1 1))') geom
  UNION ALL
  SELECT ST_GeomFromText('POLYGON((2 1, 4 2, 4 0, 2 1))') geom
), extrings AS (
SELECT 1 test, ST_ExteriorRing(geom) geom
FROM polys
UNION ALL
SELECT 2 test, ST_Union(ST_ExteriorRing(geom)) geom
FROM polys
UNION ALL
SELECT 3 test, ST_Collect(ST_ExteriorRing(geom)) geom
FROM polys
)
SELECT test, ST_AsText(ST_Polygonize(geom))
FROM extrings
GROUP BY test;

It seems to me that only test 2 (ST_Union) returns the expected result which is 
to reconstruct planar polygons from a set of linestrings however they are 
provided.

It bothers me because I would like to get the same result as test 2 without 
relying on ST_Union() which seems to sometime slightly shift vertexes in order 
to union geometries together (true or false?).

Thanks for any explanation.

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] ST_Polygonize give different results on same set of linestrings

[Pierre Racine]

Hi,

Is it normal that ST_Polygonize() gives two different results depending on if 
the linestrings are ST_Union() before calling it or not?

In the following example, test 1 (no union) and 3 (ST_Collect) give the same 
result, constructing two overlapping polygons. Test 2 (ST_Union) on the other 
side constructs 3 disjoint polygons removing overlaps. All tests provide, 
however, the same basic set of linestrings to ST_Polygonize.

WITH polys AS (
  SELECT ST_GeomFromText('POLYGON((1 1, 3 2, 3 0, 1 1))') geom
  UNION ALL
  SELECT ST_GeomFromText('POLYGON((2 1, 4 2, 4 0, 2 1))') geom
), extrings AS (
SELECT 1 test, ST_ExteriorRing(geom) geom
FROM polys
UNION ALL
SELECT 2 test, ST_Union(ST_ExteriorRing(geom)) geom
FROM polys
UNION ALL
SELECT 3 test, ST_Collect(ST_ExteriorRing(geom)) geom
FROM polys
)
SELECT test, ST_AsText(ST_Polygonize(geom))
FROM extrings
GROUP BY test;

It seems to me that only test 2 (ST_Union) returns the expected result which is 
to reconstruct planar polygons from a set of linestrings however they are 
provided.

It bothers me because I would like to get the same result as test 2 without 
relying on ST_Union() which seems to sometime slightly shift vertexes in order 
to union geometries together (true or false?).

Thanks for any explanation.

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Strange behavior with empty geometries GROUP BY

[Pierre Racine]

Much better...

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Regina Obe
Envoyé : 23 juin 2017 23:26
À : 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Objet : Re: [postgis-users] Strange behavior with empty geometries GROUP BY

Pierre,

Thanks for pointing out the issue.  Sandro has fixed the problem

https://trac.osgeo.org/postgis/ticket/3777#comment:17


Can you try with your larger workload.  If you happen to have a windows box, 
64-bit  pre-compiled binaries here:

http://postgis.net/windows_downloads/

(I have 32-bit building turned off at moment cause it errors out on the gui 
tests and haven't had a chance to troubleshoot).


If you have a patched version:

SELECT postgis_full_version();

Your version should have  r15450 (or above) for PostGIS 2.4 and r15451 (or 
above)for PostGIS 2.3.

POSTGIS="2.3.3dev r15451" GEOS="3.6.1-CAPI-1.10.1 r4317" PROJ="Rel. 4.9.1, 04 
March 2015" GDAL="GDAL 2.1.3, released 2017/20/01" LIBXML="2.7.8" 
LIBJSON="0.12" RASTER

Thanks,
Regina

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: Tuesday, June 20, 2017 5:42 PM
To: 'PostGIS Users Discussion' 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: [postgis-users] Strange behavior with empty geometries GROUP BY

Hi,

When I try to count the occurrence of a number of empty geometries like this:

SELECT ST_AsText(geom) txt, count(*)
FROM (
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 1)'::geometry geom
UNION ALL
SELECT 'LINESTRING(0 0,0 1)'::geometry geom
UNION ALL
SELECT 'GEOMETRYCOLLECTION EMPTY'::geometry geom
UNION ALL
SELECT 'POINT EMPTY'::geometry geom
) foo
GROUP BY geom;

The 'GEOMETRYCOLLECTION EMPTY' and the 'POINT EMPTY' do not aggregate and I get:

txt 
count
POINT(0 0) 2
POINT EMPTY   1
POINT(0 0) 1
LINESTRING(0 0,0 1)  1
GEOMETRYCOLLECTION EMPTY1
POINT(0 1) 1

If I remove any geometry other than the two EMPTY ones like this (I removed the 
first one 'POINT(0 0)'), the two EMPTY geometries DO aggregate:

SELECT ST_AsText(geom) txt, count(*)
FROM (
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 1)'::geometry geom
UNION ALL
SELECT 'LINESTRING(0 0,0 1)'::geometry geom
UNION ALL
SELECT 'GEOMETRYCOLLECTION EMPTY'::geometry geom
UNION ALL
SELECT 'POINT EMPTY'::geometry geom
) foo
GROUP BY geom;

Result:

txt  count
POINT(0 0)  2
LINESTRING(0 0,0 1)   1
POINT EMPTY2
POINT(0 1)  1

Any idea why empty geometries aggregation seem so dependent on what others 
geoms are in the table?

Sorry I could not build a more simple example.

Thanks,

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] PostGIS ST_EuclideanDistance() not working. Wrong upperleftx and y?

[Pierre Racine]

Try setting the ST_MakeEmptyRaster() upperlefty parameter to a negative value 
(-1).

Pierre

De : postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] De la part de 
Vitor Sapucaia
Envoyé : 21 juin 2017 04:58
À : postgis-users@lists.osgeo.org
Objet : [postgis-users] PostGIS ST_EuclideanDistance() not working. Wrong 
upperleftx and y?

Hello guys!

I need help on this question that I've posted on GIS stackexchange:
https://gis.stackexchange.com/questions/244622/postgis-st-euclideandistance-not-working-wrong-upperleftx-and-y?noredirect=1#comment383015_244622

here is the question:


I'm trying to make this ST_EuclideanDistance() function work, but no success so 
far. I'm not getting the expected result because the raster seems to be drawn 
at a wrong location.

Documentation:

https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools

Firstly, I'm considering this test page to test the function:

https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools/test

When running the second test case:

CREATE TABLE test_source_raster_1 AS

SELECT

1 AS rid, rast

  FROM

(SELECT

   ST_MakeEmptyRaster(10,10,0,0,1,1,0,0,4326) AS rast

 ) foo;

CREATE TABLE test_source_geometry_1 AS

SELECT

generate_series(1,10) AS id,

ST_RandomPoints(the_geom,10) AS the_geom

  FROM

(SELECT

   ST_SetSRID(ST_Extent(rast::geometry), 4326) AS the_geom

FROM test_source_raster_1) foo;

CREATE TABLE test_ref_raster_1 AS

SELECT

  1 AS rid, rast

  FROM

(SELECT

   ST_MakeEmptyRaster(10,10,0,0,1,1,0,0,4326) AS rast

) foo;

 CREATE TABLE test_eudist_1_a_i_1 AS

   (SELECT

  1 AS rid,

  ST_EuclideanDistance(

 rast,'32BF','public',

 'test_source_geometry_1', 'the_geom', True

   ) AS rast

 FROM test_ref_raster_1);

I'm getting this result (every raster pixel is null):

[results]

And this is the expected result I should have gotten:

[expected]
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

[postgis-users] Strange behavior with empty geometries GROUP BY

[Pierre Racine]

Hi,

When I try to count the occurrence of a number of empty geometries like this:

SELECT ST_AsText(geom) txt, count(*)
FROM (
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 1)'::geometry geom
UNION ALL
SELECT 'LINESTRING(0 0,0 1)'::geometry geom
UNION ALL
SELECT 'GEOMETRYCOLLECTION EMPTY'::geometry geom
UNION ALL
SELECT 'POINT EMPTY'::geometry geom
) foo
GROUP BY geom;

The 'GEOMETRYCOLLECTION EMPTY' and the 'POINT EMPTY' do not aggregate and I get:

txt 
count
POINT(0 0) 2
POINT EMPTY   1
POINT(0 0) 1
LINESTRING(0 0,0 1)  1
GEOMETRYCOLLECTION EMPTY1
POINT(0 1) 1

If I remove any geometry other than the two EMPTY ones like this (I removed the 
first one 'POINT(0 0)'), the two EMPTY geometries DO aggregate:

SELECT ST_AsText(geom) txt, count(*)
FROM (
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 0)'::geometry geom
UNION ALL
SELECT 'POINT(0 1)'::geometry geom
UNION ALL
SELECT 'LINESTRING(0 0,0 1)'::geometry geom
UNION ALL
SELECT 'GEOMETRYCOLLECTION EMPTY'::geometry geom
UNION ALL
SELECT 'POINT EMPTY'::geometry geom
) foo
GROUP BY geom;

Result:

txt  count
POINT(0 0)  2
LINESTRING(0 0,0 1)   1
POINT EMPTY2
POINT(0 1)  1

Any idea why empty geometries aggregation seem so dependent on what others 
geoms are in the table?

Sorry I could not build a more simple example.

Thanks,

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Fwd: PostGIS Out-DB Raster Not Behaving As Expected

[Pierre Racine]

I guess st_transform should be applied only when using the same syntax as 
shp2pgsql : -s from_srid:to_srid.

And it should be prohibited when used with –R. I don’t think we want to go into 
reprojecting the tiles in the filesystem. The only good way to do this is 
convert them to in-db and it is not what the users is expecting. So better 
prevent it. If he really want to reproject, he can do it after registering the 
out-db tiles, creating new in-db ones with SELECT ST_Transform(). At least at 
this moment it is clearer that the band will become in-db.

Make sence?

Pierre

De : Regina Obe [mailto:l...@pcorp.us]
Envoyé : 16 avril 2017 00:06
À : 'PostGIS Users Discussion' <postgis-users@lists.osgeo.org>
Cc : Pierre Racine <pierre.rac...@sbf.ulaval.ca>; 'Bborie Park' 
<dustym...@gmail.com>
Objet : RE: [postgis-users] Fwd: PostGIS Out-DB Raster Not Behaving As Expected

Thanks for bringing this to our attention.

I've ticketed it - https://trac.osgeo.org/postgis/ticket/3738

Pierre and Bborie,

Were you aware of this behavior?  This is the first time I'm hearing of it.
Though that would explain one complaint I heard from someone that their 
raster2pgsql outdb took 10 times longer than in-db
When for me it's always been really fast.

I would consider this a bug that –s could be used for both setting the SRID and 
also transforming data and would like the behavior changed in 2.4.
Probably too late to change for 2.3 since it's an api change, but we can put a 
warning on at least 2.3 as Tumasgiu suggested.

Thanks,
Regina
http://postgis.us
http://postgis.net





From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Tumasgiu Rossini
Sent: Tuesday, April 11, 2017 8:32 AM
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Fwd: PostGIS Out-DB Raster Not Behaving As Expected

Maybe a warning showing up when using -s and -R options together could be a 
nice thing to do ?
Or maybe implementing another option which would use st_setsrid instead of 
st_transform.

2017-04-11 14:25 GMT+02:00 Tumasgiu Rossini 
<rossin...@gmail.com<mailto:rossin...@gmail.com>>:
Hi,
I think I found the probelm, it is the -s option which make the sql
applying a st_transfrom on the raster to reproject it into the desired CRS
before inserting it in the table, thus creating an in-db raster.

2017-04-11 12:14 GMT+02:00 Osahon Oduware 
<osahon@gmail.com<mailto:osahon@gmail.com>>:
Hi Regina,

I would like to add that the SHOW command below returns "on":


SHOW postgis.enable_outdb_rasters;
which I believe indicates that the environment variable is SET to enable out-db 
raster. However, the result from your 1st query shows that the raster table 
isn't created with out-db. Any reason for this behaviour?


-- Forwarded message --
From: Osahon Oduware <osahon@gmail.com<mailto:osahon@gmail.com>>
Date: Tue, Apr 11, 2017 at 10:07 AM
Subject: Re: [postgis-users] PostGIS Out-DB Raster Not Behaving As Expected
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Cc: Regina Obe <l...@pcorp.us<mailto:l...@pcorp.us>>, Jorge Gustavo Rocha 
<j...@di.uminho.pt<mailto:j...@di.uminho.pt>>
Hi Regina,

I have tried the queries above. For the related raster table,

SELECT r_table_name, r_raster_column, out_db FROM raster_columns; RETURNS
"raster_table","rast","[False, False, False]"

SELECT t.rid, (md).isoutdb, (md).path
FROM nj_ned AS t, ST_BandMetaData(t.rast) AS md
limit 1; RETURNS
1,False,""

I believe from the result of the 1st query that out-db is not enabled. However, 
the environment variable POSTGIS_ENABLE_OUTDB_RASTERS is SET to 1. I have tried 
to also enable it with the following, restarted the postgres Service, and 
re-loaded the raster, but the 1st query gives the same result:

ALTER DATABASE  SET postgis.enable_outdb_rasters TO True;


Please, could you help me further on what I am not doing right to enable the 
out-db rasters.




On Mon, Apr 10, 2017 at 8:31 PM, Regina Obe 
<l...@pcorp.us<mailto:l...@pcorp.us>> wrote:
You should also check the following queries

SELECT r_table_name, r_raster_column, out_db FROM raster_columns;

Out_db should be true for all your bands for the specific table and as Jorge 
mentioned, you shouldn't see any over view tables (which start with o_..{table}.
I believe QGIS internally tries to use over view tables to speed things up, so 
it's possible it's using an overview table instead of the main outdb table if 
you are zoomed out enough.

For outdb, querying the meta data of the raster will tell you the path if it's 
an outdb.  Replace nj_ned with your table name.

SELECT t.rid, (md).isoutdb, (md).path
FROM nj_ned AS t, ST_BandMetaData(t.rast) AS md
limit 1;

Hope that helps,
R

Re: [postgis-users] Cost weighted distance raster

[Pierre Racine]

ST_Raster2WorldCoordX and ST_Raster2WorldCoordY were replaced with 
ST_RasterToWorldCoordX and RasterToWorldCoordY in newer version of PostGIS.

Just replace those names in the function and reexecute the CREATE FUNCTION…

Pierre

From: Spencer Gardner [mailto:spencergard...@gmail.com]
Sent: 16 mars 2017 15:50
To: Pierre Racine <pierre.rac...@sbf.ulaval.ca>
Cc: postgis-users@lists.osgeo.org
Subject: Re: Cost weighted distance raster

Just to follow up:

I was able to execute the CREATE FUNCTION scripts with some minor tweaking. 
When I run ST_CostDistance it chokes on making a call to ST_World2RasterCoordx. 
That's as far as I've gotten. I'll have to turn my attention to another 
solution for now.

It's a shame--having a good cost and raster routing library would be very nice 
to have.

Thanks for your help.

On Tue, Mar 7, 2017 at 11:59 AM, Spencer Gardner 
<spencergard...@gmail.com<mailto:spencergard...@gmail.com>> wrote:
Ah thanks. Didn't see the code page. I'll let you know whether these work as is 
or if they need tweaking.

On Tue, Mar 7, 2017 at 11:51 AM, Pierre Racine 
<pierre.rac...@sbf.ulaval.ca<mailto:pierre.rac...@sbf.ulaval.ca>> wrote:
In the code page :

https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools/code

There is a ST_CostDistance() function.

Pierre




___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Cost weighted distance raster

[Pierre Racine]

In the code page :

https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools/code

There is a ST_CostDistance() function.

Pierre

From: Spencer Gardner [mailto:spencergard...@gmail.com]
Sent: 7 mars 2017 12:35
To: Pierre Racine <pierre.rac...@sbf.ulaval.ca>
Cc: postgis-users@lists.osgeo.org
Subject: Re: Cost weighted distance raster

I didn't see any code for the Cost Distance function, only the Euclidean 
Distance. Did I miss something? Or is it located on a different website that I 
haven't found?

Thanks,
Spencer

On Tue, Mar 7, 2017 at 10:23 AM, Pierre Racine 
<pierre.rac...@sbf.ulaval.ca<mailto:pierre.rac...@sbf.ulaval.ca>> wrote:
There is no .sql. just copy the functions into a .sql text file and execute.

From: Spencer Gardner 
[mailto:spencergard...@gmail.com<mailto:spencergard...@gmail.com>]
Sent: 7 mars 2017 11:16
To: Pierre Racine 
<pierre.rac...@sbf.ulaval.ca<mailto:pierre.rac...@sbf.ulaval.ca>>
Cc: postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>
Subject: Re: Cost weighted distance raster

On Tue, Mar 7, 2017 at 7:07 AM, Pierre Racine 
<pierre.rac...@sbf.ulaval.ca<mailto:pierre.rac...@sbf.ulaval.ca>> wrote:
Hi Spencer,

Integrating something in the main PostGIS core is hard work. That’s why I 
started the PostGIS Addons: to quickly gather PLpgSQL-only functions into an 
easy to install .SQL file. Why I did not integrate those functions in the 
Addons yet? Lack of time, not much work involving PostGIS since a while, not 
sure the code is working properly.

Did you get time to test it? If it work fine, I might take time to add it to 
the Addons before FOSS4G 2017.

Pierre

That makes sense. I don't see a master .sql file to use for adding the 
functions to the database. The GSOC proposal I linked to previously has code 
for the Euclidean Distance function, but not for the Cost-Weighted Distance 
function. Is there a single .sql file somewhere that I can use to create both 
functions? I'd love to take a look and see if I can get them both working.


___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Cost weighted distance raster

[Pierre Racine]

There is no .sql. just copy the functions into a .sql text file and execute.

From: Spencer Gardner [mailto:spencergard...@gmail.com]
Sent: 7 mars 2017 11:16
To: Pierre Racine <pierre.rac...@sbf.ulaval.ca>
Cc: postgis-users@lists.osgeo.org
Subject: Re: Cost weighted distance raster

On Tue, Mar 7, 2017 at 7:07 AM, Pierre Racine 
<pierre.rac...@sbf.ulaval.ca<mailto:pierre.rac...@sbf.ulaval.ca>> wrote:
Hi Spencer,

Integrating something in the main PostGIS core is hard work. That’s why I 
started the PostGIS Addons: to quickly gather PLpgSQL-only functions into an 
easy to install .SQL file. Why I did not integrate those functions in the 
Addons yet? Lack of time, not much work involving PostGIS since a while, not 
sure the code is working properly.

Did you get time to test it? If it work fine, I might take time to add it to 
the Addons before FOSS4G 2017.

Pierre

That makes sense. I don't see a master .sql file to use for adding the 
functions to the database. The GSOC proposal I linked to previously has code 
for the Euclidean Distance function, but not for the Cost-Weighted Distance 
function. Is there a single .sql file somewhere that I can use to create both 
functions? I'd love to take a look and see if I can get them both working.

___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Cost weighted distance raster

[Pierre Racine]

Hi Spencer,

Integrating something in the main PostGIS core is hard work. That’s why I 
started the PostGIS Addons: to quickly gather PLpgSQL-only functions into an 
easy to install .SQL file. Why I did not integrate those functions in the 
Addons yet? Lack of time, not much work involving PostGIS since a while, not 
sure the code is working properly.

Did you get time to test it? If it work fine, I might take time to add it to 
the Addons before FOSS4G 2017.

Pierre

From: Spencer Gardner [mailto:spencergard...@gmail.com]
Sent: 6 mars 2017 18:09
To: postgis-users@lists.osgeo.org; Pierre Racine <pierre.rac...@sbf.ulaval.ca>
Subject: Re: Cost weighted distance raster


> Hi Spencer,

> I don’t know if they are still compatible with the newest PostGIS version.



> They used to work fine but don’t expect very good performance.



> Pierre
Thanks, Pierre. Did all that work stagnate? Have there been efforts to 
incorporate into the main PostGIS code? Just curious to know what happened. It 
would be a shame not to provide such useful functionality if there has already 
been substantial work done.

My alternative at the moment is GRASS, which I've not found to be terribly 
flexible for incorporating into a larger workflow.

Spencer
___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Cost weighted distance raster

[Pierre Racine]

Hi Spencer,

I don’t know if they are still compatible with the newest PostGIS version.

They used to work fine but don’t expect very good performance.

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Spencer Gardner
Sent: 3 mars 2017 15:03
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] Cost weighted distance raster

I'm looking to work with PostGIS raster to complete some network analysis using 
a cost raster surface. I came across this GSOC project from 2012 which was 
supposed to create functions to accomplish this:
https://trac.osgeo.org/postgis/wiki/PostGIS_Raster_SoC_Idea_2012/Distance_Analysis_Tools

Does anyone know the status of this work? I don't see the functions in the 
reference documentation for PostGIS so it would appear this work has not been 
folded into the main project. The bottom of the GSOC report includes code to 
create new functions for Euclidean and Cost-Weighted distance, so I'm assuming 
I can just add those to my DB. Can anyone confirm that these functions work as 
expected?

Any information is appreciated.

Thanks,
Spencer


___
postgis-users mailing list
postgis-users@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Raster time-series data

[Pierre Racine]

Sounds reasonable to me. You could try splitting the process between your 4 or 
8 processors, each importing a selection of files. Maybe this could help, maybe 
this could make the whole process even slower (depend on what is slow the CPU 
or the disk)...

If you don't necessarily need all the pixel values in the DB, try the -R option 
when importing. This is definitely faster.

Pierre

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf
> Of Tom Cook
> Sent: Wednesday, August 31, 2016 12:46 PM
> To: postgis-users@lists.osgeo.org
> Subject: [postgis-users] Raster time-series data
> 
> I'm trying to import (what I think of as) large timeseries of gridded data in
> postgis 2.0 (on PostgreSQL.  The data comes as HDF-EOS files (basically HDF4).
> Each file has the whole grid, which is 22680 points, and 24 bands, one for 
> each
> hour of the day.  Each file covers one day, and I'm trying to import 16 years
> (5,844 files).
> 
> My strategy at present is to put each day into a raster with 24 bands, really
> because this is the easiest to implement.  So for the first file:
> 
> raster2pgsql -I -t auto -c 'HDF4_EOS:EOS_GRID:"file1.hdf":EOSGRID:SWGDN'
> swgdn | psql
> psql -c 'alter table swgdn add raster_date date;
> psql -c 'update swgdn set raster_date = '2101';
> 
> And then for subsequent files:
> 
> raster2pgsql -I -t auto -a 'HDF4_EOS:EOS_GRID:"fileX.hdf":EOSGRID:SWGDN'
> swgdn | psql
> psql -c 'update swgdn set raster_date = '' where raster_date is null;
> 
> In a word, it's slow.  I've so far been running the import script for about 
> an hour
> and it's processed 146 input files.  I can't really quantify this, but it 
> feels like it's
> getting slower.
> 
> Is this a reasonable strategy for storing this data, and my performance
> expectations just unrealistic?  Or is there a better structure to use for 
> this?
> 
> Thanks for any suggestions,
> Tom
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Merge rasters with multiple bands

[Pierre Racine]

ST_Union works on overlapping rasters as well but maybe you have to provide two 
3 bands rasters. If your rasters are tiled, you also have to GROUP BY tiles 
having the same alignment:

GROUP BY ST_UpperleftX(rast), ST_UpperleftY(rast)

ST_MapAlgebra() should work also. The expression should be equal to the raster 
with values and the extenttype to the 256x256 pixel one. They shoulp overlap.

Pierre


> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf
> Of Alexander Rumanovsk
> Sent: Wednesday, August 10, 2016 12:40 PM
> To: postgis-users@lists.osgeo.org
> Subject: [postgis-users] Merge rasters with multiple bands
> 
> Hello all,
> 
> I'm trying to merge a raster with 3 bands to an empty raster. The empty 
> raster is
> created with the query bellow:
> 
> ST_AddBand(ST_MakeEmptyRaster(256,256,0,0,1), '8BUI'::text,200)
> 
> 
> This empty raster is larger than the filled one (which is 100x100). I just 
> need to
> create a raster of 256x256 with the original raster (100x100) information and
> the rest filled with empty information.
> My first idea was to use ST_Union, but found out it just merges tiles, 
> instead of
> rasters.
> Anyone have an idea how can I do that?
> 
> thanks in advance
> Rumanovsk
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Slow query ST_SummaryStatsAgg over tiles

[Pierre Racine]

From the ST_SummaryStatsAgg () doc: "By default will sample all pixels. To get 
faster response, set sample_percent to value between 0 and 1"

Pierre

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf
> Of matteo rivola
> Sent: Wednesday, August 10, 2016 11:31 AM
> To: postgis-users@lists.osgeo.org
> Subject: [postgis-users] Slow query ST_SummaryStatsAgg over tiles
> 
> Dear all,
> I'm using postigis version 2.2.1, on PostgreSQL 9.2.0.
> 
> I loaded SRTM DEM geotiff with raster2pgsql.
> I create a query in order to retrieve the global statistics of the tiles, 
> over only
> Europe zone.
> 
> The query is:
> 
> SELECT ((stats).min) as value from (SELECT
> ST_SummaryStatsAgg(ST_Clip(.rast,ST_GeomFromText(,ST_SRI
> D(.rast)),true),1,true) as stats from ) as stats;
> 
> where DB_table is the table in which are the rasters and wkt is equal to
> 'POLYGON((-16 74,39 74,39 33,-16 33,-16 74))'.
> 
> The query works fine but it takes a lot of time.
> 
> Is there a way to speed up the query?
> 
> Cheers,
> Matteo
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Raster tile size impact results

[Pierre Racine]

Sorry for the late answer. 

You have to vectorize them on the fly from within OpenJump:

SELECT (ST_DumpAsPolygons(rast).* FROM yourrastertable;

If your raster is very big and tiled, just vectorize the tiles you need:

SELECT rid, rast::geometry FROM yourrastertable;

Look at the rids and then:

SELECT (ST_DumpAsPolygons(rast).* FROM yourrastertable
WHERE rid = x or rid = y

Pierre

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf
> Of Pierre Bails
> Sent: Friday, June 24, 2016 12:32 PM
> To: PostGIS Users Discussion <postgis-users@lists.osgeo.org>
> Subject: Re: [postgis-users] Raster tile size impact results
> 
> Hello,
> 
> Unfortunatly, I don't find a way to print raster results in OpenJump. Can you
> explain me how I can do it ?
> 
> Thank you for your help !
> 
> Pierre
> 
> Le lun. 20 juin 2016 à 16:16, Pierre Racine <pierre.rac...@sbf.ulaval.ca
> <mailto:pierre.rac...@sbf.ulaval.ca> > a écrit :
> 
> 
>   I would try to understand, visually using OpenJump, how your tiles are
> clipped in the two cases... I guess sometimes 1 pixel wide tiles get clipped 
> in a
> strange way.
> 
>   Pierre
> 
>   > -Original Message-
>   > From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org
> <mailto:postgis-users-boun...@lists.osgeo.org> ] On Behalf
>   > Of Pierre Bails
>   > Sent: Wednesday, June 15, 2016 11:41 AM
>   > To: postgis-users@lists.osgeo.org <mailto:postgis-
> us...@lists.osgeo.org>
>   > Subject: [postgis-users] Raster tile size impact results
>   >
>   > Hi all,
>   >
>   > I'm looking for get all pixel's value in a raster insert thanks to
> raster2pgsql
>   > intersect by a polygon.
>   > Nevertheless, there are some things that I don't understand: depends
> on tile
>   > size, I don't have the same result... The error between 2 raster is 
> also
> correlate
>   > to the pixel size (I assume it's due to the pixel number).
>   >
>   > Can you explain me how the tile size can influence the final result ?
>   >
>   > Here is the query :
>   > SELECT DISTINCT SUM((px).val)
>   > FROM(
>   >   SELECT *
>   >   FROM  mns INNER JOIN
>   >
> st_setsrid(st_geomfromgeojson('{"type":"Polygon","coordinates":[[[x,y],[x,y],[x
>   > ,y] ]]}'), 4326) as geom on ST_Intersects(mns.rast, geom)
>   > ) AS clip, ST_PixelAsCentroids (ST_Clip(clip.rast, clip.geom),1) AS 
> px;
>   >
>   >
>   > DSM px size 70cm:
>   > raster2pgsql -a -t 250x250 -s 4326 -F pathToDSM.tiff dsmTable (result:
> 46 347
>   > 846)
>   >
>   > raster2pgsql -a -t 50x50 -s 4326 -F pathToDSM.tiff dsmTable: 46 
> 382
> 873
>   > Difference: 35027
>   >
>   >
>   > DSM px size 20cm :
>   > raster2pgsql -a -t 250x250 -s 4326 -F pathToDSM.tiff dsmTable :  567
> 505 888
>   >
>   > raster2pgsql -a -t 50x50 -s 4326 -F pathToDSM.tiff dsmTable: 596
> 897 610
>   > Difference: 29 391 722
>   >
>   >
>   > If we considere 250x250's result as reference, there is a difference 
> of
> 0.07% in
>   > the first case, and 5% in the second. If the polygon area increases, 
> the
> error rate
>   > increases.
>   >
>   > Thank you,
>   >
>   > Pierre
>   ___
>   postgis-users mailing list
>   postgis-users@lists.osgeo.org <mailto:postgis-users@lists.osgeo.org>
>   http://lists.osgeo.org/mailman/listinfo/postgis-users

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] raster question

[Pierre Racine]

Tile them!  Should be much faster.

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of Stephen Crawford
> Sent: Thursday, March 10, 2016 12:05 PM
> To: PostGIS Users Discussion
> Subject: [postgis-users] raster question
> 
> Hell all,
> 
> I have a table of rasters, daily weather data, each record is one full
> raster and an date column, all the same extent.  One of my use cases is
> to drill down to get the all the data for a month or year one grid
> cell.  Performance is slow.  Should I make tiles of these daily rasters?
> I guess I was under the impression that tiles should be used on a "one
> table per raster" basis, and that "many rasters per table" cannot have
> tiles.
> 
> Any thoughts?
> 
> Thanks,
> Steve
> 
> --
> Stephen Crawford
> Center for Environmental Informatics
> The Pennsylvania State University
> 
> 
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Raster SQL to separate high value areas

[Pierre Racine]

Or ST_ReClass() and ST_PixelAsCentroids()?

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of Regina Obe
> Sent: Wednesday, March 09, 2016 7:04 AM
> To: 'PostGIS Users Discussion'
> Subject: Re: [postgis-users] Raster SQL to separate high value areas
> 
> You probably want to use ST_Reclass  and then ST_DumpAsPolygons.
> 
> 
> 
> http://postgis.net/docs/manual-2.2/RT_ST_Reclass.html - this will reduce
> your pixels into a binary of above 70% , not above 70%
> 
> 
> 
> 
> 
> http://postgis.net/docs/manual-2.2/RT_ST_DumpAsPolygons.html would
> then give you two sets of polygons with the two new pixel values.
> 
> 
> 
> 
> 
> 
> 
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of lig...@pobox.com
> Sent: Tuesday, March 08, 2016 3:10 PM
> To: postgis-users@lists.osgeo.org
> Subject: [postgis-users] Raster SQL to separate high value areas
> 
> 
> 
> Hi,
> 
> 
> 
> I am trying to separate areas of high value of lightning probability (above
> 70%) in a Raster and calculate their centroids. I successfully loaded a 
> GeoTiff
> raster image (one band, gray scale, 982x1108 pixels) in PostgreSQL 9.5
> (Postgis 2.2.1). I see the SQL results in OenJump 1.9.
> 
> 
> 
> I used the next SQL, but the function ST_PixelAsPoints result in a geometry
> points which are of the pixel's upper-left corner. I am not sure that is the
> right way to do, because of that displacement. Can I do without this
> transformation raster-geometry ?
> 
> 
> 
> Thanks.
> 
> 
> 
> SELECT row_number() over () AS id,
> 
>   ST_NumGeometries(gc),
> 
>   gc AS geom_collection,
> 
>   ST_Centroid(gc) AS centroid
> 
> FROM (
> 
>   SELECT unnest(ST_ClusterWithin(geom, 0.1)) AS gc
> 
>   FROM (SELECT (ST_PixelAsPoints(rast, 1)).*
> 
> FROM tb_lightningprob) foo
> 
> WHERE val > 0.7
> 
> ) f;
> 
> 
> 
> 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Create raster from row-col table

[Pierre Racine]

Hi Etienne,

> Still, I can't figure how to load an arbitrary number of bands using your
> snippet. If I remove WHERE band=1 I still get a unique band.

This query will work for many bands and will aggregate all of them in a single 
raster:

WITH rows AS (
  SELECT band, colx, array_agg(v ORDER BY rowy) a
  FROM test.measure
  GROUP BY band, colx
), valarray AS (
  SELECT band, array_agg_mult(ARRAY[a] ORDER BY colx) a
  FROM rows
  GROUP BY band
), raster AS (
  SELECT rast FROM test.rasters WHERE rid=1
), newraster AS (
  SELECT ST_AddBand(NULL, array_agg(ST_SetValues(rast, 1, 1, 1, valarray.a) 
ORDER BY band), 1) As rast 
  FROM valarray, raster
)
SELECT band, (ST_PixelAsPolygons(rast, band)).*
FROM newraster, generate_series(1,2) band;

this will work only if they all have the same width/height.

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] ST_Mean4ma , ST_StdDev4ma rounding down

[Pierre Racine]

The two NOTICE are because there is actually two rasters created by in the 
dummy_rast table and the first one has no bands.

Values are rounded by ST_MapAlgebra() because, by default, ST_MapAlgebra() will 
try to produce a raster with the same pixel type as the input rasters. If you 
want the resulting raster to receive float values you must tell ST_MapAlgebra() 
by passing it '32BF' as pixel type.

Try this:

select (st_dumpvalues(st_mapalgebra(array[row(rast, 1), row(rast, 2), row(rast, 
3)]::rastbandarg[], 'st_mean4ma(double precision[], int[], 
text[])'::regprocedure, '32BF'))).* from (select st_union(rast) as rast from 
dummy_rast) as foo;

Pierre

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of Severin Thaler
> Sent: Saturday, February 20, 2016 10:26 PM
> To: PostGIS Users Discussion
> Subject: [postgis-users] ST_Mean4ma , ST_StdDev4ma rounding down
> 
> Hi all,
> 
> here’s the problem.
> after creating the dummy raster : http://postgis.net/docs/manual-
> 2.1/RT_reference.html
> i ran:
> 
> 
> select (st_dumpvalues(rast)).* from dummy_rast;
> 
> and get:
> 
> NOTICE:  Raster provided has no bands
> NOTICE:  Raster provided has no bands
>  nband |valarray
> ---+-
> 
>  1 |
> {{253,254,253,254,254},{253,254,254,253,249},{250,254,254,252,249},{25
> 1,253,254,254,253},{252,250,254,254,254}}
>  2 |
> {{78,98,122,173,209},{96,118,180,249,254},{99,112,169,245,254},{89,99,1
> 22,176,229},{79,88,97,112,135}}
>  3 |
> {{70,86,100,135,161},{80,108,162,227,250},{90,108,175,251,254},{77,86,1
> 09,164,203},{62,69,76,86,101}}
> (3 rows)
> 
> which by itself is kinda weird and that is why i unionized first, i.e.
> 
> 
> select (st_dumpvalues(st_union(rast))).* from dummy_rast;
> 
> and now get a more proper output:
> 
>  nband |valarray
> ---+-
> 
>  1 |
> {{253,254,253,254,254},{253,254,254,253,249},{250,254,254,252,249},{25
> 1,253,254,254,253},{252,250,254,254,254}}
>  2 |
> {{78,98,122,173,209},{96,118,180,249,254},{99,112,169,245,254},{89,99,1
> 22,176,229},{79,88,97,112,135}}
>  3 |
> {{70,86,100,135,161},{80,108,162,227,250},{90,108,175,251,254},{77,86,1
> 09,164,203},{62,69,76,86,101}}
> (3 rows)
> 
> however, my real concern is that when i compute the average over the 3
> bands:
> 
> select (st_dumpvalues(st_mapalgebra(array[row(rast, 1), row(rast, 2),
> row(rast, 3)]::rastbandarg[], 'st_mean4ma(double precision[], int[],
> text[])'::regprocedure))).* from (select st_union(rast) as rast from
> dummy_rast) as foo;
> 
> 
> i get
> 
> 
>  nband |valarray
> ---+-
> 
>  1 |
> {{133,146,158,187,208},{143,160,198,243,251},{146,158,199,249,252},{13
> 9,146,161,198,228},{131,135,142,150,163}}
> (1 row)
> 
> i.e. the actual means but rounded down, similarly for the standard
> deviation callback function.
> is this a bug because it shouldn’t be rounded down, we should have the
> decimals too.
> 
> thanx for feedback
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Create raster from row-col table

[Pierre Racine]

Étienne,

The columnx and rowy ST_SetValues() parameters are not references to each pixel 
values. They are the coordinates only of the first upperleft pixel to be set. 
Other pixels coordinates are automatically derived from the size of the passed 
two dimensional array (newvalueset). That means you must first aggregate your 
values into a two dimensional array. Try this:

-- Create an array of array aggregator because array_agg does not support 
arrays...
CREATE AGGREGATE array_agg_mult(anyarray) (
SFUNC = array_cat,
STYPE = anyarray,
INITCOND = '{}'
);

-- Aggregate the test.measure values into a two dimensional array and build the 
raster
WITH rows AS (
  SELECT colx, array_agg(v ORDER BY rowy) a
  FROM test.measure
  WHERE band = 1
  GROUP BY colx
), valarray AS (
  SELECT array_agg_mult(ARRAY[a] ORDER BY colx) a
  FROM rows
), raster AS (
  SELECT rast FROM test.rasters WHERE rid=1
), newraster AS (
SELECT ST_SetValues(rast, 1, 1, 1, valarray.a) rast
FROM valarray, raster
)
SELECT (ST_PixelAsPolygons(rast)).*
FROM newraster;

You might have to play a bit with the col and row ordering to get exactly what 
you want.

Pierre

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of Etienne B. Racine
> Sent: Friday, February 26, 2016 8:34 AM
> To: PostGIS Users Discussion
> Subject: [postgis-users] Create raster from row-col table
> 
> Hi list,
> 
> 
> I have raster data stored in a table row-col format as measure(rowy integer,
> colx integer, band integer, value float).
> 
> 
> I've tried loading a multidimensional array (row x col x bands) with no
> success, however building a 2d array (row x col) works but it seems
> contrived, especially once I need to stack the bands. I'm still figuring the
> best way to store the rasters, but I might end up storing 39 bands in a
> single raster. For now building an array isn't very fast, so I'm looking for
> speed improvement as well.
> 
> I wonder what's the most efficient way to populate a raster with data of
> this shape ? My actual method wouldn't work with sparse data and would
> require to pad the array with nodata values.
> 
> 
> I've put out a gist that explains my solution with sample data :
> https://gist.github.com/etiennebr/370a773029160c30d165#file-rowcol-
> raster-sql
> 
> 
> Thanks for your help,
> 
> Etienne
> 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] loading raster colour heatmap data or converting to a shape file

[Pierre Racine]

This all depends on the kind of operations you want to do... Some might be 
easier keeping rasters as raster. Others are easier if you convert everything 
to vector.

You can easily convert raster loaded in PostGIS to vector using 
ST_DumpAsPolygons(). The reverse operation is done through a mix of 
ST_AsRaster() and ST_Union().

Pierre

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of dave.po...@pinan.co.uk
> Sent: Friday, February 05, 2016 4:13 AM
> To: postgis-users@lists.osgeo.org
> Subject: [postgis-users] loading raster colour heatmap data or converting to
> a shape file
> 
> Hi list,
> 
> I have some heat map files which I want to do some spatial analyisis on.
> 
> I am wondering what the best way to proceed.  I admin that I have never
> used the raster part of postgis before.
> 
> I have two choices
> 
> 1. Convert the raster files to shape files and load the shape files in to
> postgis.
> 2. Load the raster files and have a play with some of the raster operators.
> 
> Has anybody got any suggestions as to the best way to proceed?
> 
> The heatmap are geotiff files and are correctly geoferenced and can be any
> colour that I choose.
> 
> As a side issue are there any colour combinations that  are easier to convert
> form raster to vector for  example would a black/white image digitise better
> than an image coloured in many different  shades of orange?
> 
> regards
> 
> Dave.
> 
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Storage of tile offset on a raster dataset

[Pierre Racine]

James,

> When importing a large raster dataset, it is typically tiled and the only
> reference in postgis that these tiles form a single dataset is the
> raster_columns table.

You could also load your rasters with the -f option which add a column with the 
filename of the original raster file. You would then have a quick way to 
identify tiles being part of a specific raster.

> I'd like postgis to keep track of the tile offsets (from upper left to lower
> right) as I believe this should bring about performance improvements by
> reducing geometric operations to simple scaling laws. The tile offset could
> be stored as a linear index for simplicity.
> 
> Consider trying to locate which tile in a dataset a point lies within. You
> can calculate the tile offset in x and y with the following scaling law:

Computing the relation between a simple point and a tile might be easier like 
this but PostGIS is conceived to solve much harder problem like: Is this line 
or this polygon inside the tile? Do they touch, intersects or are completely 
contained by the tile? For that you need more complex algorithms. In your 
simple case, what if the tiles are rotated?

The usual query to find whether a point is inside a tile (and the raster value 
matching this point) goes like this:

SELECT ST_Value(rast, ST_SetSRID(St_MakePoint(170.1420703, -43.594937), 4326))
FROM raster_table
WHERE ST_Intersects(rast, ST_SetSRID(St_MakePoint(170.1420703, -43.594937), 
4326))

I'm not sure the math under ST_Intersects() and ST_Value() is much more complex 
than what you are proposing. And for sure ST_Intersects() is much more 
polyvalent. What if you have 100 tiles? Does your maths benefit from a 
spacial index? ST_Intersects() does... when there is one defined.

>- Checking alignment of tiles. This is especially important when loading
> multiple datasets which should be aligned. Lets say I load 1 'base' dataset
> which gets tiled on the way in. Any subsequent tiles I load (e.g. to update
> a region of the dataset) need to be aligned. I can check that my input tile
> has the same corner coordinates of the tile it is replacing by calculating
> the tile index and pulling out the relevant information.

You can already compare raster alignment with ST_SameAlignment() and realign a 
raster with ST_Resample(), ST_Rescale(), ST_Reskew() and ST_SnapToGrid().

> - Iterating a raster dataset. Especially important when the processing
> you are doing is order-dependent (such as sequential neighbourhood
> searches
> (for e.g. terrain processing,  or imposing boundary conditions)

To find the neighbours of a tile you can intersect them with a small buffer 
around the tile of interest.

To return them always in the same sequence you can SORT them BY 
ST_UpperleftX(), ST_UpperleftY().

> To my mind, this seems quicker than the 'normal' method of doing
> something

You will have to demonstrate this affirmation... I think a lot of people would 
be amazed if you get something quicker than the actual PostGIS set of tools.

> Are there any facilities in postgis to support this already? Or am I
> overlooking some concepts about postgis that mean the above is
> unnecessary
> (i.e. the use cases can be met with no performance hit)

I think you underestimate the power of the PostgreSQL GiST index...

Good luck!

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Geoprocessing & BigData

[Pierre Racine]

> @Vincent : For ArcGIS Union, please see here.
> http://resources.esri.com/help/9.3/arcgisengine/java/gp_toolref/analysis_
> tools/union_analysis_.htm

ArcGIS Union resume to removing overlaps from the union of two tables. Look at 
the ST_DifferenceAgg() example in the PostGIS Addons:

https://github.com/pedrogit/postgisaddons/blob/master/postgis_addons.sql

This  produce a clean table (layer) without overlaps from the union of two 
(polygon) tables (layers).

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Help with SQL query?

[Pierre Racine]

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc...

Don't expect our suggestions to work blindly. I did not test this query. I'm 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said "your query should "look 
like" this"...

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I'm afraid the penny has not dropped yet :(

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure - albeit the results do not 
always please me :)

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of the data I have and I am struggling to figure out how best to 
pursue what I need to do.

I have a significant number of rasters which have double precision values.  
Without going into detail about what the rasters represent, I need to extract 
and sum values from one set of rasters in say table A based upon  values in 
another set of rasters in say table B  where the pixel value in the raster from 
Table B exceeds a threshold. Both tables are the same size (rasters are tiled) 
but I also need to figure out how I make sure the correct rasters are compared. 
 They have filenames like this rastervariable_10.tif, rastervariable_100.tif , 
presumably I need to use a logical expression to strip the numerical value (in 
this case this represents the year) and then order on that basis?

I can do this in QGIS one at a time but that is a little clumsy and rather time 
consuming.

If someone can just point me in the right dire

Re: [postgis-users] Help with SQL query?

[Pierre Racine]

Good catch Roxanne! Thanks!

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Roxanne Reid-Bennett
Sent: Monday, November 23, 2015 7:46 PM
To: postgis-users@lists.osgeo.org
Subject: Re: [postgis-users] Help with SQL query?

On 11/23/2015 12:41 PM, Darrel Maddy wrote:
Dear Pierre,

I was not looking for a total solution and I am grateful for the suggestion.

Although it may not look like it, I did consult the documentation and also I 
have Regina’s book beside me.  Unfortunately, for me at least, both documents 
assume some knowledge of SQL – which I do not have.  I am also trying to do 
this simultaneously with a large number of other things that are new to me.  I 
do not find the errors reported at all informative and consider the query I am 
trying to perform to be relatively trivial and hence I had hoped the structure 
of the query might have been more intuitive.  For others it may be.

FWIW - I don't play with rasters, but this appears to be a pure SQL thing... 
add "as rast" like below and try again.


SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated > 6 
THEN deposition ELSE NULL ENDIF ' )
as rast
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast)

Roxanne



I will try not to bother you again.

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 20:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org><mailto:postgis-users@lists.osgeo.org>
Subject: Re: [postgis-users] Help with SQL query?

The expression has to stay as it was: 'IF [rast2] > threshold THEN [rast1] ELSE 
NULL ENDIF '

Just replace the threshold value as you did.

Do not try to replace the [rast2] and [rast1]. They refer to the first and 
second raster pixel values. Read the ST_Mapalgebra doc…

Don’t expect our suggestions to work blindly. I did not test this query. I’m 
not in your context. I expect you read the doc about all the mentioned 
functions and adjust for your specific context. I said “your query should “look 
like” this”…

Pierre


From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 2:47 PM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Help with SQL query?

OK I spoke too soon.

I tried this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FROM (SELECT ST_MapAlgebra(deposition.rast, concentrated.rast, 'IF concentrated 
> 6 THEN deposition ELSE NULL ENDIF ' )
FROM mymodel.deposition, mymodel.concentrated
WHERE ST_UpperleftX(mymodel.deposition.rast) = 
ST_UpperleftX(mymodel.concentrated.rast) AND
 ST_UpperleftY(mymodel.deposition.rast) = 
ST_UpperleftY(mymodel.deposition.rast) ) foo;

And all I get is rast does not exist.

I’m afraid the penny has not dropped yet ☹

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: 23 November 2015 16:30
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Dear Pierre and Rasmus,

Many thanks for trying to help.

Rasmus: I am aware of GMT but I was looking for a solution in postgis so that I 
can keep all of the data extraction in one place.

Pierre: That is exactly what I was looking for and very many thanks for 
including the explanation.  I am a little overwhelmed with the number of 
functions offered in postgis. It is certainly a remarkable tool.  Watching the 
queries plough through my datasets is a pleasure – albeit the results do not 
always please me ☺

Hopefully I can put this to work later tonight.

Best wishes

Darrel




From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Pierre Racine
Sent: 23 November 2015 16:17
To: PostGIS Users Discussion 
<postgis-users@lists.osgeo.org<mailto:postgis-users@lists.osgeo.org>>
Subject: Re: [postgis-users] Help with SQL query?

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_Upper

Re: [postgis-users] Help with SQL query?

[Pierre Racine]

Darrel,


1)  Create a new raster selecting the right pixels with 
ST_MapAlgebra(raster, raster)

2)  Make sure only intersecting rasters are processed by using their upper 
left corner X and Y coordinates (with ST_UpperLeftX() and UpperLeftY())

3)  Sum the selected pixels with ST_SummaryStats(rast)

All in all a global query should look like this:

SELECT (ST_SummaryStats(ST_Union(rast))).sum sum
FFOM (SELECT ST_MapAlgebra(tableA.rast, tableB.rast, 'IF [rast2] > threshold 
THEN [rast1] ELSE NULL ENDIF ' )
FROM tableA, tableB
WHERE ST_UpperleftX(tableA.rast) = ST_UpperleftX(tableB.rast) AND
 ST_UpperleftY(tableA.rast) = 
ST_UpperleftY(tableB.rast) AND
 maybe some other condition here if you get time series 
e.g. tableA.year = tableB.year AND tableA.month = tableB.month) foo

If you have millions of tile you could create indexes on 
ST_UpperleftX(tableA.rast), ST_UpperleftX(tableB.rast), 
ST_UpperleftY(tableA.rast) and ST_UpperleftY(tableB.rast) to make the query 
faster.

You could also just use WHERE ST_Intersects(tableA.rast, tableB.rast) instead...

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
Darrel Maddy
Sent: Monday, November 23, 2015 7:20 AM
To: PostGIS Users Discussion
Subject: [postgis-users] Help with SQL query?

Dear all,

As you know I am relatively new to postgis and SQL and therefore  I have much 
to learn. However, I am facing a paper deadline and need to do some quick 
analysis of the data I have and I am struggling to figure out how best to 
pursue what I need to do.

I have a significant number of rasters which have double precision values.  
Without going into detail about what the rasters represent, I need to extract 
and sum values from one set of rasters in say table A based upon  values in 
another set of rasters in say table B  where the pixel value in the raster from 
Table B exceeds a threshold. Both tables are the same size (rasters are tiled) 
but I also need to figure out how I make sure the correct rasters are compared. 
 They have filenames like this rastervariable_10.tif, rastervariable_100.tif , 
presumably I need to use a logical expression to strip the numerical value (in 
this case this represents the year) and then order on that basis?

I can do this in QGIS one at a time but that is a little clumsy and rather time 
consuming.

If someone can just point me in the right direction I am sure I can figure out 
the rest for myself.

Apologies once more for asking what is probably a rather trivial question and 
yet again demonstrating my ignorance.

Many thanks

Darrel
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] Raster calculation on 3d raster (stack)

[Pierre Racine]

Doesn’t the second example (One raster, several bands) similar to your case in 
the doc?

http://postgis.net/docs/RT_ST_MapAlgebra.html

Pierre

From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of 
David Haynes
Sent: Tuesday, November 17, 2015 9:48 AM
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] Raster calculation on 3d raster (stack)

Hello,

I have a mulibanded raster and I am wondering if I can apply any of the raster 
processing functions (min, max, mean, stdev) to multibanded rasters. 
Specifically, I am wondering if these functions will be able to applied to the 
z-axis.

Here is an example, I have worked up in python starting with two 2d arrays.
x = [[1,4,5,6], [3,5,6,2]]
z = [[2,7,8,9], [8,3,5,9]]
convert to np.array
y = np.dstack((np.array(x), np.array(z)))
>>> y.sum(2)
array([[ 3, 11, 13, 15],
   [11,  8, 11, 11]])

I want to apply the raster processing functions to the z-axis and I am 
wondering if I can do it in PostGIS? It seems like the mapalgebra call back 
functions could do this, but I am having difficulty finding an example.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/postgis-users

Re: [postgis-users] SUM on rasters

[Pierre Racine]

If by "sum" you mean "merge" or "union" use ST_Union(). Be careful as the "sum" 
can exceed available memory if they are too many rasters in the table.

> -Original Message-
> From: postgis-users [mailto:postgis-users-boun...@lists.osgeo.org] On
> Behalf Of Stephen Crawford
> Sent: Monday, October 12, 2015 10:46 AM
> To: PostGIS Users Discussion
> Subject: [postgis-users] SUM on rasters
> 
> Hi,
> 
> I have a table of rasters.  I would like to SUM the entire table with
> one raster as the result. Can somebody point me in to the correct
> funtion to do this?
> 
> Thanks,
> Steve
> 
> 
> ___
> postgis-users mailing list
> postgis-users@lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster2pgsql wildcards on windows

[Pierre Racine]

Copy all your files in a single folder and only use one wildcard.

d:\postgis_in_action\ch13_data_code\vietnam\dted\allfiles\*.dt0

From: postgis-users-boun...@lists.osgeo.org 
[mailto:postgis-users-boun...@lists.osgeo.org] On Behalf Of toni hernández
Sent: Tuesday, October 06, 2015 10:49 AM
To: PostGIS Users Discussion
Subject: [postgis-users] raster2pgsql wildcards on windows

Hi all,

How can I import several raster files from several folders with raster2pgsql on 
windows?

This is what I am trying unsuccessfully.

raster2pgsql -r d:\postgis_in_action\ch13_data_code\vietnam\dted\*\*.dt0 -s 
4326 -F -t 50x50 -I vietelev > vietelev.sql

"Unable to read raster file"  is what I get.
--
Toni Hernández Vallès
Servei de Sistemes d'Informació Geogràfica i Teledetecció
-
Universitat de Girona
SIGTE
-
Pl. Ferrater Mora 1
17071 Girona
Tel +34 972 418 039 (7026 intern)
t...@sigte.udg.edu

http://www.sigte.udg.edu
Twitter http://twitter.com/SIGTE_UDG
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_Clip

[Pierre Racine]

I guess the band argument comes before the geom argument in ST_Clip(). You can 
pass an array for the bands you require.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of David Haynes
 Sent: Wednesday, July 22, 2015 4:40 PM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] ST_Clip
 
 Hello
 
 I believe I am in appropriately using the ST_Clip function on a multiband
 raster. I have a raster with 13 bands and I want to retrieve from the clip
 bands # 1,3
 
 How do I do that?
 
 I tried using:
 array[1,3] as bnd
 
 Didn't work still returns 13 bands.
 
 with bnd_num as
 (
 select array[1,3] as bnd
 )
 select p.geoid, p.label, ST_Numbands(ST_CLIP(r.rast,p.geom,b.bnd, True))
 from bnd_num b, polygon p inner join modis_igbp_stack r on
 ST_Intersects(r.rast, p.geom);
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster, stats conditioned to a set of values

[Pierre Racine]

 Adding a GROUP BY t1.rid, doesn't solve the problem, still get an error.

Which error?
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster, stats conditioned to a set of values

[Pierre Racine]

These are not error. Only notices. Are you sure the query returned nothing?

Otherwise I would restart with something more simple to double check some 
prerequired things. 

1) Make sure your two raster table actually contains pixels with values

2) Make sure some tiles actually intersects

Do you have a way to display the footprints of your rasters tiles?

Pierre

 -Original Message-
 From: juli g. pausas [mailto:juli.g.pau...@uv.es]
 Sent: Monday, July 06, 2015 9:58 AM
 To: Pierre Racine
 Cc: PostGIS Users Discussion
 Subject: Re: [postgis-users] raster, stats conditioned to a set of values
 
 
 INSERT INTO rastertmp.provanet
 SELECT t1.rid, ST_Union(ST_Intersection(t1.rast, 1, t2.rast, 1, 'BAND1'))
 FROM rastertmp.prova AS t1, rastertmp.provapositiu AS t2
 WHERE ST_Intersects(t1.rast, t2.rast)
 GROUP BY t1.rid;
 
 
 NOTICE:  The two rasters provided have no intersection.  Returning no band
 raster
 CONTEXT:  PL/pgSQL function
 st_intersection(raster,integer,raster,integer,text,double precision[]) line 20
 at assignment
 NOTICE:  Could not find raster band of index 1 when setting pixel value.
 Nodata value not set. Returning original raster
 
 
 
 
 Juli
 --
 CIDE, CSIC  |  www.uv.es/jgpausas  |
 
 
 
 On Mon, Jul 6, 2015 at 3:52 PM, Pierre Racine
 pierre.rac...@sbf.ulaval.ca wrote:
 
 
Adding a GROUP BY t1.rid, doesn't solve the problem, still get an
 error.
 
   Which error?
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster, stats conditioned to a set of values

[Pierre Racine]

 But my main problem is that I would like to do this (e.g., the query above),
 but only for the pixels in which Band2 = 0.  Any idea? any clue?

1) isolate the pixels with 0. For this you can use ST_MapAlgebra() or 
ST_Reclass(). I use ST_MapAlgebra() in the following query
2) compute the intersection with band 1
3) compute the stats as you did

SELECT region_cod, (res).* 
FROM 
  (SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1, p.geom, true)) AS res
FROM gis_wd.wd_regiones AS p, rastertmp.ndvitmp AS r
WHERE ST_Intersects(r.rast, p.geom) 
AND p.region_cod = 'PA1214'
   ) AS foo WHERE (res).value   0;


WITH bands AS ( -- reclass the second band to 0 and 1
SELECT ST_MapAlgebra(ST_Band(rast, 2), '16BSI'::text, 'CASE WHEN [rast]  0 or 
[rast]  0 THEN NULL ELSE 1 END') band2,
 ST_Band(rast, 1) band1
) rastintersect AS ( -- compute the intersection of band 1 and band 2
  SELECT ST_Intersection(band1, band2, 'BAND1') rast FROM bands
)
SELECT region_cod, (res).* 
FROM 
  (SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1, p.geom, true)) AS res
FROM gis_wd.wd_regiones AS p, rastintersect AS r
WHERE ST_Intersects(r.rast, p.geom) 
AND p.region_cod = 'PA1214'
   ) AS foo;

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Problem in Clip between geometry and raster table

[Pierre Racine]

Why is the first argument to ST_Clip() is a geometry? Should be a raster.

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Marcello Benigno
 Sent: Wednesday, July 01, 2015 7:12 AM
 To: PostGIS Users Discussion
 Subject: [postgis-users] Problem in Clip between geometry and raster table
 
 Hello guys,
 
 When attempting to perform a clip from a polygon and a raster table , the
 following error appeared:
 
 
   ERROR:  column nan does not exist
   LINE 1: SELECT (NaN)::double precision
   ^
   QUERY:  SELECT (NaN)::double precision
   CONTEXT:  PL/pgSQL function
 st_clip(raster,integer,geometry,double precision[],boolean) line 42 at
 assignment
   SQL function st_clip statement 1
   (execution time: 250 ms; total time: 485 ms)
SELECT (ST_DumpAsPolygons(ST_Clip(c.geom,
 ST_Buffer(f.geometria, 0.01), TRUE))).val AS dn,
 (ST_DumpAsPolygons(ST_Clip(c.geom, ST_Buffer(f.geometria,
 0.01), TRUE))).geom AS geom
   FROM public.propriedade_geometria f, dados.edr c
   WHERE  ST_Intersects(f.geometria, c.geom)
   AND f.tipo = 2 AND f.propriedade_id = 6063 AND c.id_uf = 2 AND
 c.cidade_id = 4681;
 
 
 I can't understand what may be happening. Can I perform a casting to solve
 this problem ?
 
 --
 
 Marcello Benigno B. de Barros Filho
 Prof. do Curso Superior de Tecnologia em Geoprocessamento - IFPB
 Mestre em Ciências Geodésicas e Tecnologias da Geoinformação - UFPE
 Doutorando em Tecnologia Ambiental e Recursos Hídricos - UFPE
 http://profmarcello.blogspot.com
 http://about.me/marcello.benigno http://about.me/marcello.benigno
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster, stats conditioned to a set of values

[Pierre Racine]

Right! Actually you need a more complex expression returning not a boolean but 
just the right value:

'CASE WHEN [rast]  = 0 THEN 0 ELSE [rast] END'

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of juli g. pausas
 Sent: Wednesday, July 01, 2015 5:34 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] raster, stats conditioned to a set of values
 
 Thanks for this.
 
 But it didn't work for me.
 
 
 These queries work
 
 
 SELECT (ST_SummaryStats(ST_MapAlgebra(rast, 1, '16BSI'::text, '[rast]+1'),
 1)).* FROM rastertmp.ndvitmp
 SELECT (ST_SummaryStats(ST_MapAlgebra(rast, 1, NULL, '[rast]+1'), 1)).*
 FROM rastertmp.ndvitmp
 
 
 But not when I use the condition [rast]  0
 
 SELECT (ST_SummaryStats(ST_MapAlgebra(rast, 1, NULL, '[rast]  0'), 1)).*
 FROM rastertmp.ndvitmp
 
 ERROR:  cannot cast type boolean to double precision
 LINE 1: SELECT ($1  0)::double precision
 
 
 Replacing NULL for '16BSI'::text, or '16BSI'::smallint, '16BSI'::double
 precision doen't solve the problem
 
 
 Any other suggestion?
 
 Thanks
 
 
 
 
 
 
 
 Juli
 --
 CIDE, CSIC  |  www.uv.es/jgpausas  |
 
 
 
 On Tue, Jun 30, 2015 at 5:29 PM, Pierre Racine
 pierre.rac...@sbf.ulaval.ca wrote:
 
 
   You can select pixels fulfilling an expression using the one-raster
 variant of ST_MapAlgebra
 
   http://postgis.net/docs/RT_ST_MapAlgebra_expr.html
 
   So just do something like:
 
   ST_MapAlgebra(rast, '8BUI'::text, '[rast]  0')
 
   before computing stats.
 
   To count the number of pixels of a certain value you can use
 ST_ValueCount().
 
   Pierre
 
 
-Original Message-
From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-
 users-
boun...@lists.osgeo.org] On Behalf Of juli g. pausas
Sent: Monday, June 29, 2015 12:49 PM
To: PostGIS Users Discussion
Subject: [postgis-users] raster, stats conditioned to a set of values
   
Hi all
   
I'm just starting to discover postgis, it is really useful.
   
I have a raster file with different bands, that I have imported to
 postgres
(raster2pgsql, without the -R option, i.e., insite the database). I
 would like to
extract information from band 1, but filtering the data using only
   
1) pixels with positive values (in that band, band 1)
   
2) and pixels in which band 2 is equal to a given value, e.g. 1
   
   
The type of analysis I'd like to do is descriptive stats but also
 intersecting
with a vector map. If my raster is rastertmp.ndvitmp, two
 examples are:
   
   
SELECT (ST_SummaryStats(rast, 1)).* FROM rastertmp.ndvitmp
 WHERE rid
= 1
   
   
   
SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1, p.geom,
 true)) AS res
  FROM gis_wd.wd_regiones AS p, rastertmp.ndvitmp AS r
  WHERE ST_Intersects(r.rast, p.geom)
  AND p.region_cod = 'PA1214';
   
   
   
This works perfectly, but how can I compute the stats  only for
 pixels with
positive values and with a given value in another band?  The idea
 would be
something like:  WHERE ST_Values(rast, 1)0 AND ST_Values(rast,
 2)=1
   
   
Thanks for any help!
   
   
Juli
--
CIDE, CSIC  |  www.uv.es/jgpausas  |
   
 
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster, stats conditioned to a set of values

[Pierre Racine]

You can select pixels fulfilling an expression using the one-raster variant of 
ST_MapAlgebra

http://postgis.net/docs/RT_ST_MapAlgebra_expr.html

So just do something like:

ST_MapAlgebra(rast, '8BUI'::text, '[rast]  0')

before computing stats.

To count the number of pixels of a certain value you can use ST_ValueCount().

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of juli g. pausas
 Sent: Monday, June 29, 2015 12:49 PM
 To: PostGIS Users Discussion
 Subject: [postgis-users] raster, stats conditioned to a set of values
 
 Hi all
 
 I'm just starting to discover postgis, it is really useful.
 
 I have a raster file with different bands, that I have imported to postgres
 (raster2pgsql, without the -R option, i.e., insite the database). I would 
 like to
 extract information from band 1, but filtering the data using only
 
 1) pixels with positive values (in that band, band 1)
 
 2) and pixels in which band 2 is equal to a given value, e.g. 1
 
 
 The type of analysis I'd like to do is descriptive stats but also intersecting
 with a vector map. If my raster is rastertmp.ndvitmp, two examples are:
 
 
 SELECT (ST_SummaryStats(rast, 1)).* FROM rastertmp.ndvitmp WHERE rid
 = 1
 
 
 
 SELECT p.region_cod, ST_ValueCount(ST_Clip(r.rast,1, p.geom, true)) AS res
   FROM gis_wd.wd_regiones AS p, rastertmp.ndvitmp AS r
   WHERE ST_Intersects(r.rast, p.geom)
   AND p.region_cod = 'PA1214';
 
 
 
 This works perfectly, but how can I compute the stats  only for pixels with
 positive values and with a given value in another band?  The idea would be
 something like:  WHERE ST_Values(rast, 1)0 AND ST_Values(rast, 2)=1
 
 
 Thanks for any help!
 
 
 Juli
 --
 CIDE, CSIC  |  www.uv.es/jgpausas  |
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Using Postgis raster for water elevation

[Pierre Racine]

Or ST_Reclass() reclassifying everything  10 to nodata...

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rémi Cura
 Sent: Wednesday, May 27, 2015 9:31 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] Using Postgis raster for water elevation
 
 Hey,
 
 
 Map Algebra?
 
 (custom function with a thresholding?)
 
 
 
 Cheers,
 
 Rémi-C
 
 
 2015-05-27 14:57 GMT+02:00 Gabriel Vatin gabriel.va...@kinaxia.fr:
 
 
   Hello all,
 
   I'd like to have some feedbacks on an operation that should be
 possible with Postgis raster, but I can't find out how to do this.
   I have a whole raster of French DEM (elevation data) stored on
 PostGis, with tiled objects. That makes a table with some 14.000 rows of
 data: france_mnt (rid, rast, filename)
   With ST_Value(france_mnt.rast, geometry), I can retrieve the
 elevation of a given coordinates in the geometry object.
 
   I'd like to create a new raster, on the fly, with a maximum value of
 elevation, such as : all pixels with value  10
   This could be used for mapping water elevation in a area, for
 instance.
 
   Anyone has already done this? I'll be happy to have some tips on
 the function to use.
   Thanks!
 
   Gabriel
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] What is the Polygon thinning function in PostGIS?

[Pierre Racine]

Read this:

https://smathermather.wordpress.com/2012/08/10/what-is-the-center-line-of-a-complex-polygon-routing-stream-and-rivers-part-two/

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Tom van Tilburg
 Sent: Wednesday, May 20, 2015 6:07 AM
 To: postgis-users@lists.osgeo.org
 Subject: Re: [postgis-users] What is the Polygon thinning function in
 PostGIS?
 
 Hi,
 
 A straight skeleton might be what you need:
 http://postgis.net/docs/ST_StraightSkeleton.html
 
 I know there are some blog posts on this but couldn't find them that
 quickly.
 
 Best,
 Tom
 
 On 20-5-2015 4:29, KhunSanAung wrote:
 
 
   Hello All,
 
   I'd like to generated line from elongated polygon in PostGIS.
   e.g. extracting center of a river polygon, center of the road polygon,
 etc.
 
   I'd to know what would be equivalence of Polygon Thinning in
 PostGIS.
 
   Thank you very much for any hints.
 
   Best regards
   --
 
   Have a nice day!
   --
 
 
   Mr. Khun San Aung
 
 
 
 
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] create table from existing raster table

[Pierre Racine]

What do you mean by proper? How is it? Your query should results in a raster 
column.

Pierre

From: postgis-users-boun...@lists.osgeo.org 
[postgis-users-boun...@lists.osgeo.org] On Behalf Of Ahmet Temiz 
[ahmettemi...@gmail.com]
Sent: Friday, May 08, 2015 8:33 AM
To: PostGIS Users Discussion; PostGIS Users Discussion
Subject: [postgis-users] create table from existing raster table

hello

I was wondering how we could create table from existing raster table ?

I tried:
create table newrastertable as select ST_Slope(rast) from existingDEMrastertable

​but I couldn't get proper raster table

how can we do that ?

regards​

--
Ahmet Temiz
Jeoloji Müh.
Afet ve Acil Durum Yönetimi Başkanlığı
Bilgi İşlem  Dairesi Başkanlığı-CBS Grubu




Ahmet Temiz
Geological Eng.
Information Systems - GIS Group
Disaster and Emergency Management
of Presidency
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Using ST_Slope in EPSG:4326

[Pierre Racine]

You have to project your raster wit ST_Transform() before calling ST_Slope().

You also have to ST_Union() your raster tiles before projecting it.

Pierre

From: postgis-users-boun...@lists.osgeo.org 
[postgis-users-boun...@lists.osgeo.org] On Behalf Of Ahmet Temiz 
[ahmettemi...@gmail.com]
Sent: Friday, May 08, 2015 5:00 AM
To: PostGIS Users Discussion; PostGIS Users Discussion
Subject: [postgis-users] Using ST_Slope in EPSG:4326

hello,

( sorry for previous untitled message)

How can we get correct values from ST_Slope, if raster table  is in EPSG:4326 ?
or
we have to transform dem table to metric projection.

please let me know

kind regards

--
Ahmet Temiz
Jeoloji Müh.
Afet ve Acil Durum Yönetimi Başkanlığı
Bilgi İşlem  Dairesi Başkanlığı-CBS Grubu




Ahmet Temiz
Geological Eng.
Information Systems - GIS Group
Disaster and Emergency Management
of Presidency
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Computation of quantiles in Polygon-Raster

[Pierre Racine]

I guess you have two options:

1) clip and union your tiles BEFORE computing the quantile

2) use the ST_Quantile variant taking table and column names. I don't know if 
these function aggregates the results though.

Pierre 

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Nick Ves
 Sent: Monday, May 04, 2015 3:53 PM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] Computation of quantiles in Polygon-Raster
 
 Hi list,
 
 I am trying to calculate the quantiles of a (tiled) raster which is
 bound by a polygon but I am stumbled.  I can calculate a desired
 quantile over a rid using ST_Quantile like so:
 
 SELECT b.rid, q.*
 FROM raster.polygon a, raster.sat1 b,
 LATERAL
 st_quantile(st_clip(b.rast,a.geom),ARRAY[0.25,0.75]) q ORDER BY quantile;
 
 (sample result)
 
 rid quantile value
 11 0,25 90
 78 0,25 90
 40 0,25 78
 79 0,25 74
 21 0,25 104
 80 0,25 74
 41 0,25 57
 
 13 0,75 152
 14 0,75 98
 15 0,75 93
 16 0,75 98
 17 0,75 93
 18 0,75 0
 19 0,75 172
 20 0,75 177
 ...
 
 but i don't know how to aggregate the results to create a
 representative result for the whole dataset.
 
 Anyone has any hint for me on how to proceed ?
 
 N
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Calculate area of raster/polygon intersections in PostGIS

[Pierre Racine]

 1 - What's the best way to do raster/polygon intersections? I read about
 two different approaches, one using ST_Intersection and the other one
 using ST_Clip. The second one apparently is faster and more efficient, but i
 don't know about its precision.

The precision depends on the sizes of your pixels comparing with the size of 
your buffers. If pixels are relatively small you should get a good precision. 
If pixels are big precision will be poor and in that case you could use the 
ST_Intersection method. But be careful at false precision... 

 2 - Assuming that i used ST_Clip to do the intersections, the result it's
 gonna be another raster. How can i calculate the area of each category of
 this raster?

Use ST_ValueCount() times the size of one pixel (make sure to be projected 
properly).

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Problem to visualize raster data from PostGIS inQGIS

[Pierre Racine]

Just

SELECT (ST_DumpAsPolygons(rast)).*
FROM deforested_areas
-–WHERE rid = 1

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Hugues François
 Sent: Tuesday, March 31, 2015 10:05 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] Problem to visualize raster data from PostGIS
 inQGIS
 
 Hello,
 
 
 
 To display raster data, you can try ST_DumpAsPolygons into Arcgis or
 OpenJump (if it is just a test, you could add a where clause on tile id).
 
 
 
 WITH a AS (select ST_DumpAsPolygons(rast) AS mydump from
 deforested_areas –where rid = 1)
 
 
 
 SELECT (mydump).val, (mydump).geom FROM a
 
 
 
 This website is a bit outdated, but you can find some useful insights :
 http://trac.osgeo.org/postgis/wiki/WKTRaster
 
 
 
 Hug
 
 
 
 
 
 De : postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] De la part de Roberto Lazarte Kaqui
 Envoyé : mardi 31 mars 2015 15:11
 À : postgis-users@lists.osgeo.org
 Objet : [postgis-users] Problem to visualize raster data from PostGIS inQGIS
 
 
 
 I have two multiband raster files and i'm trying to load these data into
 PostGIS using “raster2pgsql”.
 
 Here is the syntax that i used:
 
 
 
 raster2pgsql -c -s 4326 -I -t 100x100 *.tif  public.deforested_rasters | psql 
 -
 d r342471958
 
 
 
 It works without error, but i don't know if the command loads correctly
 both raster files. The database in PostGIS has two columns, one with the
 spatial index and the other one with raster information. I'm trying to
 visualize the resultant database in QGIS, but for some reason it doesn't
 open the raster layer.
 
 I've tried to solve this problem using this idea:
 https://duncanjg.wordpress.com/2012/11/20/the-basics-of-postgis-raster/
 https://duncanjg.wordpress.com/2012/11/20/the-basics-of-postgis-
 raster/ , transforming the raster data to geometry data through DB
 manager using this command:
 
 
 
 select rid,rast::geometry from deforested_rasters
 
 
 
 but when visualized in QGIS, only returns a rectangular area without the
 original characteristics.
 
 Any ideas of how can i confirm that the “raster2pgsql” command loads
 correctly the raster files? Or a different way to load these data into
 PostGIS?
 
 The raster files are available
 in:https://www.dropbox.com/sh/nu8nl6cigv8viz1/AADTrr_w8SCstDHOu6A
 iGeJta?dl=0
 https://www.dropbox.com/sh/nu8nl6cigv8viz1/AADTrr_w8SCstDHOu6Ai
 GeJta?dl=0
 
 
 
 Thanks,
 
 Roberto
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster2pgsql creates tables without values

[Pierre Racine]

How do you check that your rasters only nodata value?

You know, very often, raster and geometries are not displayed in the pgAdmin 
table result when they are too large.

What's the result of 

SELECT postgis_full_version()

on both machines?

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Andreas Laggner
 Sent: Monday, March 23, 2015 12:15 PM
 To: PostGis_Mailinglist
 Subject: [postgis-users] raster2pgsql creates tables without values
 
 Hi list,
 
 all my rastermaps i import are empty, no values (NaN only).
 
 Import on another machine with a different installation of PostGIS for
 this data works good.
 I also tried to split the import, first create the sql-file and than
 import with psql - NO ERRORS accurs, but no values.
 
 It must have to do with my installation of PostGIS or GDAL or something
 else. I don't know how to find/check this.
 
 I have no idea - perhaps someone else out there has
 
 hasta luegoAndreas
 
 
 
 --
 Dipl. Geoökologe Andreas Laggner
 
 Thünen-Institut für Agrarklimaschutz (AK)
 Arbeitsgruppe Emissionsinventare
 Johann Heinrich von Thünen-Institut (vTI),
 Bundesforschungsinstitut für Ländliche Räume, Wald und Fischerei
 
 Thünen Institute of Climate-Smart Agriculture
 Johann Heinrich von Thünen Institute (vTI),
 Federal Research Institute for Rural Areas, Forestry and Fisheries
 
 Bundesallee 50
 D-38116 Braunschweig
 
 Tel.: (+49) (0)531 596 2636
 Fax : (+49) (0)531 596 2645
 E-mail: andreas.lagg...@ti.bund.de
 Homepage: http://www.ti.bund.de
 
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster2pgsql: -t TILE_SIZE

[Pierre Racine]

 If you're storing the raster in the database, I have no specific suggestions.
 The reason for no suggestion is that it really depends on the characteristics
 of the raster (number of bands, pixel type of each band) and the server
 environment.

You normally want to store in-db because you want to do some kind of analysis. 
In that case, the optimal tile size is dependent on the kind of analyses you 
want to do. You tend to want to minimize the number of tiles loaded per 
analysis (bigger tiles) and at the same time minimize loading areas not 
involved in the analysis (smaller tiles). For example if you want to get an 
aggregated raster value per some kind of polygons, the process will be faster 
if the tiles size is about the mean size of your polygons.

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[postgis-users] ST_Clip on a raster without band crashes the server

[Pierre Racine]

Hi,

I have:

POSTGIS=2.1.5 r13152 GEOS=3.4.2-CAPI-1.8.2 r0 PROJ=Rel. 4.8.0, 6 March 
2012 GDAL=GDAL 1.11.1, released 2014/09/24 LIBXML=2.7.8 LIBJSON=UNKNOWN 
RASTER

and 

SELECT ST_Clip(ST_MakeEmptyRaster(42, 42, 0, 0, 1.0, 1.0, 0, 0, 4269), 
ST_MakeEnvelope(0, 0, 20, 20, 4269));

crashed my server.

Someone can confirm the bug?

Thanks,

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_Clip on a raster without band crashes the server

[Pierre Racine]

I think when the ST_Clip() crop parameter is set to TRUE, ST_Clip() could just 
change the extent of the empty raster When it is set to false, it just 
returns the raster unmodified.

I'll fill a bug.

Thanks,

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Bborie Park
 Sent: Tuesday, February 24, 2015 10:25 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] ST_Clip on a raster without band crashes the
 server
 
 File a bug. It shouldn't work for an empty raster.
 
 On Tue, Feb 24, 2015 at 6:50 AM, Pierre Racine
 pierre.rac...@sbf.ulaval.ca wrote:
 
 
   Hi,
 
   I have:
 
   POSTGIS=2.1.5 r13152 GEOS=3.4.2-CAPI-1.8.2 r0 PROJ=Rel.
 4.8.0, 6 March 2012 GDAL=GDAL 1.11.1, released 2014/09/24
 LIBXML=2.7.8 LIBJSON=UNKNOWN RASTER
 
   and
 
   SELECT ST_Clip(ST_MakeEmptyRaster(42, 42, 0, 0, 1.0, 1.0, 0, 0,
 4269), ST_MakeEnvelope(0, 0, 20, 20, 4269));
 
   crashed my server.
 
   Someone can confirm the bug?
 
   Thanks,
 
   Pierre
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Elevation profiles from rasters

[Pierre Racine]

Try 100x100...

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rasmus Aveskogh
 Sent: Wednesday, August 27, 2014 11:08 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] Elevation profiles from rasters
 
 Thanks,
 
 I'm using the technique described at http://blog.mathieu-
 leplatre.info/drape-lines-on-a-dem-with-postgis.html
 
 What would be small tiles in this context? 256x256? The raster I was testing
 against was ~1200x1200.
 
 -ra
 
 On 27 Aug 2014, at 16:19, Pierre Racine pierre.rac...@sbf.ulaval.ca
 wrote:
 
 
   Elevation profile should be quite fast depending on how you tile
 your raster coverage. Smaller tiles make computation faster. Index them
 and make sure your query is actually using the index.
 
   Pierre
 
 
 
   -Original Message-
   From: postgis-users-boun...@lists.osgeo.org
 [mailto:postgis-users-
   boun...@lists.osgeo.org] On Behalf Of Rasmus Aveskogh
   Sent: Friday, August 22, 2014 4:10 PM
   To: postgis-users@lists.osgeo.org
   Subject: [postgis-users] Elevation profiles from rasters
 
   Hi!
 
   How hard/feasible would it be to implement elevation
 profile/drape line
   functionality in PostGIS? Since PostGIS already depends on
 GDAL I assume
   making use of the GDAL API would make this quite easy,
 such as in the Zoo-
   Project service:
 
   http://www.zoo-project.org/trac/browser/trunk/zoo-
 project/zoo-
   services/gdal/profile/service.c
 
   Providing functionality that could eventually be used for
 applications like:
   http://www.zoo-
 project.org/site/ZooWebSite/Demo/GdalProfile (using the
   code above)
 
   I'm currently using an ad-hoc PostGIS-function to pull out
 elevation profiles
   from my rasters, but it's quite slow and CPU intensive and
 as such does not
   make a good fit as a backend service for a web application
 like the above.
 
   -ra
 
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Elevation profiles from rasters

[Pierre Racine]

Elevation profile should be quite fast depending on how you tile your raster 
coverage. Smaller tiles make computation faster. Index them and make sure your 
query is actually using the index.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rasmus Aveskogh
 Sent: Friday, August 22, 2014 4:10 PM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] Elevation profiles from rasters
 
 Hi!
 
 How hard/feasible would it be to implement elevation profile/drape line
 functionality in PostGIS? Since PostGIS already depends on GDAL I assume
 making use of the GDAL API would make this quite easy, such as in the Zoo-
 Project service:
 
 http://www.zoo-project.org/trac/browser/trunk/zoo-project/zoo-
 services/gdal/profile/service.c
 
 Providing functionality that could eventually be used for applications like:
 http://www.zoo-project.org/site/ZooWebSite/Demo/GdalProfile (using the
 code above)
 
 I'm currently using an ad-hoc PostGIS-function to pull out elevation profiles
 from my rasters, but it's quite slow and CPU intensive and as such does not
 make a good fit as a backend service for a web application like the above.
 
 -ra
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] PostGIS Raster Performance

[Pierre Racine]

You should use the raster stats functions instead of ST_DumpAsPolygon to get 
values or summary stats from the clipped rasters. 
http://postgis.net/docs/RT_reference.html#RasterBand_Stats

Pierre 

From: postgis-users-boun...@lists.osgeo.org 
[postgis-users-boun...@lists.osgeo.org] On Behalf Of DanielFranco 
[dfra...@computacao.ufla.br]
Sent: Wednesday, August 20, 2014 5:00 PM
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] PostGIS Raster Performance

Hi, I'm a newbie in PostGIS Raster use and GIS applications, so this may be a
dumb question. I'm a little confused by some functionalities and I don't
know if it's possible to do what i'm trying to. Basically, I'm using postgis
functions either for raster and vector forms of an image and I'm trying to
visualize parts of the raster by using the informations presented only in
the shapefile. For example, I have an TIFF file of a land coverage (one
band) about 49746 x 71819 and it's shapefile with some grid codes. I'm doing
the following query to extract pixel values from an expecific gridcode (18)
presented in the shapefile imported (vector_table) so I could infer some
grid codes by the pixel value:

SELECT (gval).val as pixel_value
FROM
( SELECT ST_DumpAsPolygons(ST_Clip(rast,1,geom)) as gval
from raster_table
join vector_table
on (st_intersects(rast,geom))
where geom in (select geom from vector_table where gridcode = 18)
) as intersection;

The query runs too slow and I don't know how to optimize it. The raster
image is big and I tried different tiles configurations (100x100, 200x200,
1000x1000) and 100x100 seemed to have better perfomance. I'm using st_clip
for a faster answer. I don't know if these pixels values would retrieve some
important information and if PostGIS raster is really suitable with area
calculations. Is there a way to optimize the query or postgis raster isn't
suitable for this case?



--
View this message in context: 
http://postgis.17.x6.nabble.com/PostGIS-Raster-Performance-tp5006902.html
Sent from the PostGIS - User mailing list archive at Nabble.com.
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Summary Statistics (vector, raster)

[Pierre Racine]

Run the whole postgis_addons.sql like a sql command in pgAdmin.

 -Original Message-
 From: Lucas Ferreira Mation [mailto:lucasmat...@gmail.com]
 Sent: Friday, June 13, 2014 2:48 PM
 To: postgis-us...@googlegroups.com
 Cc: postgis-us...@postgis.refractions.net; postgis-
 us...@postgis.refractions.net; Pierre Racine
 Subject: Re: [postgis-users] Summary Statistics (vector, raster)
 
 POSTGIS begginer here:
 
 I think I'm having a similar problem as described in this post. I tryed
 following Pierre's sucgestion (second post), but I'm not able to find or
 install the
 SummarystatsAgg
  and
 AreaWeightedSummaryStats
 functions.
 
 how do I do that?
 tks
 Lucas
 
 
 
 
 
 Em terça-feira, 6 de março de 2012 10h48min46s UTC-3, Pierre Racine
 escreveu:
 
Are the SummarystatsAgg and AreaWeightedSummaryStats
 included in the latest
trunk or should i install from plpgsql?
 
   That's what I forgot: They are not part of rtpostgis.sql... They are in
 scripts/plpgsql.
 
   Pierre
   ___
   postgis-users mailing list
   postgi...@postgis.refractions.net javascript:
   http://postgis.refractions.net/mailman/listinfo/postgis-users
 http://postgis.refractions.net/mailman/listinfo/postgis-users
 
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_SelectByValue

[Pierre Racine]

Try ST_Reclass().

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of David Haynes
 Sent: Friday, May 30, 2014 10:49 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] ST_SelectByValue
 
 Hello,
 
 I am fairly new to using rasters in postgis, but I am wondering if there is a
 faster way to extract values from categorical rasters (e.g. landcover types)
 than using map algebra. In the PostGIS in Action book, pg 404 has a
 function called ST_SelectByValue(raster|geometry, 'expression') which
 seems to incorporate map algebra, however I can't seem to tell if this
 function is finished and available in postgis as I get an unknown function
 error.
 
 My version is this.
 PostgreSQL 9.3.2 on x86_64-unknown-linux-gnu, compiled by gcc
 (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit POSTGIS=2.1.1 r12113
 GEOS=3.3.3-CAPI-1.7.4 PROJ=Rel. 4.7.1, 23 September 2009
 GDAL=GDAL 1.9.0, released 2011/12/29 LIBXML=2.7.8 LIBJSON=UN
 (...)
 
 
 And I found a webpage
 http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking03
 
 search for ST_SelectByValue at the bottom, indicating that it was finished.
 
 
 AAccomplished Objectives
 http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking03#A
 ccomplishedObjectives
 
 
 ST_ValueCount(raster, value) - integer - done see below
 
 ST_ValuePercent(raster, value) - double precision - done see below
 
 ST_Resample(raster, method, originx, originy, pixelsizex, pixelsizey) - 
 raster
 - done see below
 
 ST_SelectByValue(raster, 'expression') - same type as first argument
 
 ST_Clip(raster|geometry,geometry) - same type as first argument

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] raster support

[Pierre Racine]

These variables are system variables, not PostgreSQL variables. to know how to 
set them in your system search for set environment variable os x.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of John Payne
 Sent: Sunday, June 01, 2014 10:25 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] raster support
 
 I want to store raster data in PostGIS.  I am using Postgresql 9.2.4 with
 PostGIS 9.1.0 on a Mac Powerbook running OS X 10.8.5.  I want to use the
 command raster2pgsql but I get an error message saying that the type
 raster was not recognized.
 
 I think the problem may be (from Chapter 2.1 of the PostGIS manual):
 
 
 As of PostGIS 2.1.3, out-of-db rasters and all raster drivers are disabled by
 default. In order to re-enable these, you need to set the following
 environment variables:POSTGIS_GDAL_ENABLED_DRIVERS and
 POSTGIS_ENABLE_OUTDB_RASTERS in the server environment
 
 
 I tried adding the following to the end of my postgresql.conf file:
 POSTGIS_GDAL_ENABLED_DRIVERS=ENABLE_ALL
 POSTGIS_ENABLE_OUTDB_RASTERS=1
 However, it generated an error message saying that the variables weren't
 recognized.  Can anyone tell me how to set them?
 
 Thanks,
 
 John
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_SelectByValue

[Pierre Racine]

If you set the values you're looking for to themselves and other pixels to 
nodata, ST_Reclass is equivalent to ST_SelectByValue.

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of David Haynes
 Sent: Thursday, June 05, 2014 3:34 PM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] ST_SelectByValue
 
 reclass would allow me to change the values of the pixel.
 
 ST_SelectByValue is supposed to like a map algebra or in ArcGIS Terms
 select by attribute.
 
 
 On Thu, Jun 5, 2014 at 1:38 PM, Pierre Racine
 pierre.rac...@sbf.ulaval.ca wrote:
 
 
   Try ST_Reclass().
 
 
-Original Message-
From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-
 users-
boun...@lists.osgeo.org] On Behalf Of David Haynes
Sent: Friday, May 30, 2014 10:49 AM
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] ST_SelectByValue
   
Hello,
   
I am fairly new to using rasters in postgis, but I am wondering if
 there is a
faster way to extract values from categorical rasters (e.g.
 landcover types)
than using map algebra. In the PostGIS in Action book, pg 404 has
 a
function called ST_SelectByValue(raster|geometry, 'expression')
 which
seems to incorporate map algebra, however I can't seem to tell if
 this
function is finished and available in postgis as I get an unknown
 function
error.
   
My version is this.
PostgreSQL 9.3.2 on x86_64-unknown-linux-gnu, compiled by gcc
(Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit POSTGIS=2.1.1
 r12113
GEOS=3.3.3-CAPI-1.7.4 PROJ=Rel. 4.7.1, 23 September 2009
GDAL=GDAL 1.9.0, released 2011/12/29 LIBXML=2.7.8
 LIBJSON=UN
(...)
   
   
And I found a webpage
   
 http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking03
   
search for ST_SelectByValue at the bottom, indicating that it
 was finished.
   
   
AAccomplished Objectives
 
   
 http://trac.osgeo.org/postgis/wiki/WKTRaster/SpecificationWorking03#A
ccomplishedObjectives
 
   
   
ST_ValueCount(raster, value) - integer - done see below
   
ST_ValuePercent(raster, value) - double precision - done see
 below
   
ST_Resample(raster, method, originx, originy, pixelsizex,
 pixelsizey) - raster
- done see below
   
ST_SelectByValue(raster, 'expression') - same type as first
 argument
   
ST_Clip(raster|geometry,geometry) - same type as first
 argument
 
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] getting raster values for point from st_value

[Pierre Racine]

You said you could visualize this table in QGIS? Or were you speaking about a 
shapefile? You can visualize the rastextets table and the points table the same 
way.

Pierre

 -Original Message-
 From: Jason Mathis [mailto:jmat...@redzonesoftware.com]
 Sent: Tuesday, May 13, 2014 6:30 PM
 To: PostGIS Users Discussion; Pierre Racine
 Subject: Re: [postgis-users] getting raster values for point from st_value
 
 Hi Pierre,
 
 Thanks so much for super quick response. You will have to forgive me I am
 fairly new to postgis.
 
 Creating the table is no problem but how would I check if it overlaps with
 my points?
 
 thanks!
 
 
 
 On May 13, 2014 at 4:08:40 PM, Pierre Racine (pierre.rac...@sbf.ulaval.ca)
 wrote:
 
 
   Create a new table with the extents of the raster tiles and check if it
 overlaps with your points:
 
   CREATE TABLE rastextets AS
   SELECT rast::geometry
   FROM tbl_raster
 
   Pierre
 
-Original Message-
From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-
 users-
boun...@lists.osgeo.org] On Behalf Of Jason Mathis
Sent: Tuesday, May 13, 2014 5:56 PM
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] getting raster values for point from
 st_value
   
Hi all,
   
I have recently loaded a good amount of rasterized shape files into
 a db. I
mainly used gdal_rasterize and raster2pgsql. The shape files were
 in NAD83
so I used the below commands to rasterized and load.
   
gdal_rasterize -at -a desc -where desc='Very High' -burn 2 -tr 30
 30 -ot
byte -a_nodata -99 shape_file.shp raster.tiff
raster2pgsql -d -s 5070 -t 100x100 -F -I -C -Y raster.tiff tbl_raster 
 |
 psql -d
raster_db
   
I had two sets of shape files one set had about 10 burn values and
 another
had about 50. The issue is getting these values out of the raster. I
 believe I
should be using something from this page:
http://postgis.net/docs/RT_ST_Value.html.
   
I can visualize this file in QGIS and it looks good, the histogram
 shows
values. I can also use the postgis functions, st_summarystats or
st_histogram to see some data but I need to be able to query a
 given point
and return the raster value.
   
I have this:
   
SELECT ST_Value(rast,(ST_Point( -111.750185, 34.886948)))
FROM tbl_raster
WHERE st_intersects(rast, (ST_Point(-111.750185, 34.886948)))
   
Then i thought oh wait should I set the srid?:
   
SELECT ST_Value(rast,(ST_SetSRID(ST_Point( -111.750185,
34.886948),5070)))
FROM tbl_raster
WHERE st_intersects(rast, (ST_SetSRID(ST_Point(-111.750185,
34.886948),5070)))
   
Still I get nothing, then I thought I need to transform so I got this:
   
SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -
 111.762866,
34.874309),5070),4326)))
FROM tbl_raster
WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point( -
111.762866, 34.874309),5070),4326)))
   
But in the end nothing.
   
I did need to add the SRID to the spatial_ref_sys table:
   
INSERT into spatial_ref_sys (srid, auth_name, auth_srid,
 proj4text, srtext)
values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5
 +lat_0=23
+lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80
 +towgs84=0,0,0,0,0,0,0
+units=m +no_defs ', 'PROJCS[NAD83 / Conus
   
 Albers,GEOGCS[NAD83,DATUM[North_American_Datum_1983,SPHER
OID[GRS
   
 1980,6378137,298.257222101,AUTHORITY[EPSG,7019]],TOWGS84[0,
   
 0,0,0,0,0,0],AUTHORITY[EPSG,6269]],PRIMEM[Greenwich,0,AUTHORIT
   
 Y[EPSG,8901]],UNIT[degree,0.0174532925199433,AUTHORITY[EPSG
   
 ,9122]],AUTHORITY[EPSG,4269]],PROJECTION[Albers_Conic_Equal_
   
 Area],PARAMETER[standard_parallel_1,29.5],PARAMETER[standard_par
   
 allel_2,45.5],PARAMETER[latitude_of_center,23],PARAMETER[longitude
_of_center,-
   
 96],PARAMETER[false_easting,0],PARAMETER[false_northing,0],UNIT[
   
 metre,1,AUTHORITY[EPSG,9001]],AXIS[X,EAST],AXIS[Y,NORTH],AUT
HORITY[EPSG,5070]]’);
   
But ultimately I have nothing I can’t seem to get any data out
 using a query
with a point value. I am thinking it has to do with the srid? Since I
 can
visually inspect it and see data.
   
Anyone have any hot ideas as to why I am not seeing data?
   
Thanks!
   
   
   
This transmission contains confidential and privileged information
 intended
solely for the party identified above. If you receive this message in
 error,
you must not use it or convey it to others. Please destroy it
 immediately
and contact the sender at (303) 386-3955 or by return e-mail to
 the sender

Re: [postgis-users] add foreign-key constraint on raster import

[Pierre Racine]

You can add the filename as a new column with the -f option. Then you can 
extract whatever you want from this column to build a key.

Does this works for you?

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Max Demars
 Sent: Wednesday, May 14, 2014 10:08 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] add foreign-key constraint on raster import
 
 Hi,
 
 
 I want to import many rasters from files into a same tablewith a foreign-key
 constraint to be able to query them individually or grouped as needed.
 
 
 I was hopping that the raster2pgsql tool would be able to do it, but no
 chance.
 
  Is there a way to ask raster2pgsql to add a foreign-key during import? If
 not, what would be the best approach for this purpose?
 
 Thank you very much,
 
 
 -Max Demars
 
 
 --
 
 Stack Overflow: http://stackoverflow.com/users/1914034/burton449
 
 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449
 
 LastFm: http://www.lastfm.fr/user/burton449

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] add foreign-key constraint on raster import

[Pierre Racine]

Shouldn't it be up to your app to add (if it does not exist yet) and set a key 
column every time a new raster is loaded? Or even better change the name of the 
raster before it is loaded so it gets a unique filename based on you users 
properties? The logic of identifying the files should be done before importing 
in the DB.

Pierre

 -Original Message-
 From: Max Demars [mailto:burton449...@gmail.com]
 Sent: Wednesday, May 14, 2014 11:19 AM
 To: PostGIS Users Discussion
 Cc: Pierre Racine
 Subject: Re: [postgis-users] add foreign-key constraint on raster import
 
 I was thinking about this way to do it, but my app has many users, and in
 case two rasters with the same filename are loaded, there would be a
 problem using filename to create a foreign key. The aim of the foreign key
 here would be to specify which user is the owner of which raster.
 
 
 I suppose I could find a way to handle the possibility of two rasters with the
 same filename, but if I could create the foreign key during the import, I
 think that would be the best.
 
 
 
 On Wed, May 14, 2014 at 11:08 AM, Jason Mathis
 jmat...@redzonesoftware.com wrote:
 
 
   That what I was thinking too. Just to clarify its a capital “-F”.
 
   Its a handy option for post processing the files in the db.
 
   http://postgis.net/docs/using_raster_dataman.html#RT_Raster_Lo
 ader
 
 
 
   On May 14, 2014 at 9:04:32 AM, Pierre Racine
 (pierre.rac...@sbf.ulaval.ca) wrote:
 
 
   You can add the filename as a new column with the -f
 option. Then you can extract whatever you want from this column to build
 a key.
 
   Does this works for you?
 
   Pierre
 
-Original Message-
From: postgis-users-boun...@lists.osgeo.org
 [mailto:postgis-users-
boun...@lists.osgeo.org] On Behalf Of Max Demars
Sent: Wednesday, May 14, 2014 10:08 AM
To: postgis-users@lists.osgeo.org
Subject: [postgis-users] add foreign-key constraint on
 raster import
   
Hi,
   
   
I want to import many rasters from files into a same
 tablewith a foreign-key
constraint to be able to query them individually or
 grouped as needed.
   
   
I was hopping that the raster2pgsql tool would be able to
 do it, but no
chance.
   
Is there a way to ask raster2pgsql to add a foreign-key
 during import? If
not, what would be the best approach for this purpose?
   
Thank you very much,
   
   
-Max Demars
   
   
--
   
Stack Overflow:
 http://stackoverflow.com/users/1914034/burton449
   
GIS Overflow:
 http://gis.stackexchange.com/users/14426/burton449
   
LastFm: http://www.lastfm.fr/user/burton449
 
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-
 users
 
 
 
   This transmission contains confidential and privileged information
 intended solely for the party identified above. If you receive this message in
 error, you must not use it or convey it to others. Please destroy it
 immediately and contact the sender at (303) 386-3955
 tel:%28303%29%20386-3955  or by return e-mail to the sender.
 
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 
 
 
 
 --
 
 Stack Overflow: http://stackoverflow.com/users/1914034/burton449
 
 GIS Overflow: http://gis.stackexchange.com/users/14426/burton449
 
 LastFm: http://www.lastfm.fr/user/burton449

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Postgis raster tutorial

[Pierre Racine]

Done!

http://geospatialelucubrations.blogspot.ca/2014/05/a-guide-to-rasterization-of-vector.html

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of georgew
 Sent: Monday, April 21, 2014 6:37 PM
 To: postgis-users@lists.osgeo.org
 Subject: Re: [postgis-users] Postgis raster tutorial
 
 Thank you Pierre, your add-ons are a life saver. I'll make sure to make the
 most of them.
 And I look forward to your blog in 2000 minutes, the countdown has
 started
 already!
 
 
 
 
 --
 View this message in context: http://postgis.17.x6.nabble.com/Postgis-
 raster-tutorial-tp5006124p5006169.html
 Sent from the PostGIS - User mailing list archive at Nabble.com.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] getting raster values for point from st_value

[Pierre Racine]

Create a new table with the extents of the raster tiles and check if it 
overlaps with your points:

CREATE TABLE rastextets AS
SELECT rast::geometry
FROM tbl_raster

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Jason Mathis
 Sent: Tuesday, May 13, 2014 5:56 PM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] getting raster values for point from st_value
 
 Hi all,
 
 I have recently loaded a good amount of rasterized shape files into a db. I
 mainly used gdal_rasterize and raster2pgsql. The shape files were in NAD83
 so I used the below commands to rasterized and load.
 
 gdal_rasterize -at -a desc -where desc='Very High' -burn 2 -tr 30 30 -ot
 byte -a_nodata -99 shape_file.shp raster.tiff
 raster2pgsql -d -s 5070 -t 100x100 -F -I -C -Y raster.tiff tbl_raster | psql 
 -d
 raster_db
 
 I had two sets of shape files one set had about 10 burn values and another
 had about 50. The issue is getting these values out of the raster. I believe I
 should be using something from this page:
 http://postgis.net/docs/RT_ST_Value.html.
 
 I can visualize this file in QGIS and it looks good, the histogram shows
 values. I can also use the postgis functions, st_summarystats or
 st_histogram  to see some data but I need to be able to query a given point
 and return the raster value.
 
 I have this:
 
   SELECT ST_Value(rast,(ST_Point( -111.750185, 34.886948)))
   FROM  tbl_raster
   WHERE st_intersects(rast, (ST_Point(-111.750185, 34.886948)))
 
 Then i thought oh wait should I set the srid?:
 
   SELECT ST_Value(rast,(ST_SetSRID(ST_Point( -111.750185,
 34.886948),5070)))
   FROM  tbl_raster
   WHERE st_intersects(rast, (ST_SetSRID(ST_Point(-111.750185,
 34.886948),5070)))
 
 Still I get nothing, then I thought I need to transform so I got this:
 
   SELECT ST_Value(rast,(ST_Transform(ST_SetSRID(ST_Point( -111.762866,
 34.874309),5070),4326)))
   FROM   tbl_raster
   WHERE st_intersects(rast, (ST_Transform(ST_SetSRID(ST_Point( -
 111.762866, 34.874309),5070),4326)))
 
 But in the end nothing.
 
 I did need to add the SRID to the spatial_ref_sys table:
 
   INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext)
 values ( 5070, 'EPSG', 5070, '+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23
 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
 +units=m +no_defs ', 'PROJCS[NAD83 / Conus
 Albers,GEOGCS[NAD83,DATUM[North_American_Datum_1983,SPHER
 OID[GRS
 1980,6378137,298.257222101,AUTHORITY[EPSG,7019]],TOWGS84[0,
 0,0,0,0,0,0],AUTHORITY[EPSG,6269]],PRIMEM[Greenwich,0,AUTHORIT
 Y[EPSG,8901]],UNIT[degree,0.0174532925199433,AUTHORITY[EPSG
 ,9122]],AUTHORITY[EPSG,4269]],PROJECTION[Albers_Conic_Equal_
 Area],PARAMETER[standard_parallel_1,29.5],PARAMETER[standard_par
 allel_2,45.5],PARAMETER[latitude_of_center,23],PARAMETER[longitude
 _of_center,-
 96],PARAMETER[false_easting,0],PARAMETER[false_northing,0],UNIT[
 metre,1,AUTHORITY[EPSG,9001]],AXIS[X,EAST],AXIS[Y,NORTH],AUT
 HORITY[EPSG,5070]]’);
 
 But ultimately I have nothing I can’t seem to get any data out using a query
 with a point value. I am thinking it has to do with the srid? Since I can
 visually inspect it and see data.
 
 Anyone have any hot ideas as to why I am not seeing data?
 
 Thanks!
 
 
 
 This transmission contains confidential and privileged information intended
 solely for the party identified above. If you receive this message in error,
 you must not use it or convey it to others. Please destroy it immediately
 and contact the sender at (303) 386-3955 or by return e-mail to the sender.

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Dted elevations

[Pierre Racine]

What happened when you remove and st_value (dted_elevations.rast, f.geom) is 
not null? You get null values?

Do you have many nodata values in your raster?

How is your raster tiled? Did you build a spatial index on the raster tiles?

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Gold, Jack L (US SSA)
 Sent: Friday, May 09, 2014 1:01 PM
 To: 'postgis-users@lists.osgeo.org'
 Subject: [postgis-users] Dted elevations
 
 I know this question has already been answered a hundred times but I am
 at a remote location currently with only email access right now and I hope
 someone can help.  The following query is taking close to a second to run
 and I think it should be much faster.
 
 Select st_value from (with f as ( select st_transform ( st_SetSRID(
 st_MakePoint('106','32'), 4326) as geom) select
 st_value(dted_elevations.rast, f.geom) from dted_elevations cross join f
 where st_intersects ( dted_elevations.rast , f.from) and st_value
 (dted_elevations.rast, f.geom) is not null) as st_value;
 
 Any help is greatly appreciated.
 
 Jack Gold
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[postgis-users] Raster ST_Union memory limit

[Pierre Racine]

Hi,

What maximum size the result of a ST_Union(rast) can be if it is embedded in a 
ST_Tile() call? 

SELECT ST_Tile(ST_Union(rast), 100, 100) rast FROM ...

The RAM available to PostgreSQL or the maximum size of a PostgreSQL field (1GB)?

Thanks,

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Raster ST_Union memory limit

[Pierre Racine]

Thanks, 

The RAM limit is per row or per query? In other word, if a query produce fields 
smaller than my available RAM but the whole resulting table is bigger than my 
available RAM, are each resulting row flushed to the disk so that my whole 
query does not explode?

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Bborie Park
 Sent: Tuesday, May 06, 2014 9:52 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] Raster ST_Union memory limit
 
 RAM available.
 
 -bborie
 
 
 On Tue, May 6, 2014 at 6:35 AM, Pierre Racine
 pierre.rac...@sbf.ulaval.ca wrote:
 
 
   Hi,
 
   What maximum size the result of a ST_Union(rast) can be if it is
 embedded in a ST_Tile() call?
 
   SELECT ST_Tile(ST_Union(rast), 100, 100) rast FROM ...
 
   The RAM available to PostgreSQL or the maximum size of a
 PostgreSQL field (1GB)?
 
   Thanks,
 
   Pierre
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Raster ST_Union memory limit

[Pierre Racine]

I tend to think the opposite since I already produced 30 GB tables on a machine 
with 8 GB of RAM. I guess once finished, every row is flushed to the disk 
and its RAM share freed. 

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Bborie Park
 Sent: Tuesday, May 06, 2014 10:03 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] Raster ST_Union memory limit
 
 I think your query will explode...
 
 
 On Tue, May 6, 2014 at 7:01 AM, Pierre Racine
 pierre.rac...@sbf.ulaval.ca wrote:
 
 
   Thanks,
 
   The RAM limit is per row or per query? In other word, if a query
 produce fields smaller than my available RAM but the whole resulting table
 is bigger than my available RAM, are each resulting row flushed to the disk
 so that my whole query does not explode?
 
   Pierre
 
 
-Original Message-
From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-
 users-
boun...@lists.osgeo.org] On Behalf Of Bborie Park
Sent: Tuesday, May 06, 2014 9:52 AM
To: PostGIS Users Discussion
Subject: Re: [postgis-users] Raster ST_Union memory limit
   
RAM available.
   
-bborie
   
   
On Tue, May 6, 2014 at 6:35 AM, Pierre Racine
pierre.rac...@sbf.ulaval.ca wrote:
   
   
  Hi,
   
  What maximum size the result of a ST_Union(rast) can be if it
 is
embedded in a ST_Tile() call?
   
  SELECT ST_Tile(ST_Union(rast), 100, 100) rast FROM ...
   
  The RAM available to PostgreSQL or the maximum size of a
PostgreSQL field (1GB)?
   
  Thanks,
   
  Pierre
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
   
   
 
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Postgis raster tutorial

[Pierre Racine]

There are two strategies to convert a table composed of many geometries to a 
raster coverage:

1) ST_AsRaster each geometries, ST_Union them and then ST_Tile them.

2) Use ST_ExtractToRaster() from the PostGIS Add-ons.

I plan to blog about this as soon as I get 400 times 5 minutes...

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of georgew
 Sent: Thursday, April 17, 2014 7:25 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] Postgis raster tutorial
 
 As a new raster user I have found it extremely difficult to come to grips
 with the task of creating a raster table with postgis, in particular from an
 existing vector table, a task which is really trivial when done with some
 existing s/w packages ( i.e. SAGA, GRASS etc..). Most of the examples I
 could find were based  on the use of raster2pgsql which implies the prior
 existence of a raster file created with an external package, not on using
 postgis native capabilities to create such a raster. If you already have a
 raster file,  postgis has no problem creating a table from it, but what if
 you don't? Should you be forced to use one of the external packages?
 The postgis manual has a total of two (2) examples on how to create a
 raster
 table with st_asraster, although there are ten (10) different variants. My
 point is not a criticism of postgis raster, at the contrary, a suggestion
 that it could be much more useful if properly explained to users not
 familiar with it. In other words a good postgis raster tutorial would go a
 long way to convince users to at least try postgis raster by itself without
 recourse to external packages (except perhaps for the visualisation of the
 rasters). Specifically, how to create with postgis a raster table from a
 vector table, how to extract values from one of the vector attributes and
 store them in the raster, resulting in a raster identical to what would have
 been created with raster2pgsql. I, with my ignorance of the topic, have not
 been able to do that yet.
 If such a tutorial has already appeared somewhere I would be delighted to
 be
 pointed in that direction.
 Many thanks
 
 
 
 
 --
 View this message in context: http://postgis.17.x6.nabble.com/Postgis-
 raster-tutorial-tp5006124.html
 Sent from the PostGIS - User mailing list archive at Nabble.com.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_AsRaster vs raster2pgsql

[Pierre Racine]

It's not because the entries in raster_columns are different that the actual 
tables are different. Info in raster_columns is generated from constraints set 
by raster2pgsql. A better way to compare the two tables is to use ST_Metadata() 
and ST_BandMetadata().

You can also set the constraints on the second table using 
AddRasterConstraints(). Then have a look at raster_columns.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of georgew
 Sent: Tuesday, March 04, 2014 12:51 PM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] ST_AsRaster vs raster2pgsql
 
 Hi, I am new to raster and am having some problems creating a raster table
 from a contour vector table, PostGIS 2.1.1, Win 8.1 64bit
 
 To create the table I am using:
 
 CREATE TABLE BK30_dem2 AS
   SELECT 1 AS rid, ST_AsRaster((
SELECT
   ST_Collect(geom)
FROM contours_topo_2 where sheet='BK30'
), 5.0, 5.0 )
   AS rast;
 
 For comparison purpose I also created a grid from the same contours using
 SAGA GIS, then created a second raster table in PostGIS using:
 
  raster2pgsql -I -C -s 2193 I:\dems\pg-50k-grid\BK30\BK30_dem2.sdat
 public.bk30_saga | psql -d NZTPU;
 
  I was hoping for identical results but the output from Raster Columns
 shows
 the following:
 
 
 NZTPU;public;bk30_saga;rast;2193;5;-
 5;4801;5784;TRUE;FALSE;1;{32BF};{-
 9};{f};010320910801000580DD0F3
 A41CECF067D3C7255410080A26D3A41CECF067D3C725541008
 0A26D3A41CECF067DFE410080DD0F3A41CECF067DFE410
 080DD0F3A41CECF067D3C725541
 
 NZTPU;public;bk30_dem2;rast;0;FALSE;FALSE;
 
 Clearly I am doing something wrong, but what??
 Many thanks for any help.
 
 
 
 --
 View this message in context: http://postgis.17.x6.nabble.com/ST-
 AsRaster-vs-raster2pgsql-tp5005827.html
 Sent from the PostGIS - User mailing list archive at Nabble.com.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_AsRaster vs raster2pgsql

[Pierre Racine]

You can set the pixeltype and the extent in ST_AsRaster(). You seems also to 
have a different nodata value that you can also set in ST_AsRaster(). So a bit 
of reading about ST_AsRaster() and experiment and you'll should get what you 
want.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of georgew
 Sent: Tuesday, March 04, 2014 4:55 PM
 To: postgis-users@lists.osgeo.org
 Subject: Re: [postgis-users] ST_AsRaster vs raster2pgsql
 
 Thank you Pierre, after adding constraints I get:
 
 NZTPU;public;bk30_saga;rast;2193;5;-
 5;4801;5784;TRUE;FALSE;1;{32BF};{-
 9};{f};010320910801000580DD0F3
 A41CECF067D3C7255410080A26D3A41CECF067D3C725541008
 0A26D3A41CECF067DFE410080DD0F3A41CECF067DFE410
 080DD0F3A41CECF067D3C725541
 
 NZTPU;public;bk30_dem2;rast;2193;5;-
 5;4800;5783;TRUE;FALSE;1;{8BUI};{0};{f};0103209108010
 00500E00F3A413C725541A06D3
 A413C725541A06D3A410040FF41000
 0E00F3A410040FF41E00F3A413C725541
 
 much better! but not identical. The blocksize is slightly different, the
 pixel type and no_data value are different and the  extent geometry is also
 slightly different, I assume all due to the pixel type being generated
 differently. Using st_metadata and st_bandmetadata I get:
 
 (1707997.5,5622001.95354076,4801,5784,5,-5,0,0,2193,1);(32BF,-
 9,f,)
 saga
 (1708000,5622000,4800,5783,5,-5,0,0,2193,1);(8BUI,0,f,)
 dem2
 
 Any suggestions as to what I can do to make them identical?
 
 
 
 --
 View this message in context: http://postgis.17.x6.nabble.com/ST-
 AsRaster-vs-raster2pgsql-tp5005827p5005830.html
 Sent from the PostGIS - User mailing list archive at Nabble.com.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] import raster

[Pierre Racine]

Importing a big raster as a coverage with a huge number of small tiles always 
takes a long time. Two options:

-Increase the size of your tiles. They don't have to be that small to get good 
performance with intersections.

-Do you really need the raster to be inside the database? Try adding the -R 
raster2pgsql option. The import step will be way faster and up to now we 
haven't seen much difference in terms of intersection/extraction performance 
when the raster data is kept outside the db.

Pierre

From: postgis-users-boun...@lists.osgeo.org 
[postgis-users-boun...@lists.osgeo.org] On Behalf Of Eloi [e...@openmailbox.org]
Sent: Thursday, February 27, 2014 5:03 AM
To: Lista PostGIS users
Subject: [postgis-users] import raster

Hi all,

My SELECT PostGIS_Full_Version() is: POSTGIS=2.0.3 r11128
GEOS=3.3.8-CAPI-1.7.8 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL
1.10.0, released 2013/04/24 LIBXML=2.9.0 LIBJSON=UNKNOWN RASTER

I am importing a GeoTIFF (1 band, Type = Byte,  Not Compressed file size
= 933.5 MB, columns/rows 43200, 21600, Pixel Size =
0.008,-0.008) into PostGIS and is taking more
than one hour by now. Is that normal?

This is the command I have used:
raster2pgsql -s 4326 -d -I -C -M -t 5x5 /path/file_name.tif
schema_name.table_name|psql -d db_name

Yes, very small tiles. The goal is to intersect this raster with point
features. But the issue right now is the import it self that is taking
to much time.
I can see now that the first 1:30 were for inserts and now is starting
with the constraints. Is this the normal speed or am I missing
something?

Thanks. Cheers,

Eloi
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_MapAlgebraExpr 2 rasters problem with resultant raster

[Pierre Racine]

What is the exact query you do when trying to compute the sum or the max of the 
resulting raster?

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of g.si...@utwente.nl
 Sent: Monday, February 24, 2014 3:44 PM
 To: postgis-users@lists.osgeo.org; postgis-users-ow...@lists.osgeo.org
 Subject: [postgis-users] ST_MapAlgebraExpr 2 rasters problem with
 resultant raster
 Importance: High
 
 Hi,
 
 
 
 I am using ST_MapAlgepraExpr for the intersection of two rasters.
 
 The first raster has minimum and maximum values as 5.112e-020 and
 0.000125 respectively.
 
 The second raster has minimum and maximum value as 0.000 and
 0.000400 respectively.
 
 
 
 I carry the raster intersection using expression
 
 ST_MapAlgebraExpr(a.rast1,b.rast2, '[rast1.val] * [rast2.val]', '32BF',
 'INTERSECTION', 'NULL', 'NULL')
 
 And got a resultant raster whose value ranges from 0.000 to 1.862e-008.
 
 
 
 When I query the resultant raster for sum(ST_Summarystats(rast)) or
 max(ST_Summarystats(rast)) I get the blank cell with double precision
 datatype.
 
 
 
 I can visualize the raster and can see the sum of all cell values, max and min
 cellvalues in ILWIS software but can not calculate the same using postgis
 queries.
 
 
 
 Is there something wrong in my above map algebra expression or there is a
 problem with map algebra function?
 
 
 
 Can someone please provide any solution to this problem.
 
 
 
 Thanks in Advance.
 
 Gaurav

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Masking raster values

[Pierre Racine]

Guido,

no masks (you can create a mask as a band) means that a mask is nothing 
different than a normal band. It's just a band that is mostly used to intersect 
(mask) another band having normal values. You can do that with... ... 
ST_Intersection(raster, raster)

So first you want to ST_MapAlgebra() your BQA to create a raster with 0 when 
those bits are set and 1 when they are not. You could use a SQL expression like 
CASE [rast]  61440 THEN 0 ELSE 1 END in your ST_MapAlgebra() call.

CREATE TABLE mask AS
SELECT ST_MapAlgebraExpr(rast, bandnumber, 1BUI, CASE [rast]  61440 THEN 0 
ELSE 1 END)
FROM yourLandsatTable

You can make this query faster by implementing a small callback PL/PGSQL 
function for the ST_MapAlgebra() variant taking a callback 
(http://postgis.net/docs/manual-2.1/RT_ST_MapAlgebra.html).

Second you want to intersect this mask with you some of your data band. Look 
at ST_Intersection(raster, raster).

CREATE TABLE newLandsat AS
SELECT ST_Intersection(a.rast, BAND_NUM, b.rast, 1, 'BAND1') rast
FROM yourLandsatTable a, mask b

You can do the tow operations in a single SQL statement but it is generally 
wiser and faster to make it in two.

As your question about creating raster operators I always found this idea very 
nice in theory but in practice you generally want to/have to be able to set a 
variety of options when doing operations involving many rasters: What should be 
the pixel type of the result? What to do with nodata values? Do you want the 
result to meet the extent of the first or the second raster, the union of them 
or the intersection of them? And so on... For sure you can establish default 
behavior for all of these in order to be able to construct nice and clean 
raster expressions, but then you become very restricted. PostGIS took the side 
of providing all the flexibility possible at the cost of sometimes hard to 
write function calls.

Hope this answer your questions.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Guido LEMOINE
 Sent: Wednesday, February 12, 2014 12:19 PM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] Masking raster values
 
 Dear List,
 
 
 
 I am dabbling with some Landsat-8 imagery in PostGIS, using raster
 functionality. I am primarily interested in deriving polygon-delineated
 extracts from the various spectral bands (e.g. for a forest patch, an
 agricultural field, etc).
 
 
 
 Landsat-8 provides the so-called BQA band (see
 http://landsat.usgs.gov/L8QualityAssessmentBand.php) which is a bit-
 coded set of quality parameters, for which the cloud and haze indicators
 (bits 12-15 (right to left)) are the most important (for me). I would want to
 mask out all pixels, in the spectral band(s), for which these bits are set in
 the BQA band.
 
 
 
 Did anyone come across a solution that would address this kind of masking
 operation? The only reference I seem to be able to find is the statement
 that no masks (you can create a mask as a band) in the WKTRaster wiki
 page on osgeo, which is somewhat cryptic. I would know how to do this
 outside the database (using JAI), but it seems that masking is some core
 operation one would expect in raster functionality.
 
 
 
 As a side issue, would it be possible (is it foreseen?) to create raster
 operators that work in an equivalent arithmetic way as on single variable
 (i.e. 2*rast, rast ~ 128, rast  3) but produce rast as output (maybe through
 a mapping to ST_MapAlgebra functions)?
 
 
 
 Any pointers welcome and excuses if I have overlooked discussions on this
 topic.
 
 
 
 Guido Lemoine
 
 
 
 Scientific Officer
 
 European Commission, Joint Research Centre (JRC)
 
 Institute for the Protection and Security of the Citizen (IPSC)
 
 Global Security and Crisis Management Unit
 
 
 
 Via E.Fermi 2749, I-21027 Ispra (VA) Italy, TP 268
 
 Tel. +39 0332 786239 (direct line) Fax +39 0332 785154
 
 e-mail: guido.lemo...@jrc.ec.europa.eu
 mailto:guido.lemo...@jrc.ec.europa.eu
 
 web: http://globesec.jrc.ec.europa.eu http://globesec.jrc.ec.europa.eu/
 
 
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Removing tiny polygons

[Pierre Racine]

Give a try to ST_TrimMulti() in the PostGIS Add-ons:

http://geospatialelucubrations.blogspot.ca/2013/11/launching-postgis-add-ons-new-postgis.html

https://github.com/pedrogit/postgisaddons/releases/tag/1.21

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Astrid Bjørnerød
 Sent: Tuesday, January 28, 2014 1:19 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] Removing tiny polygons
 
 Is there anyone you have written some sql-code for removing the small
 polygons onto one of the neighbor-polygons ?
 Any stabil functions for this task?
 
 Astrid
 
 -Opprinnelig melding-
 Fra: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] På vegne av Tom McCallum
 Sendt: 27. januar 2014 20:48
 Til: mor...@sickel.net; PostGIS Users Discussion
 Emne: Re: [postgis-users] Removing tiny polygons
 
 Thanks Morten.  I wanted the first case, so that helps.
 
 On Mon, 27 Jan 2014 19:15:18 -, Morten Sickel mor...@sickel.net
 wrote:
 
  It depends on what you mean with remove - do you want to just have
  the large polygons without the area that is covered by the small ones
  or do you want to add the area of the small polygons onto one or more
  of the larges? (e.g. how one typically wants to handle shiver
  polygons)
 
  In the first case it is as easy as a query like
 
  select the_geom where ST_area(the_geom)  ...
 
  Morten
 
 
  Tom McCallum skrev:
  Hi all,
 
  If I have a polygon geometry field which seems to be made up of a
  large number of tiny polygons and some larger ones, what is the best
  way to remove polygons below a certain size?  Would I need to
  preprocess this and save in a new table/record or can I dynamically
  do it?
 
  Thanks
 
  Tom
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 
 
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_ConcaveHull giving up too soon, aka too convex

[Pierre Racine]

Try this:

SELECT ST_Buffer(ST_Union(ST_Buffer(geom, 0.0001)), -0.8) FROM points

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Ivan Price
 Sent: Tuesday, January 14, 2014 8:36 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] ST_ConcaveHull giving up too soon, aka too
 convex
 
 Thanks for your reply Remi,
 
 
 
 for now we are going with simply a union of a buffer of the points (as you
 suggested) as it makes sense in the context of the origin of the data
 anyway.
 
 
 
 cheers
 
 
 
 -i
 
 
 
 
 
 De : postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] De la part de Rémi Cura
 Envoyé : Tuesday, 14 January 2014 11:39
 À : PostGIS Users Discussion
 Objet : Re: [postgis-users] ST_ConcaveHull giving up too soon, aka too
 convex
 
 
 
 Hey,
 
 sadly I can't answer your questions,
 
 but it looks like you don't want concave hull but instead Alpha Shape (see
 SFCGAL , maybe they have integrated it).
 
 
 You have still a lot of workaround :
 if you don't mind expending the bouyndary, a union of buffers of points
 if you don't mind quantization (anyway your input looks quantized) :
 convert to raster, then any classical image processing easily available
 (watershed for example) (you can use out of the box QGIS 2.0 raster
 processing tools)
 if you want a sharp boudary : and assuming your points are regularly
 spaced : maybe something liek this (untested)
 create a segment for each point to his K nearest neighbours (K would be 2
 or 3, and you would put a max distance before drawing a line)
 
 Then pick the segments that have no  intersection with other segments, and
 either use ST_Polygonyze, or cut the bounding box polygon with the
 segments and remove the part outside
 .
 
 Maybe you can try to reduce your point precision for processing (snap to
 grid) and/or put the points in an appropriate srid (not too much digits)
 
 
 
 Cheers,
 
 Rémi-C
 
 
 
 2014/1/14 Ivan Price ivan.pr...@noveltis.fr
 
 
 
 Hi there,
 
 
 
 I am experimenting with ST_ConcaveHull, am using Postgres 9.2.2 on
 windows with postgis 2.1.1 r12113
 
 
 
 i have a table of 218 points forming roughly an arc (something like the nike
 swoosh)
 
 
 
 when i ask for the concave hull i'm expecting a poly shrink-wrapped to
 them, however this is not the case, even when i set a very low percentage
 area parameter. for example:
 
 
 
 SELECT ST_ConcaveHull(ST_Collect(geom), 0.01) as geom FROM
 points;
 
 
 
 i attach a screenshot of the points and the convext hull, as well as an sql
 dump of the points themselves.
 
 
 
 also as a side note, in researching the problem i upgraded from postgis 2.01
 to 2.1.1 and noticed a change in the processing time for the above
 statement from 1 minute to ~ 5 minutes.
 
 
 
 are there any pointers for making ConcaveHull more 'aggressive' ?
 
 
 
 thanks and regards
 
 
 
 -i
 
 
 
 
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] point to raster conversion

[Pierre Racine]

You did not paste your code...

ST_ExtracToRaster() should work on a geometry table. If you want it to work on 
point cloud table you have to add a set of methods to the 
ST_ExtractPixelValue4ma() function. Sorry that was not clear. you could start 
by converting your patches to geometries to test ST_ExtracToRaster() and then 
try to make the same query work on point patches.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rémi Cura
 Sent: Thursday, December 19, 2013 8:52 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] point to raster conversion
 
 Hey,
 thanks for the answer.
 
 I tried as you say,
 it is not working and slow (several seconds for 50*50 pixels).
 the function ST_ExtracToRaster just seems to do nothing.
 (after using the function, st_value on any pixel return NULL, st_summary
 also)
 
 
 I have 2 bands of '32BF'.
 I set the srid of both geom and raster to be the same.
 
 
 Is there any way to set multiple pixels at once (and not do
 st_set(st_set(..,
 or a band constructor taking multiple pixels (an array for example)?
 
 
 Cheers,
 
 
 Rémi-C
 
 
 Here is my test code :
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] point to raster conversion

[Pierre Racine]

MEAN_OF_VALUES_AT_PIXEL_CENTROID, which is the default ST_ExtractToRaster() 
method, search for geometries intersecting with the centroid of the pixel (a 
point). That's the method you use both calls to ExtractToRaster(). As it's 
almost impossible that a point intersects a point, this is certainly not the 
right method to get the value of points inside the pixel. A proper method would 
be MEAN_VALUE_OF_POINTS or MEAN_VALUE_OF_GEOMETRIES or 
FIRST_GEOMETRY_VALUE which are not implemented...

I suggest you try 'COUNT_OF_POINTS' for now which will set each pixel value to 
the number of point intersecting with it. If you get it to work and you're 
satisfied with the performance I could help adding a method for your specific 
need. 

Can't do anything about the performance. All this is pure pl/pgsql and it's as 
fast as it can. 5 seconds to compute 2500 pixels value using SQL is not that 
bad though if you're not on the web... Otherwise you should have a specific C 
function doing exactly what you want. Expect time and $$$.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rémi Cura
 Sent: Thursday, December 19, 2013 10:27 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] point to raster conversion
 
 
 Trying to exclude possible sources of errors here :
 it tried to convert points to 2D, no sucess.
 
 Cheers,
 Rémi-C
 
 
 
 2013/12/19 Rémi Cura remi.c...@gmail.com
 
 
   Oups,
   I figured it would be of no use.
 
   I of course works on pure geometry, raster functions only sees
 postgis points and ints.
   Here few lines of the point table
 
   qgis_id, st_astext(geom) AS geom,z,reflectance
 
   1 POINT Z (2248.0317003 20680.022999
 49.02441796875) 49024 8110
   2 POINT Z (2247.7656001 20680.251999
 49.0198828125) 49020 7890
   3 POINT Z (2248.2528996 20680.247999
 49.02616015625) 49026 8404
 
 
   Is there any function to work on pixel sets , (like read/write many at
 a time)
 
   Now the code :
   (I can make a test case out of this if you want)
 
   Thanks,
 
   Rémi-C
 
   ///
 
 
   SELECT postgis_full_version();
   --POSTGIS=2.1.0 r11822 GEOS=3.5.0dev-CAPI-1.9.0
 r3963 PROJ=Rel. 4.8.0, 6 March 2012 GDAL=GDAL 1.10.0, released
 2013/04/24 LIBXML=2.8.0 RASTER
 
 
 
 
   converting a patch to raster
   --get patch extend
   --create a raster with this extend
 
   --explode patch to points
   --set pixel value to point intersecting it.
 
 
   --preparing raster :
   --creating the raster table
   DROP TABLE IF EXISTS patch_to_raster;
   CREATE TABLE patch_to_raster(rid serial primary key
   , rast raster);
 
   --
   --DELETE  from patch_to_raster
   --inserting an empty raster
   INSERT INTO patch_to_raster (rast) VALUES (
   ST_MakeEmptyRaster(
   50
   ,50
   ,upperleftx:=0
   ,upperlefty:=0
   ,scalex:=1
   , scaley:=1
   , skewx:=0
   ,skewy:=0
   , srid:=932011
   ) --srid of translated lambert 93 to match laser referential
   );
 
   --setting the correct pixel size
   UPDATE patch_to_raster SET rast =
 ST_SetScale(rast,0.02,0.02);
 
   --checking the content:
   SELECT ST_Summary(rast)
   FROM patch_to_raster;
   --Raster of 50x50 pixels has 0 bands and extent of BOX(0 0,1
 1)
 
 
   --adding band
   --first band is Z : float
   --second band is reflectance : float
   UPDATE patch_to_raster SET rast  = ST_AddBand( --1'st
 band, Z
   rast
   ,pixeltype:='32BF'  -- '32BUI'
   , initialvalue:=NULL
   , nodataval:=NULL
   );
   UPDATE patch_to_raster SET rast  = ST_AddBand( --2nd
 band, reflectance
   rast
   ,  pixeltype:='32BF'  --'32BUI'
   , initialvalue:=NULL
   , nodataval:=NULL
   );
 
 
   --checking the content:
   SELECT ST_Summary(rast) FROM patch_to_raster;
   --Raster of 50x50 pixels has 2 bands and extent of BOX(0 0,1
 1)
   --band 1 of pixtype 32BF is in-db with no NODATA value
   --band 2 of pixtype 32BF is in-db with no NODATA value
 
   --getting a patch and creating a table with it
   DROP TABLE IF EXISTS one_patch_into_points;
   CREATE TABLE one_patch_into_points AS
   WITH patch  AS (
   SELECT gid,patch
   FROM acquisition_tmob_012013.riegl_pcpatch_space
 

Re: [postgis-users] point to raster conversion

[Pierre Racine]

 Can you please confirm if there is or isn't a method to set/get several pixels
 at a time?

I don't know any...
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] point to raster conversion

[Pierre Racine]

Try adding a rid column to your raster table. The GDAL driver was until 
recently dependent on rid.

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rémi Cura
 Sent: Thursday, December 19, 2013 12:06 PM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] point to raster conversion
 
 QGis 2.0.1.
 I tried as well to convert to UINT16
 no sucess
 
 
 2013/12/19 Rémi Cura remi.c...@gmail.com
 
 
   Is there a special trick to display the raster in QGIS?
   I do
   db manage/right click on the raster table/add to canvas
   Seems like there is nothing in it afterward.
 
   Thanks again,
 
   Rémi-C
 
 
   2013/12/19 Rémi Cura remi.c...@gmail.com
 
 
   Buffering make it works.
 
   Cheers,
   Rémi-C
 
 
   2013/12/19 Rémi Cura remi.c...@gmail.com
 
 
   Hey, thanks for the answer, and nice catch !
   I'll try to simply buffer points,
   then try th eother method to count
 
   About performance :
   this is not about plpgsql or C.
   I did quit the same processing when importing files
 of 3 millions points to split into small 1*1*1 m cubes.
   There were a lot of cubes (around 2k for 3 millions
 points), a lot of points (3 millions) with a lot of attributes (around 20
 flaot/double per point), and it was around 60 sec/ file (for the data
 processing. Data reading from disk was little longer).
 
 
   For a 50x50 pixels and a few thousand points I
 would expect around 200ms.
 
 
   I suspect it all boils down to current pixel access
 mechanism.
   Can you please confirm if there is or isn't a method
 to set/get several pixels at a time?
 
   Thanks anyway,
   I wouldn't have found =)
 
   Cheers,
   Rémi-C
 
 
   2013/12/19 Pierre Racine
 pierre.rac...@sbf.ulaval.ca
 
 
   MEAN_OF_VALUES_AT_PIXEL_CENTROID,
 which is the default ST_ExtractToRaster() method, search for geometries
 intersecting with the centroid of the pixel (a point). That's the method you
 use both calls to ExtractToRaster(). As it's almost impossible that a point
 intersects a point, this is certainly not the right method to get the value of
 points inside the pixel. A proper method would be
 MEAN_VALUE_OF_POINTS or MEAN_VALUE_OF_GEOMETRIES or
 FIRST_GEOMETRY_VALUE which are not implemented...
 
   I suggest you try 'COUNT_OF_POINTS' for
 now which will set each pixel value to the number of point intersecting with
 it. If you get it to work and you're satisfied with the performance I could
 help adding a method for your specific need.
 
   Can't do anything about the performance.
 All this is pure pl/pgsql and it's as fast as it can. 5 seconds to compute 
 2500
 pixels value using SQL is not that bad though if you're not on the web...
 Otherwise you should have a specific C function doing exactly what you
 want. Expect time and $$$.
 
 
   Pierre
 
 
 
 
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] point to raster conversion

[Pierre Racine]

I would rather use the user custom function mapalgebra technique used in 
ST_ExtracToRaster() in the PostGIS Add-ons to extract a mean value for the 
group of point intersecting with each pixel. This should work well if 1) most 
pixels have points inside so there is no need to interpolate 2) it is possible 
to ST_Intersects() with the point patches (I didn't check much the Point Cloud 
API) and extract only intersecting points.

If the density of point is so low that many pixels do not have points inside 
then we then interpolate. But we don't have any interpolation method yet other 
than filtering functions implemented with neighbourhood map algebra.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rémi Cura
 Sent: Wednesday, December 18, 2013 10:05 AM
 To: PostGIS Users Discussion
 Subject: [postgis-users] point to raster conversion
 
 Hey,
 
 not very original question :
 I would like to convert parts of a point cloud into a raster.
 
 I'm under the impression that currently this is more easily done outside
 PostGIS (with grass, gdal, R, or wathever).
 
 I this simple process a good idea?
 
 for each small part of point cloud (a patch)
 
   get spatial extend
 
   generate a raster covering this extand with givenpixel size and given
 bands
 
   for pixel which intersects a point, set pixel band values with point
 attributes
 
 
 
 Then (optionnal)
 for each raster
 
   use efficient image processing code to interpolate
 
 
 Cheers,
 
 Rémi-C
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] PostGIS ST_MapAlgebra Assistance

[Pierre Racine]

Assuming your three rasters would be in three different tables and not tiled, 
it could look like this: 

WITH 
solandcov AS (SELECT ST_MapAlgebra(a.rast, b.rast, '([rast1] = 1491 and 
[rast2] = 70)::int', 4BUI) FROM sol_gain_clipped a, can_cov_01 b),
creek AS (SELECT ST_MapAlgebra(a.rast, b.rast, '([rast1] and [rast2] = 
10)::int', 4BUI) FROM solandcov a, imp_sur_10_french_creek b)
SELECT ST_AsTiff(rast) FROM creek;

I'm using the last variant of ST_MapAlgebra from this page:

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

If your rasters ares tiled you have to add a WHERE clause to the two first 
statements to make sure only aligned tiles get involved in ST_MapAlgebra. I 
generally use 

WHERE ST_UpperLeftX(a.rast) = ST_UpperLeftX(b.rast) AND ST_UpperLeftY(a.rast) = 
ST_UpperLeftY(b.rast)

for that your tiles have to be well aligned. If not then just use 
ST_Intersects(a.rast, b.rast) but that should be slower.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Jason Coombs
 Sent: Monday, December 16, 2013 9:22 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] PostGIS ST_MapAlgebra Assistance
 
 Good Morning,
 I would like to perform the equivalent of the below ArcGIS raster calculator
 equation in PostGIS using the ST_MapAlgebra function and return it as a tif
 using ST_AsTIFF.
 
 Con((Sol_gain_clipped.tif = 1491)(can_cov_01.tif = 70)  
 (imp_sur_10_French_Creek.tif = 10),1,0)
 
 The equation uses three single band rasters with pixel value specifications
 for each one, and returns a single band raster with values of 1 for pixels
 meeting all three requirements, and 0 for pixels that do not. Any help in
 how to write the SQL, including PL/pgSQL, would be greatly appreciated.
 
 Best, Jason
 
 
 
 Jason A. Coombs
 
 Department of Environmental Conservation/US Forest Service
 
 201 Holdsworth Hall
 
 University of Massachusetts
 
 Amherst, MA 01003
 
 
 
 https://bcrc.bio.umass.edu/pedigreesoftware/
 https://bcrc.bio.umass.edu/pedigreesoftware/
 
 
 
 Phone: 413-545-1845
 
 Fax: 413-545-1860
 
 Email: jcoo...@cns.umass.edu mailto:jcoo...@cns.umass.edu
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] biggest rectangle contained within another polygon

[Pierre Racine]

You will have to define biggest more precisely. Biggest in terms of area? 
There might be many inner rectangles having the same biggest area. How to 
choose one over the other?

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of tommaso
 Sent: Wednesday, December 04, 2013 5:07 AM
 To: PostGIS Users Discussion
 Subject: [postgis-users] biggest rectangle contained within another polygon
 
 Hello, I need to calculate the biggest rectangle which is contained
 within another polygon. A sort of inner envelope or inner bounding box.
 I found this interesting article:
 http://cgm.cs.mcgill.ca/~athens/cs507/Projects/2003/DanielSud/
 Is it possible in Postgis?
 
 I have attached a picture with a example of what I need.
 
 Tom
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] biggest rectangle contained within another polygon

[Pierre Racine]

In the most simple case, take a very long rectangle and rotate it. There can be 
an infinity of horizontal rectangles inside it all having the same biggest area.

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Pierre Racine
 Sent: Wednesday, December 04, 2013 9:42 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] biggest rectangle contained within another
 polygon
 
 You will have to define biggest more precisely. Biggest in terms of area?
 There might be many inner rectangles having the same biggest area. How to
 choose one over the other?
 
 Pierre
 
  -Original Message-
  From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
  boun...@lists.osgeo.org] On Behalf Of tommaso
  Sent: Wednesday, December 04, 2013 5:07 AM
  To: PostGIS Users Discussion
  Subject: [postgis-users] biggest rectangle contained within another
 polygon
 
  Hello, I need to calculate the biggest rectangle which is contained
  within another polygon. A sort of inner envelope or inner bounding box.
  I found this interesting article:
  http://cgm.cs.mcgill.ca/~athens/cs507/Projects/2003/DanielSud/
  Is it possible in Postgis?
 
  I have attached a picture with a example of what I need.
 
  Tom
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] [postgis] Optimal tile size for Raster

[Pierre Racine]

You are getting fast result because you are not operating at the pixel level. 
Only at the raster extent level (converted to geometry). A more usual query 
would be to intersect a raster with some polygons:

SELECT (ST_Intersection(r.rast, q.geom)).*
FROM ilatlon32x32 as r, geometrytable as q
WHERE ST_Intersects(r.rast, q.geom)

Why do you want to intersect two rasters? Generally it is to mask one raster 
with another.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Heng Zhi Feng (zh...@hsr.ch)
 Sent: Wednesday, December 04, 2013 10:56 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] [postgis] Optimal tile size for Raster
 
 Hello,
 
 
 
 I am having to import Raster (.tiff) into Postgres using the raster2pgsql 
 tool.
 The raster used is 12MB in size and the Postgres version is 9.1 on Ubuntu
 13.10.
 
 
 
 The virtual machine has the following hardware specs.
 
 -  17408 Memory
 
 -  4 Processors
 
 -  16GB HDD
 
 The command used to import (with the parameters) is: raster2pgsql -s
 4326 -I -C -M -R -d -l 4 U:\Desktop\ilatlon_float.tif -F -t 8x8 
 ilatlon8x8|psql -
 d Raster
 
 And the tile size varies from 8x8 all the way up till 2048 (power of 2s).
 
 
 
 My results from the query ran showed that with larger tile size, the query
 has a faster returned timing.
 
 This trend is however opposite to some of the similar tests conducted by
 others (from the internet), where their conclusion was smaller tile size will
 have better performance.
 
 
 
 This is my one of the query I ran:
 
 EXPLAIN (ANALYZE,BUFFERS)
 
 SELECT DISTINCT ST_Intersection(r.rast::geometry,q.rast::geometry)
 
 FROM ilatlon32x32 as r, boundedilatlon32x32 as q;
 
 
 
 Is someone able to give advice on this, if this result that I obtained was
 legit? Or suggest where could the problem be that gave this result?
 
 
 
 Regards,
 
 Zhi Feng

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] [postgis] Optimal tile size for Raster

[Pierre Racine]

I guess maybe you are only looking at PostGIS for its capacity to deal with 
raster in the database. You have to know that rasters and geometries are 
located in space (georeferenced) and might very well intersects if they overlap 
in space. Geometries are generally loaded from shapefiles. To create some 
polygons that overlap your raster just open the latter in a GIS, create a new 
shapefile and add it some geometries.

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Heng Zhi Feng (zh...@hsr.ch)
 Sent: Wednesday, December 04, 2013 12:07 PM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] [postgis] Optimal tile size for Raster
 
 Hi Pierre,
 
 Thanks for the suggestion.
 
 But how could I get a geometry table in Postgres that is sure to intersect
 the raster image?
 The reason why I created another raster table is so that I am able to vary
 the query size. How can I do it with your suggestions (using polygons) since
 the polygon may not intersect with the raster.
 
 Please correct me if I am wrong, as I am still new to this concept.
 
 Regards,
 Zhi Feng
 
 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Pierre Racine
 Sent: Mittwoch, 4. Dezember 2013 17:43
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] [postgis] Optimal tile size for Raster
 
 You are getting fast result because you are not operating at the pixel level.
 Only at the raster extent level (converted to geometry). A more usual query
 would be to intersect a raster with some polygons:
 
 SELECT (ST_Intersection(r.rast, q.geom)).* FROM ilatlon32x32 as r,
 geometrytable as q WHERE ST_Intersects(r.rast, q.geom)
 
 Why do you want to intersect two rasters? Generally it is to mask one
 raster with another.
 
 Pierre
 
  -Original Message-
  From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
  boun...@lists.osgeo.org] On Behalf Of Heng Zhi Feng (zh...@hsr.ch)
  Sent: Wednesday, December 04, 2013 10:56 AM
  To: postgis-users@lists.osgeo.org
  Subject: [postgis-users] [postgis] Optimal tile size for Raster
 
  Hello,
 
 
 
  I am having to import Raster (.tiff) into Postgres using the raster2pgsql
 tool.
  The raster used is 12MB in size and the Postgres version is 9.1 on
  Ubuntu 13.10.
 
 
 
  The virtual machine has the following hardware specs.
 
  -  17408 Memory
 
  -  4 Processors
 
  -  16GB HDD
 
  The command used to import (with the parameters) is: raster2pgsql -s
  4326 -I -C -M -R -d -l 4 U:\Desktop\ilatlon_float.tif -F -t 8x8
  ilatlon8x8|psql - d Raster
 
  And the tile size varies from 8x8 all the way up till 2048 (power of 2s).
 
 
 
  My results from the query ran showed that with larger tile size, the
  query has a faster returned timing.
 
  This trend is however opposite to some of the similar tests conducted
  by others (from the internet), where their conclusion was smaller tile
  size will have better performance.
 
 
 
  This is my one of the query I ran:
 
  EXPLAIN (ANALYZE,BUFFERS)
 
  SELECT DISTINCT ST_Intersection(r.rast::geometry,q.rast::geometry)
 
  FROM ilatlon32x32 as r, boundedilatlon32x32 as q;
 
 
 
  Is someone able to give advice on this, if this result that I obtained
  was legit? Or suggest where could the problem be that gave this result?
 
 
 
  Regards,
 
  Zhi Feng
 
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] Postgis raster : jfif? jp2?

[Pierre Racine]

Laurent,

The list of available drivers is dependent on your installation. Normally if 
you can gdalinfo a raster, raster2psgql should be able to load it. So the first 
thing is: Does gdalinfo works fine? If not then you should refer to the gdal 
list. I yes then it might be a raster2pgsql problem which would be surprising.

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of laurent
 Sent: Monday, November 25, 2013 9:23 AM
 To: postgis-users@lists.osgeo.org
 Subject: [postgis-users] Postgis raster : jfif? jp2?
 
 Dear all,
 
 I would like import a world raster DEM in my postgis db. I would like use
 raster2pgsql.
 When i try to convert the jp2 file or jfif file to the .sql, i have an error
 message saying Unable to read raster.
 
 When i type: raster2pgsql -G , the jpeg jfif seems to be in the list. Why i
 have this error message? Except tiff, what raster format can be imported
 into postgis db?
 
 
 
 --
 View this message in context: http://postgis.17.x6.nabble.com/Postgis-
 raster-jfif-jp2-tp5004960.html
 Sent from the PostGIS - User mailing list archive at Nabble.com.
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

Not really the same thing as I said... ;-)

I would prefer You just DO WHAT THE FUCK YOU WANT TO as long as it stay 
open...

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Mateusz Loskot
 Sent: Wednesday, November 20, 2013 9:52 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] The first release of the PostGIS Add-ons is out!
 
  I would have gone with a license saying You can modify and redistribute
 as long as the derived work is also under an open license, not necessarily
 GPL. Does that make sense? Does that exist?
 
 http://en.wikipedia.org/wiki/WTFPL
 
 Best regards,
 --
 Mateusz  Loskot, http://mateusz.loskot.net
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] ST_SetValue problem

[Pierre Racine]

 The result is that I change raster values in only 1 of the 81 polygon cells 
 that I have in
 schema.data_vector. What am I missing here? Thanks for your help! Keep
 up the great work!

UPDATE schema.data_raster SET rast = ST_SetValue(rast,1, geom, 100)
FROM schema.data_vector

I don't fully understand the behavior of UPDATE in that case but, or it join 
with only one row of schema.data_vector, or it join with all of them but only 
the last one is used to update schema.data_raster. Or actually probably 
schema.data_raster is updated many times but always from the original value. 
Not incrementally.

Hugues suggestion:

SELECT rid,  ST_SetValue(rast,1, geom, 100) AS RAST
FROM schema.data_raster, schema.data_vector

Would produce as many raster rows as you have geometries...

There are two better ways to do this...

The first one is to ST_Union the rasterizization of the geometries and to union 
it with the existing raster:

SELECT ST_Union(rast, 'LAST') rast
FROM 
(SELECT rast FROM schema.data_raster
UNION ALL
SELECT ST_Union(ST_AsRaster(geom, rast, '16BSI', elevation, -1)) rast
FROM schema.data_raster, schema.data_vector
) foo

You can probably assign this single result as an UPDATE to your table..

The second one is to install the PostGIS Add-ons from 
https://github.com/pedrogit/postgisaddons/releases
and to do something like:

SELECT ST_ExtractToRaster(ST_AddBand(ST_MakeEmptyRaster(rast), '16BSI'), 
'schema', 'geomtable', 'data_vector', 'elevation', 
'MEAN_OF_VALUES_AT_PIXEL_CENTROID') rast FROM schema.data_raster;

You can probably also assign this single result as an UPDATE to your table...

Pierre

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

I knew we would fall into a license nightmare... So far we've got suggestions 
for No rights reserved, BSD, Creative Common Zero, GPL, LGPL, MIT...

Any of these license would make it hard to move a function from the Add-ons to 
PostGIS core? That would be a good candidate for elimination.

Why people do not like the PostGIS GPL?

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Mathieu Basille
 Sent: Monday, November 18, 2013 2:04 PM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] The first release of the PostGIS Add-ons is out!
 
 Le 11/18/2013 12:25 PM, Stephen Woodbridge a écrit :
  On 11/18/2013 12:20 PM, Mateusz Loskot wrote:
  On 18 November 2013 17:13, Pierre Racine
 pierre.rac...@sbf.ulaval.ca
  wrote:
  I'm a license ignorant.
 
  Release it into the public domain and
  include a statement that you release all rights
  adding No rights reserved.
 
  Best regards,
 
+1 on this because it is the most friendly and can be used by everyone
  regardless of license that they are using. It absolutely has no license
  conflicts.
 
 Although I'm a very big fan of public domain, it has its own limitations.
 For instance, some countries (randomly picked: France) do not allow one to
 declare its creation in the public domain. This is something that's granted
 from its nature (e.g. a representative speech or a math formula) or that is
 gained after a couple of decades after the death of the author (which is
 quite unlikely for PostGIS).
 
 Instead, I would much favor explicit licenses, such as the much simple
 (2-clauses) BSD [1] or the Creative Common Zero [2]. Both of them give the
 user the maximum flexibility and make sure there is no license conflict
 afterwards. The BSD, my favorite, is also super easy to read and understand.
 
 Sincerely,
 Mathieu.
 
 
 [1] http://opensource.org/licenses/BSD-2-Clause
 
 [2] https://creativecommons.org/publicdomain/zero/1.0/
 
 
 
  Best,
 -Steve
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 --
 
 ~$ 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
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

I guess the fact that PL/pgSQL is not compiled and hence is open by nature 
should have an incidence on the choice of the license. No?

I bend toward GPL in any case...

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Pierre Racine
 Sent: Tuesday, November 19, 2013 1:45 PM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] The first release of the PostGIS Add-ons is out!
 
 I knew we would fall into a license nightmare... So far we've got suggestions
 for No rights reserved, BSD, Creative Common Zero, GPL, LGPL, MIT...
 
 Any of these license would make it hard to move a function from the Add-
 ons to PostGIS core? That would be a good candidate for elimination.
 
 Why people do not like the PostGIS GPL?
 
 Pierre
 
  -Original Message-
  From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
  boun...@lists.osgeo.org] On Behalf Of Mathieu Basille
  Sent: Monday, November 18, 2013 2:04 PM
  To: PostGIS Users Discussion
  Subject: Re: [postgis-users] The first release of the PostGIS Add-ons is 
  out!
 
  Le 11/18/2013 12:25 PM, Stephen Woodbridge a écrit :
   On 11/18/2013 12:20 PM, Mateusz Loskot wrote:
   On 18 November 2013 17:13, Pierre Racine
  pierre.rac...@sbf.ulaval.ca
   wrote:
   I'm a license ignorant.
  
   Release it into the public domain and
   include a statement that you release all rights
   adding No rights reserved.
  
   Best regards,
  
 +1 on this because it is the most friendly and can be used by everyone
   regardless of license that they are using. It absolutely has no license
   conflicts.
 
  Although I'm a very big fan of public domain, it has its own limitations.
  For instance, some countries (randomly picked: France) do not allow one
 to
  declare its creation in the public domain. This is something that's granted
  from its nature (e.g. a representative speech or a math formula) or that is
  gained after a couple of decades after the death of the author (which is
  quite unlikely for PostGIS).
 
  Instead, I would much favor explicit licenses, such as the much simple
  (2-clauses) BSD [1] or the Creative Common Zero [2]. Both of them give
 the
  user the maximum flexibility and make sure there is no license conflict
  afterwards. The BSD, my favorite, is also super easy to read and
 understand.
 
  Sincerely,
  Mathieu.
 
 
  [1] http://opensource.org/licenses/BSD-2-Clause
 
  [2] https://creativecommons.org/publicdomain/zero/1.0/
 
 
  
   Best,
  -Steve
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
  --
 
  ~$ 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
  ___
  postgis-users mailing list
  postgis-users@lists.osgeo.org
  http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 ___
 postgis-users mailing list
 postgis-users@lists.osgeo.org
 http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

 IMHO, It's fine (actually wonderful) that postgis is in the GPL, it protects 
 the
 project and does not really limit its use in any substantive way.  It's big
 enough to stand on its own.  But re: #2 GPL would limit the wide use of
 these scripts.

How would GPL limit the wide use of these scripts? I find it's ok that any 
derivative be licensed the same way. It doesn't force a using software to be 
completely under GPL, just the derivative scripts of those scripts. It seems to 
me that the fact that PL/pgSQL is open by nature makes it very easy to 
redistribute these files or any derivative with the same license without 
limiting the rest of the software to be under any other kind of license. Does 
that make sense?

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

[postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

Hey!

I'm glad today to launch the PostGIS Add-ons. This new PostGIS extension aims 
at providing a quick, more Agile way to release pure PL/pgSQL functions 
developed by the PostGIS user community. It includes 15 new functions. Among 
them, two functions to aggregate raster/vector overlay statistics, one to 
extract values from a vector coverage to a raster, one to generate random 
points inside a polygon and two new aggregates functions helping removing 
overlaps in a table of geometry.

Read more about the PostGIS Add-ons at 
http://geospatialelucubrations.blogspot.ca/2013/11/launching-postgis-add-ons-new-postgis.html

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

Good point. I'm such not a license guy Should be the same as PostGIS (GPL). 
I will add license notice in the next release (very soon!).

Pierre

 -Original Message-
 From: postgis-users-boun...@lists.osgeo.org [mailto:postgis-users-
 boun...@lists.osgeo.org] On Behalf Of Rémi Cura
 Sent: Monday, November 18, 2013 11:49 AM
 To: PostGIS Users Discussion
 Subject: Re: [postgis-users] The first release of the PostGIS Add-ons is out!
 
 Very cool idea.
 I have a lot of plpgsql functions.
 
 One more detail : what about licenses?
 
 Cheers,
 Rémi-C
 
 
 2013/11/18 Pierre Racine pierre.rac...@sbf.ulaval.ca
 
 
   Hey!
 
   I'm glad today to launch the PostGIS Add-ons. This new PostGIS
 extension aims at providing a quick, more Agile way to release pure
 PL/pgSQL functions developed by the PostGIS user community. It includes
 15 new functions. Among them, two functions to aggregate raster/vector
 overlay statistics, one to extract values from a vector coverage to a raster,
 one to generate random points inside a polygon and two new aggregates
 functions helping removing overlaps in a table of geometry.
 
   Read more about the PostGIS Add-ons at
 http://geospatialelucubrations.blogspot.ca/2013/11/launching-postgis-
 add-ons-new-postgis.html
 
   Pierre
   ___
   postgis-users mailing list
   postgis-users@lists.osgeo.org
   http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
 
 

___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

Humm.. Lots of requirements...

 Sorry for the license, I needed a clear answer :-) I personnaly would prefere
 LGPL, like postgres.

Doesn't seems like a good idea to make it different than PostGIS. How is LGPL 
different from LPG? I'm a license ignorant.
 
 can you consider to use a template doc (wiki in github?, or pure markdown)
 , more detailled than just in code comment? I feel the current PostGIS doc
 per function should be a lower limit of documentation.

There is a quick list of available functions at the beginning of the file. I 
would rather go for less maintenance as possible. If you want to copy/remove 
the doc from the file to the wiki. Feel free to do it.

 Also important : you may quickely go to several dozens of functions, why
 not decide now for an architecture (classifying the contrib), maybe with
 label, maybe with classical folders? Like raster, util, rewrite, etc 
 etc.

A small step at a time...
 
 Last : what is the process for a contrib to be accepted? Surely it should be
 tested by another personn than dev before being added?

I wrote some criterias at the beginning of the file. It goes like this:

-your function is written in pure PL/pgSQL (no C!),
-your function is generic enough to be useful to other PostGIS users, 
-you follow functions and variables naming and indentation conventions already 
in use in the files, 
-you document your function like the ones already provided, 
-you provide some tests in the postgis_addons_test.sql file,
-you add the necessary DROP statements in the postgis_addons_uninstall.sql file.

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users

Re: [postgis-users] The first release of the PostGIS Add-ons is out!

[Pierre Racine]

Rémi,

 About license :
 if you use GPL function, your code becomes GPL, this might be a problem if
 using it inside company.
 If you use LGPL, you can do basically what you want.

I think about going like Mateusz and Stephen suggested. What do you think?

 Using PostGIS everyday, I already find it hard to find the functions I need. I
 feel there is no point to provide additionnal functions if you have to know
 them to use them.

Every function is preceded with a description of the parameters and a 
self-contained and a typical example. The question is: Should this 
documentation be: 1) embedded in the .sql file 2) in the readme 3) in the wiki? 
for now I thought that embedding it in the file, beside the code, with a quick 
list at the beginning like the one in the PostGIS doc, was the most simple way.

 Say I'm developping and I need a function, how would I know you have a
 function I might use, by searching trough all the functions name?

You look in the quick list at the beginning.

 I understand your maintenance concern, yet it could be mandatory to
 document new function in markdown or in wiki when someone contribute,
 thus eanbling a google search.

It could be. I tried to put as less barrier as possible to new contributors.

 About test : Your idea being more for users than regular developpers, how
 to ensure quality.

There is a .sql test file to be executed manually. Please give it a try.

 I may be totally wrong, but I fear it will become like many open source
 graveyard of non-maintained poorly written functions (I include myself in
 people whose code have to be reviewed).

I engage to be the guardian of quality for a while.  Hopefully you're wrong and 
I suggest you be the first contributor... ;-)

For myself I had most of those functions written since a while and they were 
spread in many files. I had to clean them so they are easier to find and 
install for my users.

Pierre
___
postgis-users mailing list
postgis-users@lists.osgeo.org
http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users