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.