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:
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!
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
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.
There is also https://github.com/JuliaGeo/Shapefile.jl/pull/57
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…
@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.
@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.
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.
