Plotting thousands of polygons

General Usage Performance
,
1

Would you recommend some specific configuration for efficiently plotting a large number (say 20k) polygons?

Now we use a custom Plots recipe with a seriestype := :shape and a loop like

...
    for i in indices # loop over polygons
        # get the list of vertices of this polygon
        vlist = hcat(vertices_list(sol.Xk[i])...).' 
        @series (x, y) = vlist[:, 1], vlist[:, 2]
    end

and the GR backend (it is sligtly faster than other backends such as PyPlot but still this is negligible since both take several minutes). This is all good, but I was wondering if this task can be accelerated by using an alternative setup.

2

GMT.jl?

3

Does it work to use a single list and insert NaNs where the breaks between polygons occur? This works for line segments as far as I remember.

4

Alright, i’ve installed GMT.jl and can reproduce the rectangles examples :+1:

Do you have an example for my use case? I tried plot!(vlist, region=R, show=false, fill="blue") written in a for loop (as above). The plots looks ok for N=1000 sets, but the problem is that it pops up a lot of windows, one after the other, so i can’t try higher instances. (FWIW, in Plots you would do something like my_plot = plot(...); plot!(my_plot, ...).) I expect to get one single graph with all the plots.

I also tried a single command plot(vi, region=R, fill="blue", show=false) where vi is the concatenation of all the polygons (vertically), and the plot looks strange. Is there a special argument to separate the different shapes? (like a NaN).

In any case thanks for the pointer to GMT, that I would like to explore further :slight_smile:

5

To answer my own question, yes you can use NaNs to separate the polygons with Plots.jl:

x = [1, 2, 3, NaN, 4, 5, 6]
y = [3, 2, 7, NaN, 1, 9, 2]

using Plots; gr()

plot(x, y, seriestype=:shape)  # plot polygons
2 Likes
6

Indeed, I remember that is a valid input from reading the docs, but didn’t take it seriously as a performance factor… yet it is :smile: For a N = 10000 instance I get a 30x speedup (!).

2 Likes
7

Ah, sorry you found a small bug. With this calls it’s not checking the value of the show keyword and it takes it always as true. So the solution meanwhile is to not use at all (when not needed). Here is an working example that creates a png image:

julia> p1 = [0 0; 1 0; 1 1];

julia> p2 = [3 2; 4 2; 4 3];

julia> p3 = [2 4; 2 5; 3 4];

julia> plot(p1, region=[0 5 0 5], fill=:blue, frame="a")

julia> plot!(p2, fill=:blue)

julia> plot!(p3, fill=:blue, fmt=:png, show=true)

However, this would be terribly inefficient for 1000’s of polygons. So yes you can use NaNs.

julia> pp = [p1; NaN NaN; p2; NaN NaN; p3];

julia> plot(pp, region=[0 5 0 5], fill=:blue, fmt=:png, frame="a", show=true)

There is even a more versatile way where one can specify different line widths, fill colors etc by creating GMTDataset type an filling in the segment headers but will need to consult the code to see if I had already implemented the interface in GMT.jl that allows calling this feature of GMT (no the wrapper but the package itself).

1 Like
8

If you update with Pkg.checkout("GMT") you will be able to this too:

# Create an array of GMTdataset types
DSarr = [GMTdataset() for i = 1:3]

# Same example polygons
p1 = [0 0; 1 0; 1 1]; p2 = [3 2; 4 2; 4 3]; p3 = [2 4; 2 5; 3 4];

# Fill the DSarr array setting different colors to each polygon (Sorry, here wevhave to use pure-and-dure GMT syntax for this)
DSarr[1].data = p1;  DSarr[1].header = "-W1 -Gred";
DSarr[2].data = p2;  DSarr[2].header = "-W1 -Ggreen";
DSarr[3].data = p3;  DSarr[3].header = "-W1 -Gblue";

# Do the plot
plot(DSarr, region=[0 5 0 5], fmt=:png, frame="a", show=true)

You can extend this principle to plot thousands of polygons.

2 Likes