[postgis-users] point to raster conversion

Thu Dec 19 07:26:33 PST 2013

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.cura at 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.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/1b9300a0/attachment.html>