[postgis-users] point to raster conversion

[Rémi Cura]

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

[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

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

[Rémi Cura]

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 :

[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-us...@lists.osgeo.org [mailto:postgis-users-
> bou...@lists.osgeo.org] On Behalf Of Rémi Cura
> Sent: Thursday, December 19, 2013 8:52 AM
> To: PostGIS Users Discussion

[Rémi Cura]

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.03170000003 20680.0229999999 49.02441796875) 49024 8110

2 POINT Z (2247.76560000001 20680.2519999999 49.0198828125) 49020 7890

3 POINT Z (2248.25289999996 20680.2479999999 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

WHERE gid = 87263

)

,point AS (

SELECT  PC_Explode(patch) AS point

FROM patch

)

SELECT row_number() over() AS qgis_id, ST_SetSRID(point::geometry,932011) AS geom, @(PC_Get(point, 'z')*1000)::int AS z, @((PC_Get(point, 'reflectance')+50)*200)::int aS reflectance

FROM point;

--2746 points

--checking

SELECT *

FROM one_patch_into_points

--1 01010000A0AB380E00DC7EFB3A1090A140CDFDD4780132D440A69BC42020834840 49024 8110

--2 POINT Z (2247.76560000001 20680.2519999999 49.0198828125) 49020 7890

--updating the raster to set it on the patch place.

--we get the upper left coordinate of patch

WITH patch  AS (

SELECT gid,patch, ST_AsText(patch::geometry)

FROM acquisition_tmob_012013.riegl_pcpatch_space

WHERE gid = 87263

)

,centroid AS (

SELECT ST_Centroid(patch::geometry) AS centroid

FROM patch

)

,upperleft AS (

SELECT round(ST_X(centroid)) -0.5 AS upper_left_x, round(ST_y(centroid))-0.5 AS upper_left_y

FROM centroid

)--updating the raster with it

UPDATE patch_to_raster SET rast =  ST_SetUpperLeft(rast,upper_left_x,upper_left_y)

FROM upperleft;

--checking the raster

SELECT ST_Summary(rast) FROM patch_to_raster;

--Raster of 50x50 pixels has 2 bands and extent of BOX(2247.5 20679.5,2248.5 20680.5) 

--note : the patch bbox : POLYGON((2247.51 20679.50,2247.51 20680.49,2248.49 20680.49,2248.49 20679.50,2247.51 20679.50)), some pixels will be empty

--updating the raster with pixel value

--using postgis_addons function : https://raw.github.com/pedrogit/postgisaddons/master/postgis_addons.sql

--

UPDATE patch_to_raster SET rast = 

ST_ExtractToRaster(

rast 

,band:=1

,schemaname:= 'test_raster'

,tablename:= 'one_patch_into_points'

,geomrastcolumnname:= 'geom'

,valuecolumnname:='z'

);

--5 sec : way too long!

UPDATE patch_to_raster SET rast = 

ST_ExtractToRaster(

rast 

,band:=2

,schemaname:= 'test_raster'

,tablename:= 'one_patch_into_points'

,geomrastcolumnname:= 'geom'

,valuecolumnname:='reflectance'

,method:='MEAN_OF_VALUES_AT_PIXEL_CENTROID'

);

--6 sec, way too long !

--checking the value of pixel in both band (no shortcut to get several bands at a time??)

SELECT  s1 as x, s2 as y, ST_Value(rast, 1,s1, s2, exclude_nodata_value:=true) AS value_band_1, ST_Value(rast, 2,s1, s2, exclude_nodata_value:=true) AS value_band_2

FROM   generate_series(1,50) as s1,   generate_series(1,50) as s2, patch_to_raster;

--values are NULL!  

--checking the histogram

SELECT ST_Histogram( rast, 1 ) --, boolean exclude_nodata_value=true

FROM patch_to_raster;

--no line returned

--checking result 

SELECT ST_Summary(rast) FROM patch_to_raster;

--Raster of 50x50 pixels has 2 bands and extent of BOX(2247.5 20679.5,2248.5 20680.5)

--  band 1 of pixtype 32BF is in-db with NODATA value of -3.40282346638529e+38

--  band 2 of pixtype 32BF is in-db with NODATA value of -3.40282346638529e+38

////////////////////////////////

[Rémi Cura]

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...@gmail.com>

[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-us...@lists.osgeo.org [mailto:postgis-users-
> bou...@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
>
>

[Rémi Cura]

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...@sbf.ulaval.ca>

[Rémi Cura]

Buffering make it works.

Cheers,

Rémi-C

2013/12/19 Rémi Cura <remi...@gmail.com>

Hey, thanks for the answer, and nice catch !

[Rémi Cura]

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

[Rémi Cura]

QGis 2.0.1.

I tried as well to convert to UINT16

no sucess

2013/12/19 Rémi Cura <remi...@gmail.com>

[Rémi Cura]

For the record :

when outputting the image to disk then loading it into postgis, image shows fine (band 1 and band2)

When 'add canvas' on the image in QGIS, a nex layer is created, but empty.

Do you know other way to do it?

Cheers,

Rémi-C

[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...

[Pierre Racine]

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

> -----Original Message-----
> From: postgis-us...@lists.osgeo.org [mailto:postgis-users-
> bou...@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
>

[Jan Hartmann]

Hi all,

I have written a simple OpenLayers-application to edit a PostGIS polygon layer: read one polygon from the server as GeoJSON, adapt it interactively, and write it back again via GeoJSON. Points are only moved, not created or deleted.

The polygon layer is topologically correct: all adjoining polygons share the same points. With QGIS I can choose "topological editing", to automatically move corresponding points from adjoining polygons. With the GeoJSON application that is not possible: whenever I move a boundary point of one polygon, the corresponding points on adjoining polygons remain where they are.

Is it possible with the new topology methods in PostGIS to write the following update procedure: whenever a changed polygon is written back into the database, all corresponding points on neighbouring polygons will be looked up and changed as well? I guess I would need some sort of pointer for all boundary points to all other boundary points that share the same location.

Jan

[Rémi Cura]

Hey,

you seem to mix 2 different things.

_You can do what you want using the regular postgis function

You would write a trigger on your polygon column, os that on update, you also update other touching polygons (you would need to do polygon breanking into lines, which may become difficult if you have inner ring).

Sorry to say that this may not be a good idea to use this as a topological model (a lot of data duplication, non standart topology model)

_You can do what you want (with a tweak) with postgis topology

In postgis topology you have shared lines and points, then polygons are defined upon those (like : edge 23 then edge 43 then edge 46 form the polygon 1)

So out of the box, when modifying a shared line, all surfaces are implicitly modified.

The small problem you may have is that the surfaces are computed each time you need it, which is slower than just store it like a polygon.

If you need some precision about the data model behind postgis topology : http://trac.osgeo.org/postgis/wiki/PostgisTopology_Data_Model

Cheers,

Rémi-C

2014/1/3 Jan Hartmann <j.l.h.h...@uva.nl>

[Paul Ramsey]

Maybe on OL3?

http://vimeo.com/78298228

P

[Jan Hartmann]

Thanks Paul, that was indeed what I was looking for.

I'm not sure if it would be wise to switch already to OL3 at this moment. What would you suggest would be a reasonable time to wait for the definitive version for a production environment? And would it make sense to work in the meantime on my original idea, implementing this functionality within PostGIS-topology?

Jan

[Rémi Cura]

If you want to modify postgis topology maybe it could be in the way of a new column in the "face" table,

with explicit geometry.

it would be keep in sync with triggers.

You have to keep in mind that is would be duplicating data.

You can get example triggers in other topology table, and I think there already is a function to get a polygon  geom from a face id.

Cheers,

Rémi-C

2014/1/6 Jan Hartmann <j.l.h.h...@uva.nl>