Best Way of Handling Shapefiles in Makie

Hi all,

What is the best way handle shapefiles with Makie?

For example, I have a shapefile that defines the states for each state in the United States. DataFrames can handle displaying the shapefile as follows:

julia>states
56×10 DataFrame
 Row │ geometry               STATEFP  STATENS   AFFGEOID     GEOID   STUSPS  NAME                               LSAD    ALAND          AWATER
     │ Polygon…?              String   String    String       String  String  String                             String  Int64          Int64
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ Polygon(1516 Points)   31       01779792  0400000US31  31      NE      Nebraska                           00       198956658395    1371829134
   2 │ Polygon(1847 Points)   53       01779804  0400000US53  53      WA      Washington                         00       172112588220   12559278850
   3 │ Polygon(434 Points)    35       00897535  0400000US35  35      NM      New Mexico                         00       314196306401     728776523
   4 │ Polygon(1260 Points)   46       01785534  0400000US46  46      SD      South Dakota                       00       196346981786    3382720225
   5 │ Polygon(3551 Points)   48       01779801  0400000US48  48      TX      Texas                              00       676653171537   19006305260   

Currently, if I try to plot polygons representing each shape using the poly function using the geometry column, I get this error:

ERROR: `Makie.convert_arguments` for the plot type Scatter{ArgType} where ArgType and its conversion trait PointBased() was unsuccessful.

The signature that could not be converted was:
::Shapefile.Point

Makie needs to convert all plot input arguments to types that can be consumed by the backends (typically Arrays with Float32 elements).
You can define a method for `Makie.convert_arguments` (a type recipe) for these types or their supertypes to make this set of arguments convertible (See http://makie.juliap
lots.org/stable/documentation/recipes/index.html).

Alternatively, you can define `Makie.convert_single_argument` for single arguments which have types that are unknown to Makie but which can be converted to known types and
fed back to the conversion pipeline.      

The way I can get around this is by doing this conversion:

states = joinpath(data_path, shape_file) |> Shapefile.Table |> DataFrame
state = states[1, :].geometry.points
points = [(state[i].x, state[i].y) for i in 1:length(state)]
points = Point2f0[points...]

f = Figure()
Axis(f[1, 1])
poly!(points)

Is this the easiest way to handle Shapefiles in Makie? Are there better ways of doing this? Thanks!

~ tcp :deciduous_tree:

I think this way you would lose possible holes in your polygons. A correct translation would be to a Polygon type from GeometryBasics, not a vector of points (which can’t have holes). But I’m not super well versed in the GeometryBasics infrastructure either.

Here I tried something out, the conversion method is inefficient because the GeoInterface seems to create lots of little arrays. But hopefully it’s correct:

using Shapefile
using CairoMakie
using DataFrames
using GeometryBasics: Polygon

# I downloaded this manually and extracted the shape file
url = "https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_state_5m.zip"

shapefile = "Downloads/cb_2018_us_state_5m/cb_2018_us_state_5m.shp"
table = DataFrame(Shapefile.Table(shapefile))
filter!(:NAME => (
    x -> x ∉ [
        "Alaska",
        "American Samoa",
        "Guam",
        "Commonwealth of the Northern Mariana Islands",
        "United States Virgin Islands",
        "Hawaii",
        "Puerto Rico",
        ]),
    table)

function Makie.convert_arguments(::Type{<:Poly}, p::Shapefile.Polygon)
    # this is inefficient because it creates an array for each point
    polys = Shapefile.GeoInterface.coordinates(p)
    ps = map(polys) do pol
        Polygon(
            Point2f0.(pol[1]), # interior
            map(x -> Point2f.(x), pol[2:end]))
    end
    (ps,)
end

let
    f = Figure()
    ax = Axis(f[1, 1])
    foreach(table.geometry) do geo
        poly!(ax, geo)
    end
    f
end

3 Likes

Alternatively, you can rely on the infrastructure surrounding the GeoStats.jl project. You can load the shapefile using GeoTables.jl and it will produce geometries from the Meshes.jl project. Then you can use the MeshViz.jl recipes for Makie to visualize any kind of geospatial domain. The good thing about this solution is that it has a lot of options to color the boundaries of states, the polygons, vertices, etc. Special attention must be given to the decimation option, which is still needed currently to simplify polygons with thousands of vertices that are available on the internet, such as the polygons from the GADM project.

2 Likes

There is also [sd/geometrybasics] Simplify Polygon code and add tests by greimel · Pull Request #57 · JuliaGeo/Shapefile.jl · GitHub
It has been open forever, since it’s not as sure anymore if we should switch Shapefile to GeometryBasics…
Maybe we should just put that in a registered fork or so…

1 Like

@jules - thanks for that very comprehensive answer.
It didn’t even dawn on me that such gapping could occur in my polygons.
Marking your post as a solution for now.

@juliohm - this looks like a very robust pipeline that I wasn’t aware of and very straightforward.
Will experiment with this.

Also, @visr commented that GeoMakie.jl might be worth looking into here too as well.

@juliohm – when I apply the following:

using GeoTables

NZL = GeoTables.gadm("NZL", children = true)
viz(NZL.geometry, decimation = 0.02)

obtained from here

I am returning a stream error as:

image

Do I need to reference GeoStats.jl, if so, what components do
I need to abstract to perform the instructions above?

Thanks,

@YummyPampers2 are you sure the first line of code downloaded the geometries from the web successfully? You have to confirm with y in the prompt if this is the first time you are trying to visualize NZL. GeoTables.jl will download the data for the country of interest if it is not already on disk.

I did run the two lines you shared and it worked fine after confirmation. Please pay attention to the messages that are shown on the terminal at each step. Are you running this on the REPL or Pluto? Maybe Pluto is hidding these messages from you?

Another option is to set the environment variable ENV["DATADEPS_ALWAYS_ACCEPT"] = true at the top of your code so that the confirmation is automatic.

@YummyPampers2 I’ve updated the underlying GADM.jl package to avoid this prompt from DataDeps.jl. I know many users will try this package on Pluto and asking beginners to set env variables is not ideal.

The updated version should be out in a couple of minutes: New version: GADM v0.2.3 by JuliaRegistrator · Pull Request #50735 · JuliaRegistries/General · GitHub

1 Like

@juliohm Thank you for your effort.

When I am using:

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
NZL = GeoTables.gadm("NZL", children = true)
viz(NZL.geometry, decimation = 0.02)

I am getting an error regarding the viz
method as:

image

You need to import the packages:

using GeoTables
using MeshViz

# choose a Makie backend
import GLMakie as Mke

@juliohm Good Day Dr. H:

I was able to produce the following:

I have since modified the code to
include:

viz(NZL.geometry, decimation = 0.10, color = 1:nelements(NZL.geometry), 
    showboundary=true, colormap = :viridis)

I am noticing it is taking up a lot of RAM, everything is
is lagging when I am running the kernels.

My intention is to produce something
similar to a question I presented a
colleague days ago THREAD.

Also are there any wrappers that allow for
some interactive display rather than the
static one we rendered?

The goal is two-fold:

  1. Highlight regions of interest.

  2. Add annotations at centroids
    within those regions to explain certain
    characteristics of those regions.

    • An annotation in a card would be
      better for presentation purposes if
      possible.

Please start a new thread on Discourse. Also these questions have answers that you can find by yourself by searching and experimenting a bit more with the code.

1 Like