Re: [postgis-users] point to raster conversion

2013-12-19 Thread 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

2013-12-19 Thread 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.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
 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 0101A0AB380E00DC7EFB3A1090A140CDFDD4780132D440A69BC42020834840
 49024 8110
 --2 POINT Z (2247.7656001 20680.251999 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 )

Re: [postgis-users] point to raster conversion

2013-12-19 Thread 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.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
 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 0101A0AB380E00DC7EFB3A1090A140CDFDD4780132D440A69BC42020834840
 49024 8110
  --2 POINT Z (2247.7656001 20680.251999 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,

Re: [postgis-users] point to raster conversion

2013-12-19 Thread 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

2013-12-19 Thread 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.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

2013-12-19 Thread Rémi Cura 
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

2013-12-19 Thread 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


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

 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

2013-12-19 Thread 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

2013-12-19 Thread 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

2013-12-18 Thread 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