<div dir="ltr">Oups,<div>I figured it would be of no use.</div><div><br></div><div>I of course works on pure geometry, raster functions only sees postgis points and ints.<br>Here few lines of the point table</div><div><br>
</div><div>qgis_id, st_astext(geom) AS geom,z,reflectance<br></div><div><div>1<span class="" style="white-space:pre"> </span>POINT Z (2248.03170000003 20680.0229999999 49.02441796875)<span class="" style="white-space:pre"> </span>49024<span class="" style="white-space:pre"> </span>8110</div>
<div>2<span class="" style="white-space:pre"> </span>POINT Z (2247.76560000001 20680.2519999999 49.0198828125)<span class="" style="white-space:pre"> </span>49020<span class="" style="white-space:pre"> </span>7890</div><div>
3<span class="" style="white-space:pre"> </span>POINT Z (2248.25289999996 20680.2479999999 49.02616015625)<span class="" style="white-space:pre"> </span>49026<span class="" style="white-space:pre"> </span>8404</div></div>
<div><br></div><div><br></div><div>Is there any function to work on pixel sets , (like read/write many at a time) </div><div><br></div><div>Now the code : </div><div>(I can make a test case out of this if you want)</div><div>
<br></div><div>Thanks,</div><div><br></div><div>Rémi-C</div><div><br></div><div>///////////////////////////////</div><div class="gmail_extra"><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div class=""><div class=""><br></div><div class="">SELECT postgis_full_version();</div><div class="">--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</div>
<div class=""><br></div><div class=""><br></div><div class=""><br></div><div class=""><br></div><div class="">----converting a patch to raster</div><div class=""><span class="" style="white-space:pre"> </span>--get patch extend</div>
<div class=""><span class="" style="white-space:pre"> </span>--create a raster with this extend</div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--explode patch to points</div><div class="">
<span class="" style="white-space:pre"> </span>--set pixel value to point intersecting it.</div><div class=""><br></div><div class=""><br></div><div class="">--preparing raster :</div><div class=""><span class="" style="white-space:pre"> </span>--creating the raster table</div>
<div class=""><span class="" style="white-space:pre"> </span>DROP TABLE IF EXISTS patch_to_raster;</div><div class=""><span class="" style="white-space:pre"> </span>CREATE TABLE patch_to_raster(rid serial primary key</div>
<div class=""><span class="" style="white-space:pre"> </span>, rast raster);</div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--</div><div class=""><span class="" style="white-space:pre"> </span>--DELETE from patch_to_raster</div>
<div class=""><span class="" style="white-space:pre"> </span>--inserting an empty raster</div><div class=""><span class="" style="white-space:pre"> </span>INSERT INTO patch_to_raster (rast) VALUES (</div><div class=""><span class="" style="white-space:pre"> </span>ST_MakeEmptyRaster(</div>
<div class=""><span class="" style="white-space:pre"> </span>50</div><div class=""><span class="" style="white-space:pre"> </span>,50</div><div class=""><span class="" style="white-space:pre"> </span>,upperleftx:=0</div>
<div class=""><span class="" style="white-space:pre"> </span>,upperlefty:=0</div><div class=""><span class="" style="white-space:pre"> </span>,scalex:=1</div><div class=""><span class="" style="white-space:pre"> </span>, scaley:=1</div>
<div class=""><span class="" style="white-space:pre"> </span>, skewx:=0</div><div class=""><span class="" style="white-space:pre"> </span>,skewy:=0</div><div class=""><span class="" style="white-space:pre"> </span>, srid:=932011</div>
<div class=""><span class="" style="white-space:pre"> </span>) --srid of translated lambert 93 to match laser referential</div><div class=""><span class="" style="white-space:pre"> </span>);</div><div class=""><span class="" style="white-space:pre"> </span></div>
<div class=""><span class="" style="white-space:pre"> </span>--setting the correct pixel size</div><div class=""><span class="" style="white-space:pre"> </span>UPDATE patch_to_raster SET rast = ST_SetScale(rast,0.02,0.02);</div>
<div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking the content:</div><div class=""><span class="" style="white-space:pre"> </span>SELECT ST_Summary(rast)</div><div class=""><span class="" style="white-space:pre"> </span>FROM patch_to_raster;</div>
<div class=""><span class="" style="white-space:pre"> </span>--Raster of 50x50 pixels has 0 bands and extent of BOX(0 0,1 1)</div><div class=""><br></div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--adding band </div>
<div class=""><span class="" style="white-space:pre"> </span>--first band is Z : float</div><div class=""><span class="" style="white-space:pre"> </span>--second band is reflectance : float</div><div class=""><span class="" style="white-space:pre"> </span>UPDATE patch_to_raster SET rast = ST_AddBand( --1'st band, Z</div>
<div class=""><span class="" style="white-space:pre"> </span>rast</div><div class=""><span class="" style="white-space:pre"> </span>,pixeltype:='32BF' -- '32BUI'</div><div class=""><span class="" style="white-space:pre"> </span>, initialvalue:=NULL</div>
<div class=""><span class="" style="white-space:pre"> </span>, nodataval:=NULL</div><div class=""><span class="" style="white-space:pre"> </span>);</div><div class=""><span class="" style="white-space:pre"> </span>UPDATE patch_to_raster SET rast = ST_AddBand( --2nd band, reflectance</div>
<div class=""><span class="" style="white-space:pre"> </span>rast</div><div class=""><span class="" style="white-space:pre"> </span>, pixeltype:='32BF' --'32BUI' </div><div class=""><span class="" style="white-space:pre"> </span>, initialvalue:=NULL</div>
<div class=""><span class="" style="white-space:pre"> </span>, nodataval:=NULL</div><div class=""><span class="" style="white-space:pre"> </span>);</div><div class=""><br></div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking the content:</div>
<div class=""><span class="" style="white-space:pre"> </span>SELECT ST_Summary(rast) FROM patch_to_raster;</div><div class=""><span class="" style="white-space:pre"> </span>--Raster of 50x50 pixels has 2 bands and extent of BOX(0 0,1 1)</div>
<div class=""><span class="" style="white-space:pre"> </span>--band 1 of pixtype 32BF is in-db with no NODATA value</div><div class=""><span class="" style="white-space:pre"> </span>--band 2 of pixtype 32BF is in-db with no NODATA value</div>
<div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--getting a patch and creating a table with it </div><div class=""><span class="" style="white-space:pre"> </span>DROP TABLE IF EXISTS one_patch_into_points;</div>
<div class=""><span class="" style="white-space:pre"> </span>CREATE TABLE one_patch_into_points AS </div><div class=""><span class="" style="white-space:pre"> </span>WITH patch AS (</div><div class=""><span class="" style="white-space:pre"> </span>SELECT gid,patch</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM acquisition_tmob_012013.riegl_pcpatch_space</div><div class=""><span class="" style="white-space:pre"> </span>WHERE gid = 87263</div><div class=""><span class="" style="white-space:pre"> </span>)</div>
<div class=""><span class="" style="white-space:pre"> </span>,point AS (</div><div class=""><span class="" style="white-space:pre"> </span>SELECT PC_Explode(patch) AS point</div><div class=""><span class="" style="white-space:pre"> </span>FROM patch</div>
<div class=""><span class="" style="white-space:pre"> </span>)</div><div class=""><span class="" style="white-space:pre"> </span>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</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM point;</div><div class=""><span class="" style="white-space:pre"> </span>--2746 points</div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking</div>
<div class=""><span class="" style="white-space:pre"> </span>SELECT *</div><div class=""><span class="" style="white-space:pre"> </span>FROM one_patch_into_points</div><div class=""><span class="" style="white-space:pre"> </span>--1<span class="" style="white-space:pre"> </span>01010000A0AB380E00DC7EFB3A1090A140CDFDD4780132D440A69BC42020834840<span class="" style="white-space:pre"> </span>49024<span class="" style="white-space:pre"> </span>8110</div>
<div class=""><span class="" style="white-space:pre"> </span>--2<span class="" style="white-space:pre"> </span>POINT Z (2247.76560000001 20680.2519999999 49.0198828125)<span class="" style="white-space:pre"> </span>49020<span class="" style="white-space:pre"> </span>7890</div>
<div class=""><span class="" style="white-space:pre"> </span></div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--updating the raster to set it on the patch place.</div><div class="">
<span class="" style="white-space:pre"> </span>--we get the upper left coordinate of patch</div><div class=""><span class="" style="white-space:pre"> </span>WITH patch AS (</div><div class=""><span class="" style="white-space:pre"> </span>SELECT gid,patch, ST_AsText(patch::geometry)</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM acquisition_tmob_012013.riegl_pcpatch_space</div><div class=""><span class="" style="white-space:pre"> </span>WHERE gid = 87263</div><div class=""><span class="" style="white-space:pre"> </span>)</div>
<div class=""><span class="" style="white-space:pre"> </span>,centroid AS (</div><div class=""><span class="" style="white-space:pre"> </span>SELECT ST_Centroid(patch::geometry) AS centroid</div><div class=""><span class="" style="white-space:pre"> </span>FROM patch</div>
<div class=""><span class="" style="white-space:pre"> </span>)</div><div class=""><span class="" style="white-space:pre"> </span>,upperleft AS (</div><div class=""><span class="" style="white-space:pre"> </span>SELECT round(ST_X(centroid)) -0.5 AS upper_left_x, round(ST_y(centroid))-0.5 AS upper_left_y</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM centroid</div><div class=""><span class="" style="white-space:pre"> </span>)--updating the raster with it</div><div class=""><span class="" style="white-space:pre"> </span>UPDATE patch_to_raster SET rast = ST_SetUpperLeft(rast,upper_left_x,upper_left_y)</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM upperleft;</div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking the raster</div><div class=""><span class="" style="white-space:pre"> </span>SELECT ST_Summary(rast) FROM patch_to_raster;</div>
<div class=""><span class="" style="white-space:pre"> </span>--Raster of 50x50 pixels has 2 bands and extent of BOX(2247.5 20679.5,2248.5 20680.5) </div><div class=""><span class="" style="white-space:pre"> </span>--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</div>
<div class=""><br></div><div class=""><br></div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--updating the raster with pixel value</div><div class=""><span class="" style="white-space:pre"> </span>--using postgis_addons function : <a href="https://raw.github.com/pedrogit/postgisaddons/master/postgis_addons.sql">https://raw.github.com/pedrogit/postgisaddons/master/postgis_addons.sql</a></div>
<div class=""><span class="" style="white-space:pre"> </span>--</div><div class=""><span class="" style="white-space:pre"> </span>UPDATE patch_to_raster SET rast = </div><div class=""><span class="" style="white-space:pre"> </span>ST_ExtractToRaster(</div>
<div class=""><span class="" style="white-space:pre"> </span>rast </div><div class=""><span class="" style="white-space:pre"> </span>,band:=1</div><div class=""><span class="" style="white-space:pre"> </span>,schemaname:= 'test_raster'</div>
<div class=""><span class="" style="white-space:pre"> </span>,tablename:= 'one_patch_into_points'</div><div class=""><span class="" style="white-space:pre"> </span>,geomrastcolumnname:= 'geom'</div><div class="">
<span class="" style="white-space:pre"> </span>,valuecolumnname:='z'</div><div class=""><span class="" style="white-space:pre"> </span>);</div><div class=""><span class="" style="white-space:pre"> </span>--5 sec : way too long!</div>
<div class=""><span class="" style="white-space:pre"> </span></div><div class=""><span class="" style="white-space:pre"> </span>UPDATE patch_to_raster SET rast = </div><div class=""><span class="" style="white-space:pre"> </span>ST_ExtractToRaster(</div>
<div class=""><span class="" style="white-space:pre"> </span>rast </div><div class=""><span class="" style="white-space:pre"> </span>,band:=2</div><div class=""><span class="" style="white-space:pre"> </span>,schemaname:= 'test_raster'</div>
<div class=""><span class="" style="white-space:pre"> </span>,tablename:= 'one_patch_into_points'</div><div class=""><span class="" style="white-space:pre"> </span>,geomrastcolumnname:= 'geom'</div><div class="">
<span class="" style="white-space:pre"> </span>,valuecolumnname:='reflectance'</div><div class=""><span class="" style="white-space:pre"> </span>,method:='MEAN_OF_VALUES_AT_PIXEL_CENTROID'</div><div class="">
<span class="" style="white-space:pre"> </span>);</div><div class=""><span class="" style="white-space:pre"> </span>--6 sec, way too long !</div><div class=""><br></div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking the value of pixel in both band (no shortcut to get several bands at a time??)</div>
<div class=""><span class="" style="white-space:pre"> </span>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</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM generate_series(1,50) as s1, generate_series(1,50) as s2, patch_to_raster;</div><div class=""><span class="" style="white-space:pre"> </span>--values are NULL! </div>
<div class=""><span class="" style="white-space:pre"> </span></div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking the histogram</div><div class=""><span class="" style="white-space:pre"> </span>SELECT ST_Histogram( rast, 1 ) --, boolean exclude_nodata_value=true</div>
<div class=""><span class="" style="white-space:pre"> </span>FROM patch_to_raster;</div><div class=""><span class="" style="white-space:pre"> </span>--no line returned</div><div class=""><br></div><div class=""><span class="" style="white-space:pre"> </span>--checking result </div>
<div class=""><span class="" style="white-space:pre"> </span>SELECT ST_Summary(rast) FROM patch_to_raster;</div><div class=""><span class="" style="white-space:pre"> </span>--Raster of 50x50 pixels has 2 bands and extent of BOX(2247.5 20679.5,2248.5 20680.5)</div>
<div class=""><span class="" style="white-space:pre"> </span>-- band 1 of pixtype 32BF is in-db with NODATA value of -3.40282346638529e+38</div><div class=""><span class="" style="white-space:pre"> </span>-- band 2 of pixtype 32BF is in-db with NODATA value of -3.40282346638529e+38</div>
<div class="">////////////////////////////////</div></div></blockquote></div><br></div></div>