[postgis-users] point to raster conversion

[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
> ////////////////////////////////
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20131219/c817adb3/attachment.html>