[postgis-users] ST_Buffer + grid problem
Ed Linde
edolinde at gmail.com
Fri Mar 23 03:57:03 PDT 2012
Hi Pierre,
I tried the following based on your advice last night. I was unable to
upgrade to GDAL 1.9, so I went down the empty raster approach.
Below is another piece of code that uses "generate_series" to make a grid
but it doesn't work in my case as I have a float cell dimension i.e.
0.000036.
Can you please let me know how I can do this? So I did the make empty
raster -> map algebra -> pixel as polygons approach!
Since there is no dependency on data now, is it possible for you to maybe
outline the steps you performed at your end to get this working?
MapAlgebra needed 2 rasters, so I inputted the same raster twice. So please
help! :)
/*
CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_SetSRID(ST_POINT(x,y),ST_SRID($1))) FROM
generate_series(floor(st_xmin($1))::int,
ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x
,generate_series(floor(st_ymin($1))::int,
ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y
where st_intersects($1,ST_SetSRID(ST_POINT(x,y),ST_SRID($1)))'
LANGUAGE sql
*/
--SELECT (md).*
--FROM (SELECT ST_MetaData(rast) As md
-- FROM (select ST_MakeEmptyRaster( 8, 4, 8.07734039737749,
57.7505109647578, 0.000036,0.000036,0,0, 4326 ) rast) foo) foo2;
create table vector_grid as
select (ST_PixelAsPolygons(rastfin)).geom ,
(ST_PixelAsPolygons(rastfin)).val
from
(
select ST_MapAlgebraExpr(rast, rast, '([rast.x] - 1) * ' ||
ST_Width(rast)::text || ' + [rast.y]' ) rastfin
from (select ST_MakeEmptyRaster( 8, 4, 8.07734039737749, 57.7505109647578,
0.000036,0.000036,0,0, 4326 ) rast) foo) foo2;
*NOTICE: The two rasters provided do not have the respectively specified
band indices. Returning no band raster of correct extent
NOTICE: Raster do not have band 1. Returning null values
NOTICE: The two rasters provided do not have the respectively specified
band indices. Returning no band raster of correct extent
NOTICE: Raster do not have band 1. Returning null values
*Query returned successfully: 32 rows affected, 272 ms execution time.
select count(*) from vector_grid;
> 32
On Thu, Mar 22, 2012 at 10:17 PM, Pierre Racine <Pierre.Racine at sbf.ulaval.ca
> wrote:
> 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
> >
> >
> >
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at postgis.refractions.net
> http://postgis.refractions.net/mailman/listinfo/postgis-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20120323/550679b5/attachment.html>
More information about the postgis-users
mailing list