Subscribe to this thread
Home - General / All posts - creating multi-ring, non-overlapping buffers in Manifold 9 SQL
artlembo


3,400 post(s)
#26-Jun-18 17:19

I just posted a new video on how to use SQL in Manifold 9 to create multi-ring, non-overlapping buffers. In this example, larger buffers are trimmed to the edges of the smaller buffers.

The video is 16 minutes long, but I show you the answer at the beginning. If you want, you can watch the rest of the video to see the logic in structuring the SQL. It not only gives you a very useful tool for your GIS toolbox, but helps you think about the logic of how to create a more sophisticated SQL query, one step at a time.

You can see the video here.

StanNWT
196 post(s)
#20-Sep-18 19:38

Hi Art,

I just watched your video, great to see. However, because of the resolution of the video it's hard to see many parts. Regardless of this, I'm a beginner to Manifold SQL9, so I am wondering how I can take my a complex country polygon and create multi-ring non-overlapping buffers to it at the following distances: (500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, 7500, 8000, 8500, 9000, 9500, 10000)? Distances are in meters. I like the little 0.1 interval at the start so that when you have the tolerance of 1, that it wipes away the 0-0.1 buffer completely. I tried to use the buffer transform and modify it at least generate the buffers before removing the overlap, but I get a cannot parse query statement. I know I jumped into the middle with this, trying to use an example of someone's code to apply it to my data, but I figured it was worth a shot. Below is the code I tried to write in the command window:

The source for the data is still in an Esri File GeoDatabase.

-- $manifold$

-- Auto-generated

-- Transform - Buffer - updateField

PRAGMA ( ' progress.percentnext ' = ' 100 ' ) ;

UPDATE (

SELECT [OBJECTID] ,

[shape ] ,

GeomBuffer( [shape ], result, O) AS gl, result AS bufdist

INTO bufdist

FROM [Esri_2018_2019_General_File_GeoDatabase] :[AC_lM_Bndry_NWT_onlyDiss_NoEskLak_NoHrbrEntr_Is_union2_Dissolved2Drawing] , (VALUES (0.1), (500), (1000), (1500), (2000), (2500), (3000), (3500), (4000), (4500), (5000), (5500), (6000), (6500), (7000), (7500), (8000), (8500), (9000), (9500), (10000));

THREADS systemcpucount ()

) SET [shape] = [n_shape] ;

It's quite likely I'm light years away from getting this right, but since you made a nice video for this, I was wondering if you can point out where I went wrong and what I might need to do to get it right?

The reason for my attempt at this is to assess the time it takes to do these multi-ring buffers in Manifold 9 vs. ArcGIS Desktop 10.5.1 vs. ArcGIS Pro 2.2.2. I already suspect who the winner is. I just like hard data.

I think a great modification to the Buffer transform would be an option for adding rows in the transform for the buffer distances, with a radio or checkbox for internal or external or both. This would take care of all the stuff you've written in SQL, would it not?

tjhb
10,094 post(s)
#20-Sep-18 23:56

--SQL9

VALUE @drawing TABLE = [Drawing];

VALUE @table TABLE = [Drawing Table];

VALUE @w FLOAT64 = 500; -- ring width

VALUE @max FLOAT64 = 10000; -- max radius

