Hi Pierre,<br>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.<br>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.<br>
Can you please let me know how I can do this? So I did the make empty raster -> map algebra -> pixel as polygons approach! <br>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?<br>
MapAlgebra needed 2 rasters, so I inputted the same raster twice. So please help! :)<br><br><br>/*<br>CREATE OR REPLACE FUNCTION makegrid(geometry, integer)<br>RETURNS geometry AS<br>'SELECT ST_Collect(ST_SetSRID(ST_POINT(x,y),ST_SRID($1))) FROM <br>
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x<br>,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y <br>where st_intersects($1,ST_SetSRID(ST_POINT(x,y),ST_SRID($1)))'<br>
LANGUAGE sql<br>*/<br><br><br>--SELECT (md).*<br>--FROM (SELECT ST_MetaData(rast) As md <br>-- FROM (select ST_MakeEmptyRaster( 8, 4, 8.07734039737749, 57.7505109647578, 0.000036,0.000036,0,0, 4326 ) rast) foo) foo2;<br>
<br><br>create table vector_grid as <br>select (ST_PixelAsPolygons(rastfin)).geom , (ST_PixelAsPolygons(rastfin)).val <br>from<br>(<br>select ST_MapAlgebraExpr(rast, rast, '([rast.x] - 1) * ' || ST_Width(rast)::text || ' + [rast.y]' ) rastfin<br>
from (select ST_MakeEmptyRaster( 8, 4, 8.07734039737749, 57.7505109647578, 0.000036,0.000036,0,0, 4326 ) rast) foo) foo2;<br><br><b>NOTICE: The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent<br>
NOTICE: Raster do not have band 1. Returning null values<br>NOTICE: The two rasters provided do not have the respectively specified band indices. Returning no band raster of correct extent<br>NOTICE: Raster do not have band 1. Returning null values<br>
</b>Query returned successfully: 32 rows affected, 272 ms execution time.<br><br><br>select count(*) from vector_grid;<br>> 32 <br><br><div class="gmail_quote">On Thu, Mar 22, 2012 at 10:17 PM, Pierre Racine <span dir="ltr"><<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">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.<br>
<br>
Then convert each pixel to a table of polygons with (ST_PixelAsPolygons(rast)).geom & (ST_PixelAsPolygons(rast)).val<br>
<br>
Then proceed with a normal geometry/geometry intersection.<br>
<br>
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.<br>
<span class="HOEnZb"><font color="#888888"><br>
Pierre<br>
</font></span><div class="im HOEnZb"><br>
> -----Original Message-----<br>
> From: <a href="mailto:postgis-users-bounces@postgis.refractions.net">postgis-users-bounces@postgis.refractions.net</a> [mailto:<a href="mailto:postgis-users-">postgis-users-</a><br>
> <a href="mailto:bounces@postgis.refractions.net">bounces@postgis.refractions.net</a>] On Behalf Of Ed Linde<br>
</div><div class="HOEnZb"><div class="h5">> Sent: Thursday, March 22, 2012 5:03 PM<br>
> To: PostGIS Users Discussion<br>
> Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
><br>
> Oh and if you think there is a smarter way to do this, like maybe not making a<br>
> raster like I am trying to etc then please let me know. The idea is that the road<br>
> geometry is way way smaller than the point cloud which is in Terabytes, so I am<br>
> making a grid index in some way around the road buffers, to speed up the<br>
> computation.<br>
><br>
><br>
> On Thu, Mar 22, 2012 at 10:01 PM, Ed Linde <<a href="mailto:edolinde@gmail.com">edolinde@gmail.com</a>> wrote:<br>
><br>
><br>
> Hi Pierre,<br>
> The idea is that I want to build this uniform grid on the road geometry<br>
> and like you mention it has dimensions 1020x24798. I want to assign each cell a<br>
> unique cell ID based on some formula like you had earlier ... x*width + y to get a<br>
> 1D cell ID from (x,y). Hopefully I can do this on just (lat.long) of srid 4326<br>
> directly. Then I need to<br>
> intersect it to the buffers around the roads which I think should be easy<br>
> once I have a raster, so this is just a normal intersect. Now I should know each<br>
> cell and what road buffer it belongs to! Then on my second data set which is a<br>
> massive point cloud with srid = 4326 as well, I just want to pass through it ONCE<br>
> and compute in which cell it falls and hence associate a point with a road buffer.<br>
> Hope that made sense :)<br>
><br>
> Cheers,<br>
> Ed<br>
><br>
> On Thu, Mar 22, 2012 at 9:53 PM, Pierre Racine<br>
> <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>> wrote:<br>
><br>
><br>
> 1020x24798 so more than 25 000 000 polygons... or pixels...<br>
><br>
> So if I understand well you want to assign some values to each<br>
> of those cells based on a vector coverage (of how many polygons)?<br>
><br>
><br>
> > -----Original Message-----<br>
> > From: <a href="mailto:postgis-users-bounces@postgis.refractions.net">postgis-users-bounces@postgis.refractions.net</a><br>
> [mailto:<a href="mailto:postgis-users-">postgis-users-</a><br>
> > <a href="mailto:bounces@postgis.refractions.net">bounces@postgis.refractions.net</a>] On Behalf Of Ed Linde<br>
><br>
> > Sent: Thursday, March 22, 2012 4:43 PM<br>
> > To: PostGIS Users Discussion<br>
> > Subject: Re: [postgis-users] ST_Buffer + grid problem<br>
> ><br>
><br>
> > So far no proper raster, because the following failed, but<br>
> maybe the polygon<br>
> > gives us an idea? This polygon is the extent of all the road<br>
> geometries in<br>
> > Denmark. Also I remember that when you did the raster at<br>
> your end you said it<br>
> > worked. Is 0.000036 degrees = 4m correct if I want to get a<br>
> 4m by 4m cell sized<br>
> > uniform raster grid?<br>
> ><br>
> > select ST_AsRaster(<br>
> ST_GeomFromText('POLYGON((8.07734039737749<br>
> > 54.4984986588244,8.07734039737749<br>
> 57.7505109647578,15.1919565742587<br>
> > 57.7505109647578,15.1919565742587<br>
> 54.4984986588244,8.07734039737749<br>
> > 54.4984986588244))'), 0.000036, 0.000036);<br>
> ><br>
> ><br>
> > On Thu, Mar 22, 2012 at 9:39 PM, Pierre Racine<br>
> <<a href="mailto:Pierre.Racine@sbf.ulaval.ca">Pierre.Racine@sbf.ulaval.ca</a>><br>
> > wrote:<br>
> ><br>
> ><br>
> > > Thanks. The raster I need to visualise is a 4m by 4m grid<br>
> on the entire<br>
> > map of<br>
> > > Denmark! :) So do you classify that as a large raster? If so<br>
> is there a<br>
> > way to see a<br>
> > > portion of it or something?<br>
> > > I just want to manually check for a few areas on the map<br>
> that the<br>
> > intersection<br>
> > > indeed works as expected and things haven't gone awry<br>
> thanks to<br>
> > projection<br>
> > > differences or that I had degrees instead of meters or<br>
> some such<br>
> > thing. Well, I<br>
> > > still have to first test if this GDAL upgrade will fix things<br>
> and if make<br>
> > empty raster<br>
> > > .. makes a difference.<br>
> ><br>
> ><br>
> > That must be big but Danemark is a small country ;-) Only<br>
> width and<br>
> > height tell you if a raster is big.<br>
> ><br>
> > In the database, your grid is a set of rectangular geometry<br>
> (how many?)<br>
> > or just a big raster now (width & height)?<br>
> ><br>
> > _______________________________________________<br>
> > postgis-users mailing list<br>
> > <a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
> > <a href="http://postgis.refractions.net/mailman/listinfo/postgis-" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-</a><br>
> users<br>
> ><br>
> ><br>
><br>
> _______________________________________________<br>
> postgis-users mailing list<br>
> <a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
> <a href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
><br>
><br>
><br>
<br>
_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@postgis.refractions.net">postgis-users@postgis.refractions.net</a><br>
<a href="http://postgis.refractions.net/mailman/listinfo/postgis-users" target="_blank">http://postgis.refractions.net/mailman/listinfo/postgis-users</a><br>
</div></div></blockquote></div><br>