Help for improvement for drawing shapefile using Gadfly

Specific Domains Visualization
,
1

I use Shapefile.jl parse *.shp file into DataFrame, with each row having some region information and a Polygon struct specifies the shape of that region. Everything works smoothly. Big thanks :kissing_smiling_eyes:
After that, I needs to plot those Polygons(since every region has its own polygon, here take Australia as example, it has 11 regions, so 11 Polygons to plot) using Gadfly

# this is the data structure of Polygon parsing by Shapefile
struct Polygon
  # the box boundary contains that polygon 
  MBR::Tuple{Float64}
  # each part specify the start of next shape, my understanding 
  parts::Array{Int, 1}
  # Point with x, y field for latitude and longitude
  points::Array{Point, 1}

/* toy example 
polygon.parts == [0, 3, 7] 
polygon.points == [p1, p2, ... , p10]

parts is 0-indexed instead of 1-index
it means 
  * point 0 to 2 defines a shape
  * 3 to 6 is another
  * lastly 7 to 10 
convert to 1-index 1:3, 4:7, 8:10
*/

Remember, for Australia, we have 11 such Polygons, so total #shapes is 6708. My understanding is that I need to plot 6708 polygons. And this is the code

# plot single Polygon and return Array{layer, 1} with each layer represents a shape 
# in the toy example case, that's 3 layers
function plotPolygon(poly)
  starts = poly.parts .+ 1
  ends = length(poly.parts) == 1 ? [length(poly.points)] : vcat(poly.parts[2:end], length(poly.points))
  xs = [p.x for p in poly.points]
  ys = [p.y for p in poly.points]; layers = []
  for part in 1:length(poly.parts)
    s = starts[part]
    e = ends[part]
    l = layer(x = xs[s:e], y = ys[s:e], Geom.polygon(preserve_order = true))
    push!(layers, l)
    end 
  layers
end

# in toy example, one needs  to plot 11 Polygons in one plot
function plotPolygons(polygons::Array{Polygon, 1})  
  p = plot()
  for poly in polygons 
    for layer = plotPolygon(poly)
      push!(p, layer)
    end 
  end 
  p
end 

plotPolygons(df[!, :polygon])

When that runs in repl, it automatically displays to browser, with long delay(and noise from my laptop). But when I try to export that to svg, thinking maybe the delay comes from browser(it turns out not), it StackOverFlow.

plotPolygons(df[!, :polygon]) |> SVG("Australia.svg", 20cm, 20cm)

Is 6708 polygons large? Or is it 1426114 points?
That must be a better way, right?

Iā€™m aware the existence of VegaLight.jl and OpenStreetMapXPlot.jl. But I just want to keep things simple and do it in Gadfly. So Any idea for improvement?(yep, a new laptop could do that)

2

Iā€™m not aware how to do this in Gadfly. But with Plots you could do

plot(myshapefile.shapes...)

GeoMakie also has good support for Shapefiles I believe. I read that you want to do it in Gadfly and donā€™t want to advice you to to anything else, it was just your request to ā€œkeep it simpleā€ that prompted the answer.

2 Likes
3

Welcome to the Julia community!

In Gadfly, drawing many layers is expensive. In my example there is only 1 layer, containing 10 polygons:
Give your polygons an id variable.

using DataFrames, Gadfly

# Construct some polygons for this example 
npolys = 10 # 6708 also works fine
z = [range(-Ļ€/4, stop=6Ļ€/4, length=6)+ 0.2*randn(6) for i in 1:npolys]
xc, yc = 10*rand(npolys), 10*rand(npolys)

D = reduce(vcat, [DataFrame(x=cos.(a).+x, y=sin.(a).+y, id=i)  
        for (i, a, x, y) in zip(1:npolys, z, xc, yc)])
first(D, 4)

p = plot(D, x=:x, y=:y, group=:id,
    Geom.polygon(preserve_order=true),
)

polygons

4

Many thanks! For simplicity I mean to avoid compiling another plotting package just to achieve the same functionality and no offense. :wink:

5

Got it! :kissing_smiling_eyes:
Iā€™ve thought layer is nothing but another element, since with same order 0, they will be rendered just iteratively.