SELECT

    t.[mfd_id] AS [source_id],

    s.[Value] AS [radius],

    CASE

        WHEN s.[Value] > 0 THEN

            GeomMergeAreasPair(

                GeomBuffer(t.[Geom], s.[Value], 0),

                GeomBuffer(t.[Geom], s.[Value] + @w, 0)

                )

        ELSE -- first ring

            CASE

                WHEN GeomIsArea([Geom]THEN

                    GeomClip(

                        GeomBuffer(t.[Geom], @w, 0), 

                        t.[Geom],

                        FALSE, -- subtract

                        0)

                ELSE -- line or point

                    GeomBuffer(t.[Geom], @w, 0)

            END

    END AS [Ring]

INTO [Rings Table]

FROM

    @table AS t

    CROSS JOIN

    (TABLE CALL ValueSequence(0, @max, @w)) AS s

THREADS SystemCpuCount() / 2 + 1

;

ALTER TABLE [Rings Table] (

    ADD [mfd_id] INT64-- filled automatically

    ADD INDEX [mfd_id_x] BTREE ([mfd_id]),

    ADD INDEX [Ring_x] RTREE ([Ring]),

    ADD PROPERTY 'FieldCoordSystem.Ring' ComponentCoordSystem(@drawing)

    );

CREATE DRAWING [Rings] (

    PROPERTY 'Table' '[Rings Table]',

    PROPERTY 'FieldGeom' 'Ring'

    );

You will need to adjust it:

(1) Change the first two constants to match your source drawing

VALUE @drawing TABLE = [Esri_2018_2019_General_File_GeoDatabase] :[AC_lM_Bndry_NWT_onlyDiss_NoEskLak_NoHrbrEntr_Is_union2_Dissolved2Drawing];

and its table.

(2) Replace t.[mfd_id] with t.[OBJECTID.

(2) Replace t.[Geom] with t.[shape ].

That should be it.

Attachments:
concentric buffers.sql

tjhb
10,094 post(s)
#21-Sep-18 00:21

There is a missing colon in the drawing path above. Should be

--SQL9

VALUE @drawing TABLE = [Esri_2018_2019_General_File_GeoDatabase]::[AC_lM_Bndry_NWT_onlyDiss_NoEskLak_NoHrbrEntr_Is_union2_Dissolved2Drawing];

StanNWT
196 post(s)
#21-Sep-18 16:45

Thanks very much!

I decided to copy and past the data from the Esri File GeoDatabase into the (.map) file, so that I wasn't running the operation accessing an Esri dataset but a data set directly inside Manifold 9. I'm using 9.0.168.1 right now to process the multi-ring buffers. I renamed all the drawing, table and Geom and mfd_id bits in the code. It took a couple tries, as I was missing things here and there, but once I tracked down where I missed some change that needed to be made, it's now running. I copied the data from my Drobo 5D which gets 250MB/s write and 225 MB/s read, to my Samsung 960 Pro 2TB <C:> drive which has 3GB/s read adn 2.7GB/s write. I hope it will run faster. The data is in a large Manifold (.map) file, that's 152GB in size, as it has a large DEM in it. I supposed I could have done it with a (.map) file that only has that single polygon of the coastline that needs the multi-ring buffer created in that project, however, I tend to work in projects that have lots of data in them. I'm not sure if that will negatively impact the processing time?

I'll let you know how it goes, if it's successful and the processing time.

tjhb
10,094 post(s)
#21-Sep-18 18:47

Since you only have a single, very large area (which I admit I hadn't appreciated), it would be faster to comment out the THREADS directive (and just use one thread).

Alternatively, you could split the country area to shapes (assuming there are offshore islands) and leave THREADS in. But most of the work would still be for one big area, and that will not parallelize efficiently.

A better approach for data like this would be to calculate all buffers first, in a multi-threaded subquery, then assemble in an outer layer.

That would also avoid a redundancy, where each buffer (except the last) is currently calculated twice, which is poor design on my part (and matters most for one big area).

Once the buffers are calculated there is little work left to do, since only the first ring needs to be clipped. The other rings are simply combined in pairs--Manifold assumes that the smaller ring is a hole (which is true).

Do let me know how you get on with this code, then we can improve it. Attach the code after your edits if you like, so I can match them.

There is no need to worry about large projects affecting performance (or productivity) in Manifold 9. On the other hand, splitting data up is also more productive in Manifold 9 than in 8, since projects can be infinitely nested as modules (and recombined and reused). Shared modules can also be separately updated. It's pretty cool.

tjhb
10,094 post(s)
#21-Sep-18 19:32

My first para above is wrong. I wasn't quite awake yet.

Keep using multiple threads, as written. The 20 buffer pairs will be divided among the thread pool. That's fine.

But the redundancy is an issue that needs fixing.

StanNWT
196 post(s)
#21-Sep-18 21:22

Here's my code as executed:

--SQL9

VALUE @drawing TABLE = [AC_1M_Bndry_NWT_OnlyDiss_NoEskLak_NoHrbrEntr_Is_Union2_Dissolved2 Drawing];

VALUE @table TABLE = [AC_1M_Bndry_NWT_OnlyDiss_NoEskLak_NoHrbrEntr_Is_Union2_Dissolved2];

VALUE @w FLOAT64 = 500; -- ring width

VALUE @max FLOAT64 = 10000; -- max radius

SELECT

t.[OBJECTID] AS [source_id],

s.[Value] AS [radius],

CASE

WHEN s.[Value] > 0 THEN

GeomMergeAreasPair(

GeomBuffer(t.[Shape], s.[Value], 0),

GeomBuffer(t.[Shape], s.[Value] + @w, 0)

)

ELSE -- first ring

CASE

WHEN GeomIsArea([Shape]) THEN

GeomClip(

GeomBuffer(t.[Shape], @w, 0),

t.[Shape],

FALSE, -- subtract

0)

ELSE -- line or point

GeomBuffer(t.[Shape], @w, 0)

END

END AS [Ring]

INTO [Rings Table]

FROM

@table AS t

CROSS JOIN

(TABLE CALL ValueSequence(0, @max, @w)) AS s

THREADS SystemCpuCount() / 2 + 1

;

ALTER TABLE [Rings Table] (

ADD [mfd_id] INT64, -- filled automatically

ADD INDEX [mfd_id_x] BTREE ([mfd_id]),

ADD INDEX [Ring_x] RTREE ([Ring]),

ADD PROPERTY 'FieldCoordSystem.Ring' ComponentCoordSystem(@drawing)

);

CREATE DRAWING [Rings] (

PROPERTY 'Table' '[Rings Table]',

PROPERTY 'FieldGeom' 'Ring'

);

Here's the log:

2018-09-21 11:27:47 -- Manifold System 9.0.168.2 Beta

2018-09-21 11:27:47 -- Starting up

2018-09-21 11:27:47 Log file: C:\Users\stan_johnston\AppData\Local\Manifold\v9.0\20180921.log

2018-09-21 11:27:50 -- Startup complete (2.770 sec)

2018-09-21 11:27:50 -- Create: (New Project) (0.014 sec)

2018-09-21 11:28:05 -- Close: (New Project) (0.006 sec)

2018-09-21 11:28:08 -- Open: K:\USGS_NED_Fused_1_and_2_arc_second_DEM_with_Contours_MultiringBuffers_Coastline_Sept20_2018.map (3.036 sec)

2018-09-21 13:40:35 -- Query: (ad-hoc) (7919.729 sec)

2018-09-21 14:03:59 Render: [Rings] (0.587 sec)

2018-09-21 14:04:04 Render: [Rings] (0.492 sec)

2018-09-21 14:04:06 Render: [Rings] (0.331 sec)

2018-09-21 14:05:00 Render: [Rings] (0.343 sec)

2018-09-21 14:05:49 Render: [Rings] (0.333 sec)

2018-09-21 14:05:52 Render: [Rings] (0.343 sec)

Sorry the formatting didn't seem to come across properly in the code that I copied and pasted from Manifold's command window. But it's your code with the parameters changed.

You can see the results are exactly what I wanted.

Attachments:
Mainfodl_Multiring_Buffer_Result1.jpg
Mainfodl_Multiring_Buffer_Result2.jpg
Mainfodl_Multiring_Buffer_Result3.jpg
Mainfodl_Multiring_Buffer_Result4.jpg
Mainfodl_Multiring_Buffer_Result5.jpg

tjhb
10,094 post(s)
#22-Sep-18 00:04

If you have time, could you test this version as well? It should be more efficient, unless I've cocked up.

--SQL9

VALUE @drawing TABLE = [AC_1M_Bndry_NWT_OnlyDiss_NoEskLak_NoHrbrEntr_Is_Union2_Dissolved2 Drawing];

VALUE @table TABLE =  [AC_1M_Bndry_NWT_OnlyDiss_NoEskLak_NoHrbrEntr_Is_Union2_Dissolved2];

VALUE @w FLOAT64 = 500; -- ring width

VALUE @max FLOAT64 = 10000; -- max radius

VALUE @buffers TABLE =

    (

    SELECT

        t.[OBJECTID] AS [source_id],

        t.[Shape],

        s.[Value] AS [radius],

        CASE s.[Value]

            WHEN 0 THEN NULL

            ELSE GeomBuffer(t.[Shape], s.[Value], 0)

        END AS [Buffer]

    FROM

        @table AS t

        CROSS JOIN

        (TABLE CALL ValueSequence(0, @max, @w)) AS s

    THREADS SystemCpuCount() / 2 + 1

    BATCH 1 -- assign one source record at a time per thread

    );

SELECT

    b1.[source_id],

    b2.[radius],

    CASE

        WHEN b1.[radius] > 0 THEN

            GeomMergeAreasPair(b1.[Buffer], b2.[Buffer])

        ELSE -- first ring

            CASE

                WHEN GeomIsArea(b1.[Shape]THEN

                    GeomClip(b2.[Buffer], b1.[Shape], FALSE, 0)

                    -- FALSE -> subtract

                ELSE -- line or point

                    b2.[Buffer]

            END

    END AS [Ring]

INTO [Rings Table]

FROM

    @buffers AS b1 INNER JOIN @buffers AS b2

    ON b1.[radius] = b2.[radius] - @w

    -- last row from b1 is excluded (unmatched)

THREADS SystemCpuCount() / 2 + 1

BATCH 2

;

ALTER TABLE [Rings Table] (

    ADD [mfd_id] INT64-- filled automatically

    ADD INDEX [mfd_id_x] BTREE ([mfd_id]),

    ADD INDEX [Ring_x] RTREE ([Ring]),

    ADD PROPERTY 'FieldCoordSystem.Ring' ComponentCoordSystem(@drawing)

    );

CREATE DRAWING [Rings] (

    PROPERTY 'Table' '[Rings Table]',

    PROPERTY 'FieldGeom' 'Ring'

    );

Attachments:
Draw concentric buffers b forum.sql

adamw


10,447 post(s)
#24-Sep-18 16:09

You can perhaps speed things up by not doing GeomClip - and doing GeomMergeAreas in place of it, no? (The clipping is between the smallest buffer and the area it has been built on, or am I mistaken?)

But in general, our buffer function needs a rework. It is robust and safe, but it is not fast, we have long wanted to redo it - and will do so (hopefully in the near future, because we need a fast buffer in tons of places).

tjhb
10,094 post(s)
#25-Sep-18 03:12

Thanks for that suggestion Adam. I had thought of it, but for some reason I wasn't as confident doing that with the source area as with the subsequent buffer pairs. As you say, it works fine. I've incorporated it below.

I'll try to provide some decent timings of buffers between Manifold 8 and Manifold 9. I don't know in advance which will be faster. For complex geometry these multi-ring buffers certainly don't look like a trivial amount of work. I have probably got too used to Manifold 9 performance seeming miraculous.

tjhb
10,094 post(s)
#22-Sep-18 00:22

By the way, the previous version was creating one ring too many.

They went from radius 0 to radius 10000, where the radius described the inner edge.

In the revised version, rings run from radius 500 to radius 10000, where the radius describes the outer edge of the ring.

StanNWT
196 post(s)
#22-Sep-18 05:10

I actually did want buffer rings with a width of 500 m from the coastline out to 10,000 m. So the 0-500 m buffer ring is required.

I have a dual hexacore with hyperthreading so I'm wondething if the thread count setting is restricting speed or if given the nature of what is being processed it's irrelevant?

This multi-ring buffer drawing is used to create a halo effect where each buffer ring will have a different transparency integer value so that I can create a glow around the coastline. Is great when you have a bathymetric DEM that may or may not have some issues when conspired to a much higher resolution DEM for the land. With the halo you can create a nice cartographic effect at the land water transition. I can show you an example that I have used in ArcGIS Desktop for years and more recently in ArcGIS Pro.

tjhb
10,094 post(s)
#22-Sep-18 06:22

Both queries produce the 0-500 ring, and all intermediate rings out to 9500-10000.

My first query also produced a ring 10000-10500, which it should not have.

The main difference between the two queries is that the second does not create duplicate buffers. That ought to be a significant saving. If you have time to run it and report, I'd be grateful.

By all means experiment with more cores (or fewer). Given hyperthreading, and with a task that can use all available cores, my rule of thumb is to assign threads for half the number of logical cores (= all physical cores) plus one. On my four physical/eight virtual core machines, adding threads beyond five normally slows things down progressively. (This shows, by the way, how well the Manifold scheduler is written. Intel can find hardly any slack to redeploy.) So with 6 physical cores, I would start with 7--but only for a task that can use all available cores. There is a lot of detail to that. (And you'll see that for the second query I have restricted the number of records fed to each thread to 1 and 2 respectively. That's unusual, but I think it works well in this case.)

That is task-dependent, data-dependent, and of course hardware-dependent. Definitely worth experimenting.

tjhb
10,094 post(s)
#22-Sep-18 19:44

My second query has a major bug, really stupid.

Don't test it--I will.

StanNWT
196 post(s)
#24-Sep-18 16:56

Thanks for your help on this, it's greatly appreciated.

My workstation is a dual hexacore so I have 12 physical cores, 24 logical cores. I could easily assign 80% of my cores to work on this. I've got 128GB of RAM, my data is on a M.2 SSD that reads at 1400 MB/s, writes at 300 MB/s, my <C:> drive is an M.2 SSD that reads at 3000 MB/s writes at 2700 MB/s. I have put the manifold project on the slower drive because it's not backed, up and I don't want to be backing up a 152 GB file if I don't have to. I have a copy of the project on a slower RAID.

What did you think of my idea for check boxes and extra rows to allow the variables for the extra rings in the ring buffer in the transform itself? Or perhaps having a min, max and equal interval values option in the transform? I'm just thinking for those of us that don't yet know how to write the SQL or scripts necessary.

My goal for this data set is to create a glowing halo effect around polygons, either external or internal multi-ring buffers. I use these for cartographic effect, where I want to create a gradient fill effect similar to what you can do in vector graphic design applications, in GIS applications that don't support it.

Do you want this polygon to test things for yourself?

Here's an example map that shows how I use this.

I've attached a screen grab from MS Excel, which shows an internal multi-ring buffer example, with the resulting transparency curve that I apply to each interval so that I get a fading effect to each successive buffer interval to get the kind of gradient fill transparency effect I'm looking for. I do this all the time in ArcGIS, but would like to be able to do it in Manifold, and I also want and know it can be faster in Manifold. ArcGIS Desktop often crashes when the geometry is too large and complex, and doing it in ArcGIS Pro can be problematic trying to create feature classes that I can use in Desktop. Basically, in ArcGIS, I create the multi-ring buffers, then create a transparency field according to the values I've specified in a MS Excel spreadsheet an example of which is attached. Then the transparency setting for the layer properties allows me to have different transparency values for each buffer ring. I can also apply a secondary transparency for the whole layer if I want to fine tune the overall transparency.

Attachments:
Example_of_Glowing_Multirng_Buffer_Effect_Sept24_2018.jpg
Transparency_Curve_Values_Internal_Buffer_Spet24_2018.jpg

tjhb
10,094 post(s)
#25-Sep-18 00:34

That's really exemplary cartography, if you don't mind me saying so. (Looking at the PDF.) Beautiful terrain shading and hypso, brilliant choice of overlay colours (helped by matching fill styles in places)--overlays bright and clear, without obscuring the underlying terrain at all--and the custom projection makes great use of page space and looks natural (with equally clear inlay maps)... really the works. A beautiful, informative map.

I have been travelling for a day.

The bug in my second query was really dopey--I forgot to join back on the source ID field, as well as on radius. So assuming you have more than one source area, every one would get multiple copies of every buffer ring.

Easily fixed, just one line. But then in my own testing (using the New Zealand coastline, which is hardly Canada, but our countries are anyway good friends), I was really surprised how slow it proved in practice.

Adam has helped substantially on that point above. It seems this will change when the buffering code is further optimized.

I don't want to labour that point but I think it is worth posting some direct timings comparing 8 with 9 (also for the somewhat related thread on clipping).

So I will come back soon with data, code for 8, and timings.

Thanks for your offer to provide your own data Stan. I don't need it (New Zealand is probably an adequate proxy for testing), but if you do provide it then I will test with that too.

In the meantime I have attached the second query above with the bug corrected, and adding in Adam's suggestion to use GeomMergeAreasPair, instead of GeomClip, for the first ring as for subsequent rings. That also allows one CASE expression to be removed.

The SELECT INTO query is now like this. (There are no other material changes.)

--SQL9

--...

SELECT

    b1.[source_id],

    b2.[radius],

    CASE

        WHEN GeomIsArea(b1.[Shape]THEN

            GeomMergeAreasPair(b1.[Buffer], b2.[Buffer])

            -- inner and outer areas combined

        ELSE

            [b2].[Buffer]

            -- for first buffer only

            -- iff source object is line or point

    END AS [Ring]

INTO [Rings Table]

FROM

    @buffers AS b1 INNER JOIN @buffers AS b2

    ON b1.[source_id] = b2.[source_id]

    AND b1.[radius] = b2.[radius] - @w

    -- last row from b1 is excluded (unmatched)

THREADS SystemCpuCount() / 2 + 1

BATCH 2

;

--...

Attachments:
Draw concentric buffers d forum.sql

StanNWT
196 post(s)
#27-Sep-18 22:33

Thanks for the feedback! I've been doing GIS and obviously cartography since 1992/3. I still feel it's critically important to be able to effectively and unambiguously visual information in a map. Far too few people have a sound background in geography and graphic design. It seems like people favour programming almost exclusively, much to my chagrin.

The (.mbx) file is 152GB with the buffering I ran. Do you want a (.mbx) file that is just the outline only? In either case it will be larger than I can upload on here directly. Do you have a way I can share a dropbox link with you?

tjhb
10,094 post(s)
#27-Sep-18 22:42

Yes, the outline only would be great. The email link in my profile is written backwards (a pathetic attempt to stop bots) if you want to share a link that way. Or I can share a Dropbox folder with you.

You might also like to have a look at the website of Roger Smith at www.geographx.co.nz, who taught me a tiny slice of what he knows about cartography. He would like your work and I think you would like his. (Oh--and I strongly suspect he would agree with your comments and chagrin.)

Dimitri


7,413 post(s)
#28-Sep-18 08:51

Thanks for sharing Roger Smith's site. That's really beautiful work!

StanNWT
196 post(s)
#11-Oct-18 19:36

I can send you some links for the other larger DEM data sets as well so you have some bigger data sets to play with?

tjhb
10,094 post(s)
#12-Oct-18 02:51

Well I'd certainly be interested, if you can be bothered.

Thanks very much for the landmass/buffers file. If I get time I'll see if I can speed it up somehow.

StanNWT
196 post(s)
#19-Oct-18 21:15

Just sent you two more files to download.

StanNWT
196 post(s)
#25-Oct-18 16:55

Tjhb,

Did you get the other two projects?

tjhb
10,094 post(s)
#25-Oct-18 19:48

Thanks for checking. Not yet, since I've been away for a few days with no internet access. I'll get them today.

StanNWT
196 post(s)
#30-Oct-18 20:10

Did you download them yet?

tjhb
10,094 post(s)
#30-Oct-18 22:58

Sorry for the delay. One file was successful. The 49 GB file failed at 35.9 GB. Retrying hasn't worked so far.

tjhb
10,094 post(s)
#30-Oct-18 23:14

...But I have saved that file to my own Dropbox account, successfully, so I can retry from there, and you can delete both files if you like. Thanks very much for them.

tjhb
10,094 post(s)
#31-Oct-18 00:17

The download failed from my account too. Never mind.

StanNWT
196 post(s)
#01-Nov-18 17:36

Want me to try and upload a (map) file instead of the (mbx)? They're considerably larger. Did you try and download the file again? Is it the download failing or timing out or is the file corrupt?

StanNWT
196 post(s)
#06-Nov-18 17:27

Do you want me to try and upload the (map) file for the one you couldn't get to download correctly as an (mbx), or do you want me to try and replace the potentially broken (mbx), with a new one?

Manifold User Community Use Agreement Copyright (C) 2007-2021 Manifold Software Limited. All rights reserved.