[postgis-users] point to raster conversion
Pierre Racine
Pierre.Racine at sbf.ulaval.ca
Thu Dec 19 08:15:10 PST 2013
"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-bounces at lists.osgeo.org [mailto:postgis-users-
> bounces at 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.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
> 01010000A0AB380E00DC7EFB3A1090A140CDFDD4780132D440A69BC420
> 20834840 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.s
> ql
> --
> 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
> ////////////////////////////////
>
>
More information about the postgis-users
mailing list