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.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 >> //////////////////////////////// >> > >
_______________________________________________ postgis-users mailing list postgis-users@lists.osgeo.org http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users