[postgis-users] ST_Buffer + grid problem

[Pierre Racine]

The raster was just a quick way to make a vector index grid. If you can't get ST_AsRaster to work by upgrading to GDAL 1.9 use ST_MakeEmptyRaster and then ST_MapalgebraExpr with the '([rast.x] - 1) * '  || ST_Width(rast)::text || ' + [rast.y]' expression.

Then convert each pixel to a table of polygons with (ST_PixelAsPolygons(rast)).geom  & (ST_PixelAsPolygons(rast)).val

Then proceed with a normal geometry/geometry intersection.

We don't have a very fast way to do that in raster mode right now. The fastest approach with what we have would probably be to use ST_MapAlgebraFct and write a custom function that generate the values by doing a point/geometry instersection with the road buffer table.

Pierre

> -----Original Message-----
> From: postgis-users-bounces at postgis.refractions.net [mailto:postgis-users-
> bounces at postgis.refractions.net] On Behalf Of Ed Linde
> Sent: Thursday, March 22, 2012 5:03 PM
> To: PostGIS Users Discussion
> Subject: Re: [postgis-users] ST_Buffer + grid problem
> 
> Oh and if you think there is a smarter way to do this, like maybe not making a
> raster like I am trying to etc then please let me know. The idea is that the road
> geometry is way way smaller than the point cloud which is in Terabytes, so I am
> making a grid index in some way around the road buffers, to speed up the
> computation.
> 
> 
> On Thu, Mar 22, 2012 at 10:01 PM, Ed Linde <edolinde at gmail.com> wrote:
> 
> 
> 	Hi Pierre,
> 	The idea is that I want to build this uniform grid on the road geometry
> and like you mention it has dimensions 1020x24798. I want to assign each cell a
> unique cell ID based on some formula like you had earlier ...  x*width + y to get a
> 1D cell ID from (x,y). Hopefully I can do this on just (lat.long) of srid 4326
> directly. Then I need to
> 	intersect it to the buffers around the roads which I think should be easy
> once I have a raster, so this is just a normal intersect. Now I should know each
> cell and what road buffer it belongs to! Then on my second data set which is a
> massive point cloud with srid = 4326 as well, I just want to pass through it ONCE
> and compute in which cell it falls and hence associate a point with a road buffer.
> Hope that made sense :)
> 
> 	Cheers,
> 	Ed
> 
> 	On Thu, Mar 22, 2012 at 9:53 PM, Pierre Racine
> <Pierre.Racine at sbf.ulaval.ca> wrote:
> 
> 
> 		1020x24798 so more than 25 000 000 polygons... or pixels...
> 
> 		So if I understand well you want to assign some values to each
> of those cells based on a vector coverage (of how many polygons)?
> 
> 
> 		> -----Original Message-----
> 		> From: postgis-users-bounces at postgis.refractions.net
> [mailto:postgis-users-
> 		> bounces at postgis.refractions.net] On Behalf Of Ed Linde
> 
> 		> Sent: Thursday, March 22, 2012 4:43 PM
> 		> To: PostGIS Users Discussion
> 		> Subject: Re: [postgis-users] ST_Buffer + grid problem
> 		>
> 
> 		> So far no proper raster, because the following failed, but
> maybe the polygon
> 		> gives us an idea? This polygon is the extent of all the road
> geometries in
> 		> Denmark. Also I remember that when you did the raster at
> your end you said it
> 		> worked. Is 0.000036 degrees = 4m correct if I want to get a
> 4m by 4m cell sized
> 		> uniform raster grid?
> 		>
> 		> select ST_AsRaster(
> ST_GeomFromText('POLYGON((8.07734039737749
> 		> 54.4984986588244,8.07734039737749
> 57.7505109647578,15.1919565742587
> 		> 57.7505109647578,15.1919565742587
> 54.4984986588244,8.07734039737749
> 		> 54.4984986588244))'), 0.000036, 0.000036);
> 		>
> 		>
> 		> On Thu, Mar 22, 2012 at 9:39 PM, Pierre Racine
> <Pierre.Racine at sbf.ulaval.ca>
> 		> wrote:
> 		>
> 		>
> 		>       > Thanks. The raster I need to visualise is a 4m by 4m grid
> on the entire
> 		> map of
> 		>       > Denmark! :) So do you classify that as a large raster? If so
> is there a
> 		> way to see a
> 		>       > portion of it or something?
> 		>       > I just want to manually check for a few areas on the map
> that the
> 		> intersection
> 		>       > indeed works as expected and things haven't gone awry
> thanks to
> 		> projection
> 		>       > differences or that I had degrees instead of meters or
> some such
> 		> thing. Well, I
> 		>       > still have to first test if this GDAL upgrade will fix things
> and if make
> 		> empty raster
> 		>       > .. makes a difference.
> 		>
> 		>
> 		>       That must be big but Danemark is a small country ;-) Only
> width and
> 		> height tell you if a raster is big.
> 		>
> 		>       In the database, your grid is a set of rectangular geometry
> (how many?)
> 		> or just a big raster now (width & height)?
> 		>
> 		>       _______________________________________________
> 		>       postgis-users mailing list
> 		>       postgis-users at postgis.refractions.net
> 		>       http://postgis.refractions.net/mailman/listinfo/postgis-
> users
> 		>
> 		>
> 
> 		_______________________________________________
> 		postgis-users mailing list
> 		postgis-users at postgis.refractions.net
> 		http://postgis.refractions.net/mailman/listinfo/postgis-users
> 
> 
>