However as a newbie I might be unable to understand how to do this. I have a GeoTable of points and I want to plot the points coloured by a column in the GeoTable.
I can do this using: Makie.viz(data.geometry, colour=data.values). However this doesn’t have equal x and y coordinates which is needed in the UTM CRS I’m using.
I can get the desired result as follows:
fig = Mke.Figure()
ax1 = Mke.Axis(fig[1, 1], aspect = Mke.DataAspect())
x = [coords(pt).x for pt in data.geometry]
y = [coords(pt).y for pt in data.geometry]
Mke.scatter!(ax1, x, y, color = data.values)
fig
However this seems very clunky for something that in R can be done (using the great sf and tmap packages):
tm_shape(data) + tm_dots(col = "values")
I imagine the nicest way to achive this would be some way to just do Mke.geoscatter!(ax1, data.geometry, color = data.values) without the need to extract the x and y values manually.
Assume a dataframe df with :geometry and :is_trauma_center, typeof IGeometry and Bool. It was read from shapefile with GeoDataTables.read()
using CairoMakie
using GeoMakie
using Geometry Basics
# Function to safely calculate centroids for polygons
function safe_centroid(geom)
try
# Try ArchGDAL centroid first
cent = ArchGDAL.centroid(geom)
return (ArchGDAL.getx(cent, 0), ArchGDAL.gety(cent, 0))
catch e
try
# Fallback: use GeoInterface to get coordinates and calculate centroid manually
coords = GeoInterface.coordinates(geom)
if !isempty(coords) && !isempty(coords[1]) && !isempty(coords[1][1])
# For polygons, use the first ring (exterior)
ring_coords = coords[1][1]
if length(ring_coords) > 0
# Calculate centroid as mean of coordinates
x_sum = sum(coord[1] for coord in ring_coords)
y_sum = sum(coord[2] for coord in ring_coords)
n = length(ring_coords)
return (x_sum/n, y_sum/n)
end
end
catch e2
println("Warning: Could not calculate centroid for geometry: $e2")
end
return (0.0, 0.0) # Fallback coordinates
end
end
# Create the figure
f = Figure(size = (1400, 1000))
# Main map for CONUS
ga = GeoAxis(f[1, 1:3]; dest=conus_crs)
# Plot all counties in a neutral color (light gray)
poly!(ga, conus.geometry, color = :lightgray, strokecolor = :white, strokewidth = 0.5)
# Add dots for counties with Level 1 trauma centers
trauma_counties = subset(conus, :is_trauma_center => ByRow(x -> x === true))
if nrow(trauma_counties) > 0
# Calculate centroids for trauma center counties
trauma_centroids = [safe_centroid(geom) for geom in trauma_counties.geometry]
trauma_x = [coord[1] for coord in trauma_centroids]
trauma_y = [coord[2] for coord in trauma_centroids]
# Filter out invalid coordinates
valid_indices = [i for i in 1:length(trauma_x) if trauma_x[i] != 0.0 || trauma_y[i] != 0.0]
if !isempty(valid_indices)
valid_x = trauma_x[valid_indices]
valid_y = trauma_y[valid_indices]
# Plot dots at trauma center locations
scatter!(ga, valid_x, valid_y, color = :red, markersize = 8, marker = :circle)
end
end
hidedecorations!(ga)
Our philosophy is that you shouldn’t be forced to memorize different plotting functions to visualize the data. You just need to inform the desired geometries to the viz and it will do the job.
Also, I would simply ignore the fallacy that GeoStats.jl is heavyweight.
People don’t realize how heavyweight GDAL and ArchGDAL.jl are. And how problematic their installation can be across platforms.
fig = Mke.Figure()
ax1 = Mke.Axis(fig[1, 1], aspect = Mke.DataAspect(), title = "b_1")
x = [coords(pt).x for pt in test.geometry]
y = [coords(pt).y for pt in test.geometry]
# for (dat, ax) in zip([test.b1, test.b2], [ax1, ax2])
for (dat, ax) in zip([test.b1], [ax1])
Mke.scatter!(ax, x, y, color = dat, markersize = 15)
end
fig
Mke.save("test2.png", fig)
There is an additional issue regarding units: Makie.jl doesn’t have full support for unitful axis across recipes, so we strip out the units of the coordinates before calling the built-in Makie.jl recipe (e.g., scatter) under the hood. At some point in the future, when all Makie.jl recipes become more consistent units-wise, we can review viz to avoid striping units.
Think of viz and viz! as the same thing as Mke.scatter or Mke.scatter!, but more general in the sense that it recognizes the types of geometries it receives.
Highly recommend reading the chapter of the book I shared above.
Ah yeah that makes sense, and works great. Thanks a lot! I’ve read that chapter a few times but I’m still pretty ignorant with this type of plotting (spoiled by ggplot2) so thanks for the guidance.
Once again, thanks a lot for your help and patience!
I will consider improving the chapter to include examples where the axis is manually constructed and populated with viz! calls. That might be useful to future readers.
You’ve been so helpful, I’m glad I was able to talk to you about the issues I was having. Feeling great about using Julia going forward since it clearly has a fantastic community
