<div dir="ltr"><div dir="ltr"><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, 9 Nov 2024 at 11:52, Regina Obe <<a href="mailto:lr@pcorp.us">lr@pcorp.us</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg-8107322994282190670"><div lang="EN-US" style="overflow-wrap: break-word;"><div class="m_-8107322994282190670WordSection1"><p class="MsoNormal"><span style="font-size:11pt">The main overhead I see with postgis raster are the following<u></u><u></u></span></p><p class="MsoNormal"><span style="font-size:11pt"><u></u> <u></u></span></p><ol style="margin-top:0in" start="1" type="1"><li class="m_-8107322994282190670MsoListParagraph" style="margin-left:0in"><span style="font-size:11pt">Limitation of raster sizes of about 1GB as I recall<u></u><u></u></span></li><li class="m_-8107322994282190670MsoListParagraph" style="margin-left:0in"><span style="font-size:11pt">Detoasting of data which is almost always needed  <a href="https://www.postgresql.org/docs/current/storage-toast.html" target="_blank">https://www.postgresql.org/docs/current/storage-toast.html</a><u></u><u></u></span></li><li class="m_-8107322994282190670MsoListParagraph" style="margin-left:0in"><span style="font-size:11pt">Mem copying of data, probably more so than using say GDAL directly<u></u><u></u></span></li><li class="m_-8107322994282190670MsoListParagraph" style="margin-left:0in"><span style="font-size:11pt">Probably easier to do parallel operations on command line (though admittedly I haven’t tried)</span></li></ol></div></div></div></blockquote><div><br></div><div>In a recent experience where I have to join vector data with lot of raster data (stored on S3, with resolutions from 0.25 to 0.0027°, some of them about 200GB), I found the most efficient way was to keep rasters as files (COG, ZSTD-compressed) and to perform the raster query using gdallocationinfo (that can process millions of coords in seconds (last run: COPY 1 333 259  3.983s total).</div><div><br></div><div>The process consist of computing the pixels under each vector geometry (and only under geoms, with st_squareGrid and st_snapToGrid for polygons, and just st_segmentize and st_snapToGrid for lines) using raster metadata (resolution, origin), and to feed these pixel centroids to gdallocationinfo, storing the result back to postgis.<br></div><div>It is very easy to parallelize the process with bash and gnu parallel and I was able to process ~40M polygon buildings against ~270 TIFF files at ~30m resolution containing flood data in about 8 hours on a macbook pro connecting to S3 for raster data. (complete process intersects polygons with each raster pixel and performs global statistics)<br></div><div><br></div><div>An example of command line calling gdallocationinfo from vector table could be, using a result table like this:</div><div><br></div><div><div style="background-color:rgb(43,43,43);color:rgb(169,183,198)"><pre style="font-family:"Consolas",monospace;font-size:10.5pt"><span style="color:rgb(204,120,50)">drop table if exists </span>gdal_res<span style="color:rgb(204,120,50)">;<br></span><span style="color:rgb(204,120,50)">create unlogged table </span>gdal_res (<br>    <span style="color:rgb(152,118,170)">val </span><span style="color:rgb(204,120,50)">float,<br></span><span style="color:rgb(204,120,50)">    </span><span style="color:rgb(152,118,170)">fid </span><span style="color:rgb(204,120,50)">int,<br></span><span style="color:rgb(204,120,50)">    </span><span style="color:rgb(152,118,170)">cgrid_id </span><span style="color:rgb(204,120,50)">int primary key<br></span>)<span style="color:rgb(204,120,50)">;<br></span></pre></div></div><div>Then in bash</div><div><br></div><div><div style="background-color:rgb(43,43,43);color:rgb(169,183,198)"><pre style="font-family:"Consolas",monospace;font-size:10.5pt"><span style="color:rgb(197,118,51)">time </span>psql -X --no-align --quiet --tuples-only --field-separator <span style="color:rgb(106,135,89)">" " </span>--pset footer=off -d mydb -c <span style="color:rgb(106,135,89)">"<br></span><span style="color:rgb(106,135,89)">    select st_x(st_centroid(geom)), st_y(st_centroid(geom)), fid, id<br></span><span style="color:rgb(106,135,89)">    from cgrid" </span>| <span style="color:rgb(197,118,51)">gdallocationinfo </span>-field_sep <span style="color:rgb(106,135,89)">" " </span>-valonly -geoloc \<br>    S3:/bucket/dem09.tif | <span style="color:rgb(197,118,51)">psql </span>-X -d mydb -c <span style="color:rgb(106,135,89)">"COPY gdal_res from STDIN NULL '' delimiter ' '"<br></span></pre></div></div><div><br></div><div>In this way, no intermediate data is written to disk.</div><div>Pixels seen as polygon or bboxes can then be processed against vector data.</div><div>gdallocationinfo in echo mode (-E option) allows you to pass extra information aside from x-y coordinates, like feature identifiers. </div><div>These ids can be used in Postgis to JOIN against vector features very easily/efficiently.<br></div><div>Compared to python approach reading data with Rasterio for instance, performance was waaayyy better.</div><div><br></div><div>Nicolas<br></div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class="msg-8107322994282190670"><div lang="EN-US" style="overflow-wrap: break-word;"><div class="m_-8107322994282190670WordSection1"><p class="MsoNormal"><span style="font-size:11pt">Admittedly my experience with using GDAL and python directly is very low compared to doing everything in the database so I really don’t know the benefits of doing it all in GDAL.<u></u><u></u></span></p><p class="MsoNormal"><span style="font-size:11pt"><u></u> <u></u></span></p><p class="MsoNormal"><span style="font-size:11pt">I still prefer to use the database hammer over having to learn yet another tool.  My sense is also that the amount of SQL/plpgsql code you have to write would be much less than trying to pull your stuff out of the DB and stuff it in python/GDAL to do the same thing.  This is all to say if the performance is good enough for you, I say go for that hammer even if you think you are just dealing with a nail.<u></u><u></u></span></p><p class="MsoNormal"><span style="font-size:11pt"><u></u> <u></u></span></p><p class="MsoNormal"><span style="font-size:11pt"><u></u> <u></u></span></p><p class="MsoNormal"><span style="font-size:11pt"><u></u> <u></u></span></p><div style="border-width:medium medium medium 1.5pt;border-style:none none none solid;border-color:currentcolor currentcolor currentcolor blue;padding:0in 0in 0in 4pt"><div><div style="border-width:1pt medium medium;border-style:solid none none;border-color:rgb(225,225,225) currentcolor currentcolor;padding:3pt 0in 0in"><p class="MsoNormal"><b><span style="font-size:11pt;font-family:"Calibri",sans-serif">From:</span></b><span style="font-size:11pt;font-family:"Calibri",sans-serif"> Thiemo Kellner <<a href="mailto:thiemo@gelassene-pferde.biz" target="_blank">thiemo@gelassene-pferde.biz</a>> <br><b>Sent:</b> Friday, November 8, 2024 12:16 PM<br><b>To:</b> PostGIS Users <<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a>><br><b>Subject:</b> Re: Composing raster tiles?<u></u><u></u></span></p></div></div><p class="MsoNormal"><u></u> <u></u></p><p class="MsoNormal">Thanks for sharing your experience. Is that saying that PostGIS has considerable overhead? Sorry, if that is naiv. I might be suffering from the hammer nail deficiency. My background is very database heavy, so too much might look like a database to me. <u></u><u></u></p><div><div><p>08.11.2024 16:15:32 Vera Green <<a href="mailto:vera.green.ca@gmail.com" target="_blank">vera.green.ca@gmail.com</a>>:<u></u><u></u></p></div><blockquote style="border-width:medium medium medium 2.25pt;border-style:none none none solid;border-color:currentcolor currentcolor currentcolor rgb(204,204,204);padding:0in 0in 0in 8pt;margin-left:0in;margin-right:0in"><p>We use command line GDAL for all our Easter processes, if your data is large I recommend you look into that option.<u></u><u></u></p><p class="MsoNormal"><u></u> <u></u></p><div><div><p class="MsoNormal">On Fri, Nov 8, 2024, 6:12<span style="font-family:"Arial",sans-serif"> </span>AM David Haynes <<a href="mailto:haynesd2@gmail.com" target="_blank">haynesd2@gmail.com</a>> wrote: <u></u><u></u></p></div><blockquote style="border-width:medium medium medium 1pt;border-style:none none none solid;border-color:currentcolor currentcolor currentcolor rgb(204,204,204);padding:0in 0in 0in 6pt;margin-left:4.8pt;margin-right:0in"><div><p class="MsoNormal">I'll try and help you with A & B <u></u><u></u></p><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">a) Is it more efficient to convert the raster to vector data and calculate on the those than to calculate directly on the raster? <u></u><u></u></p></div><div><p class="MsoNormal">>> I don't think it would be faster to convert to vector because that would be dumping the raster into polygons. I think for specifically calculating slope, you are better off staying in raster.  <br><br>b) To my understanding, if I calculate the slope on a raster tile, the slope,… of the borders will have accuracy problems. My I idea, was to "stitch" a tile with its direct neighbours, calculate on the composed  tile, and either save a cropped calculated composed tile to its  original dimension or save the calculated composed tile as is,  probably the latter. <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">>> Overall correct, some small adjustments. For any raster operation, the tiles are operated on independently, so you need to "stitch" them together. I've provided a couple of ways in pseduo code <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Way 1 <u></u><u></u></p></div><div><p class="MsoNormal">1) Use ST_Union and make a big tile,  <u></u><u></u></p></div><div><p class="MsoNormal">2) Do the spatial operation, <u></u><u></u></p></div><div><p class="MsoNormal">3) Break it up using ST_Tile() <br>The downside is you might run out of memory doing this. Also consider ST_MemUnion() <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">SELECT ST_Tile(ST_Slope(ST_Union(r1.ras)), 350,350) <u></u><u></u></p></div><div><p class="MsoNormal">FROM raster_table <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Way 2) <u></u><u></u></p></div><div><p class="MsoNormal">A second option is to basically create an aggregate, which is likely faster.  <u></u><u></u></p></div><div><p class="MsoNormal">1) Union the tiles based on ST_Touch  - make mega_tiles <u></u><u></u></p></div><div><p class="MsoNormal">2) Do the spatial operation on the mega_tiles <u></u><u></u></p></div><div><p class="MsoNormal">3) Clip the mega_tiles by the old tiles bounding box <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">WITH rtest as <br>( <br>SELECT r1.ras, ST_Union (r1.ras) as megatile <br>FROM raster_table r1 <br>LEFT JOIN raster_table r2 <br>ON ST_Touches (r1.geom, r2.geom) <br>GROUP BY r1.ras <br>) <br>SELECT ST_CLIP(ST_SLOPE(megatile), ST_Envelop(r1.rast) ) as ras <br>FROM rtest <u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Maybe someone wants to make an aggregate for the raster functions? <u></u><u></u></p></div></div><p class="MsoNormal"><u></u> <u></u></p><div><div><p class="MsoNormal">On Thu, Nov 7, 2024 at 3:31<span style="font-family:"Arial",sans-serif"> </span>PM <<a href="mailto:thiemo@gelassene-pferde.biz" target="_blank">thiemo@gelassene-pferde.biz</a>> wrote: <u></u><u></u></p></div><blockquote style="border-width:medium medium medium 1pt;border-style:none none none solid;border-color:currentcolor currentcolor currentcolor rgb(204,204,204);padding:0in 0in 0in 6pt;margin-left:4.8pt;margin-right:0in"><p class="MsoNormal" style="margin-bottom:12pt">Hi <br><br>In my project  <br><a href="https://sourceforge.net/p/treintaytres/code/HEAD/tree/trunk/code_files/data_storage/" target="_blank">https://sourceforge.net/p/treintaytres/code/HEAD/tree/trunk/code_files/data_storage/</a> I have the  <br>table <br><br>TABLE_SCHEMA    TABLE_NAME      DATA_TYPE       TYPE_NAME       COLUMN_NAME <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    1111            uuid            id <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    93              timestamptz     entry_pit <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    1111            uuid            source_id <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    12              text            file_name <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    1111            raster          tile <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    93              timestamptz     file_creation_pit <br>treintaytres    topo_files<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>t    12              text            file_hash <br><br>TILE contains topographical height raster data from OpenTopography.  <br>They are from different regions, let's say, some tiles cover  <br>Switzerland, some cover New Zealand. I want to create slope and other  <br>data from the height data and I have some questions I hope you can  <br>answer or point me to answers. <br><br>a) Is it more efficient to convert the raster to vector data and  <br>calculate on the those than to calculate directly on the raster? <br><br>b) To my understanding, if I calculate the slope on a raster tile, the  <br>slope,… of the borders will have accuracy problems. My I idea, was to  <br>"stitch" a tile with its direct neighbours, calculate on the composed  <br>tile, and either save a cropped calculated composed tile to its  <br>original dimension or save the calculated composed tile as is,  <br>probably the latter. <br>Can I compose as follows? <br>with RASTER_NEIGHBORS as (         select R1.TILE   as CURRENT_TILE <br>                                          ,R2.TILE   as NEIGHBOR_TILE <br>                                          ,<a href="http://R1.ID" target="_blank">R1.ID</a>     as CURRENT_ID <br>                                      from TOPO_FILES<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>T R1 <br>                           left outer join TOPO_FILES<span style="font-family:"Segoe UI Symbol",sans-serif">␟</span>T R2 <br>                                        on ST_Touches(R1.TILE <br>                                                     ,R2.TILE) <br>                                        or ST_Intersects(R1.TILE <br>                                                        ,R2.TILE) <br>                                     where TRUE <br>                                       --and <a href="http://R1.ID" target="_blank">R1.ID</a> =  <br>'6b8ca53a-bb5f-4c2b-a9c9-94b6a706e9b0' <br>                                       and TRUE) <br>    ,NION as (select CURRENT_TILE as TILE <br>                    ,CURRENT_ID <br>                from RASTER_NEIGHBORS <br>              union <br>              select NEIGHBOR_TILE as TILE <br>                    ,CURRENT_ID <br>                from RASTER_NEIGHBORS) <br>   select ST_Union(TILE) as COMPOSED_TILE <br>         ,CURRENT_ID <br>     from NION <br>group by CURRENT_ID; <br><br>c) Finally, I want to select all the areas where slope, TRI,… conform  <br>certain criteria and have a minium surface size. Do I do it this  <br>better on vector data and do I need to do this on data composed of all  <br>the contiguous areas? <br><br>I would be grateful for any nudge into the right direction. Maybe URLs  <br>with samples. <br><br>Kind regards <br><br>Thiemo <u></u><u></u></p></blockquote></div></blockquote></div></blockquote></div></div></div></div></div></blockquote></div></div>