# Plotting in GeoStats

**URL:** <https://discourse.julialang.org/t/plotting-in-geostats/130770>\
**Category:** Geo\
**Created:** [July 16, 2025, 8:19am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770 "2025-07-16T08:19:33Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![beppi](https://avatars.discourse-cdn.com/v4/letter/b/b77776/32.png) [@beppi](https://discourse.julialang.org/u/beppi)\
**Post date:** [July 16, 2025, 8:19am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/1 "2025-07-16T08:19:33Z")

</div>

Hey,

I’ve seen the page here: [Visualization · GeoStats.jl](https://juliaearth.github.io/GeoStatsDocs/stable/visualization/)

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:

```julia
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):

```julia
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.

Thanks so much for any advice you guys can give!

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 16, 2025, 10:09am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/2 "2025-07-16T10:09:42Z")

</div>

Could you please share the result you are getting with viz? Why is it not matching the x and y coordinates?

---

<div class="post-metadata">

**Author:** ![technocrat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/technocrat/32/220947_2.png) [@technocrat](https://discourse.julialang.org/u/technocrat)\
**Post date:** [July 16, 2025, 10:18am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/3 "2025-07-16T10:18:57Z")

</div>

Assume a dataframe `df` with :geometry and :is\_trauma\_center, typeof IGeometry and Bool. It was read from shapefile with GeoDataTables.read()

```julia
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)

```

 ![dots](https://global.discourse-cdn.com/julialang/original/3X/0/9/09f173e436331e1c925e7000813c78632a03e957.jpeg)

With all respect to Dr. Hoffimann, GeoStats may be too heavyweight where its full capabilities aren’t needed.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 16, 2025, 10:47am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/4 "2025-07-16T10:47:55Z")

</div>

@beppi if all you need is plot the centroids, avoid all the mess of the previous answer and do:

```julia
viz(centroid.(geotable.geometry))

```

and if you need to display the geometries and their centroids:

```julia
viz(geotable.geometry)
viz!(centroid.(geotable.geometry), color="red")

```

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](https://github.com/JuliaEarth/GeoIO.jl/issues/167) their installation can be across platforms.

---

<div class="post-metadata">

**Author:** ![beppi](https://avatars.discourse-cdn.com/v4/letter/b/b77776/32.png) [@beppi](https://discourse.julialang.org/u/beppi)\
**Post date:** [July 17, 2025, 2:11am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/5 "2025-07-17T02:11:20Z")

</div>

This is the warped version using:

```julia
test = viz(
    test.geometry,
    color = test.b1,
    pointsize = 30
)
Mke.save("../test.png", test)

```

 ![Warped](https://global.discourse-cdn.com/julialang/original/3X/5/f/5f5569fd3d3af1dbffe76883e818e2f4238631db.png)

And this is using:

```julia
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)

```

 ![Not warped](https://global.discourse-cdn.com/julialang/original/3X/9/8/9895e9aa0106a7dee702e41c24fee16fa56e0fc5.png)

This would be solved in R/ggplot2 using `coord_equal`.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 17, 2025, 2:21am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/6 "2025-07-17T02:21:06Z")

</div>

I believe you just need to pass the `Mke.DataAspect()` the same way you passed it with the manual `Mke.scatter!` call.

The `viz` command is a normal Makie.jl recipe like `Mke.scatter`, so you can create the `Mke.Axis` as usual with the data aspect and call `viz` directly.

Explained in this chapter:

> **[2  Scientific visualization – Geospatial Data Science with Julia](https://juliaearth.github.io/geospatial-data-science-with-julia/02-geoviz.html)**

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 17, 2025, 2:24am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/7 "2025-07-17T02:24:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![beppi](https://avatars.discourse-cdn.com/v4/letter/b/b77776/32.png) [@beppi](https://discourse.julialang.org/u/beppi)\
**Post date:** [July 17, 2025, 2:36am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/8 "2025-07-17T02:36:11Z")

</div>

Thanks a lot, but I’m sorry I’m not sure how to pass the `Makie.DataAspect()` or `Makie.Axis()` instances to `viz`.

None of its arguments appear to be related to aspect ratio, or axes

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 17, 2025, 2:43am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/9 "2025-07-17T02:43:35Z")

</div>

Try this (adapted from your code):

```julia
ax = Mke.Axis(fig[1, 1], aspect = Mke.DataAspect(), title = "b_1")

viz!(ax, test.geometry, color=test.b1)

```

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.

---

<div class="post-metadata">

**Author:** ![beppi](https://avatars.discourse-cdn.com/v4/letter/b/b77776/32.png) [@beppi](https://discourse.julialang.org/u/beppi)\
**Post date:** [July 17, 2025, 2:48am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/10 "2025-07-17T02:48:40Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 17, 2025, 2:54am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/11 "2025-07-17T02:54:10Z")

</div>

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.

Thanks for raising the question!

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 17, 2025, 3:06pm UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/12 "2025-07-17T15:06:44Z")

</div>

@beppi added a new section to the chapter with more details about axis customization:

> **[2  Scientific visualization – Geospatial Data Science with Julia](https://juliaearth.github.io/geospatial-data-science-with-julia/02-geoviz.html#axis-customization)**

---

<div class="post-metadata">

**Author:** ![beppi](https://avatars.discourse-cdn.com/v4/letter/b/b77776/32.png) [@beppi](https://discourse.julialang.org/u/beppi)\
**Post date:** [July 18, 2025, 2:02am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/13 "2025-07-18T02:02:04Z")

</div>

Amazing, thanks so much!

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

---

<div class="post-metadata">

**Author:** ![technocrat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/technocrat/32/220947_2.png) [@technocrat](https://discourse.julialang.org/u/technocrat)\
**Post date:** [July 26, 2025, 3:22am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/14 "2025-07-26T03:22:10Z")

</div>

Julio, I truly meant no disrespect. I’ve no doubt that the backend of the GeoStats ecosystems improves greatly over GDAL. I had more in mind that among all the wealth I have found much in the way of advanced techniques, but I still haven’t found the geoscience _Hello, World_, a simple recipe to plot a choropleth map given one of the most widespread formats in existence, shapefiles. I might wish that Makie had a recipe for IGeometry, but it requires conversion to GeometryBasic first. I’m not yet up to rolling my own Makie convert recipe to do it without the routine above. But it may come to that, may the gods protect me.

So far, the only user-lightweight workaround that I have found is a PostGIS back-end to serve WKBs. In fact, it was that which mislead me to miss the conversion from IGeometry issues earlier.

What would cause me to jump feet-first into GeoStats is a short how-to of taking shapefiles and GeoJSON, to start, plot them using `vis` and embellish them, a bit, in Makie.

I’m old, not to put too fine a point on it. My degrees in geology and regional planning are over 50 years old. But I’m not convinced that flat maps have seen yet their obsolescence. I think GeoStats could find a broad audience for that use case.

As an inducement, I have a book in draft _Thematic Mapping with Julia_ aimed at a general audience of data scientists, journalists and general academic users who are assumed to have neither GIS nor programming backgrounds. My aim is to provide the minimum Julia necessary to produce high-quality thematic maps, such as choropleths. Anything that I can find to ease the learning experience on that is what I will feature.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 26, 2025, 7:00am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/15 "2025-07-26T07:00:15Z")

</div>

> [@technocrat](#):
>
> What would cause me to jump feet-first into GeoStats is a short how-to of taking shapefiles and GeoJSON, to start, plot them using `vis` and embellish them, a bit, in Makie.

```julia
using GeoStats
using GeoIO

import GLMakie as Mke

geotable = GeoIO.load("file.shp")

viewer(geotable)

```

?

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [July 26, 2025, 7:00am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/16 "2025-07-26T07:00:58Z")

</div>

One of the first things I did in Julia when I started as a complete novice was to make a chloropleth map. I recall the basics being pretty straightforward. I loaded shape files with geoIO, manipulated the GeoTable a bit and used Makie to plot. I had no skill, so it must have been comparatively easy!

I’m travelling, so can’t easily offer an mwe. However, this thread may help: [Scaling viz plots in GeoStats - #19 by TimG](https://discourse.julialang.org/t/scaling-viz-plots-in-geostats/112242/19)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 26, 2025, 7:49am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/17 "2025-07-26T07:49:14Z")

</div>

> [@technocrat](#):
>
> As an inducement, I have a book in draft _Thematic Mapping with Julia_ aimed at a general audience of data scientists, journalists and general academic users who are assumed to have neither GIS nor programming backgrounds

Notice that potencial readers coming from Python and R communities have no reason whatsoever to switch to Julia for GIS work if all they get are the same old GDAL+PROJ scripts and their inherent limitations.

Your book could target a modern audience, i.e., it could promote native Julia packages that bring something new to the mapping world. Our treatment of coordinate systems and geospatial domains is quite unique across all programming languages and you might be missing an opportunity here.

Good luck with the writing.

---

<div class="post-metadata">

**Author:** ![technocrat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/technocrat/32/220947_2.png) [@technocrat](https://discourse.julialang.org/u/technocrat)\
**Post date:** [July 26, 2025, 7:58am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/18 "2025-07-26T07:58:06Z")

</div>

![Screenshot 2025-07-26 at 12.43.40 AM](https://global.discourse-cdn.com/julialang/original/3X/0/f/0f6715351aa800e71ae0968ab847a60d99588bb4.png)

```julia-auto
julia> @time viewer(geotable)
102.084344 seconds (684.73 M allocations: 47.473 GiB, 5.39% gc time)

```

from that script in a clean environment on an M3 Ultra

using [https://www2.census.gov/geo/tiger/GENZ2024/shp/cb\_2024\_us\_state\_500k.zip](https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip)

I hadn’t found viewer(), but I’m still dissatisfied—it shouldn’t impute a z-axis missing from the data.

With the annoyance of the intermediate step of converting IGeometry to GeometryBasic, I get

 ![Screenshot 2025-07-26 at 12.54.24 AM](https://global.discourse-cdn.com/julialang/original/3X/6/c/6c3d033cfd04194484e3939352c6288be3def38e.jpeg)

with

```julia-auto
julia> @time poly!(ga, all_polygons, color= :white, strokecolor= :black, strokewidth=0.5)
  0.088717 seconds (17.18 k allocations: 17.771 MiB)
Poly{Tuple{Vector{Polygon{2, Float32, P, L, V} where {P<:AbstractPoint{2, Float32}, L<:(AbstractVector{<:GeometryBasics.Ngon{2, Float32, 2, P}}), V<:AbstractVector{L}}}}}

```

I’m sure it can be done, and I’m even sure that _I_ can do it, but I’m not sure that I will clearly be able to explain it to my readers.

If you like, we could take this offline?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 26, 2025, 8:57am UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/19 "2025-07-26T08:57:46Z")

</div>

If you read the GDSJL book I asked you to read a couple of times you will understand what is happening.

> [@technocrat](#):
>
> If you like, we could take this offline?

Only after you read all chapters of the book.

---

<div class="post-metadata">

**Author:** ![TimG](https://avatars.discourse-cdn.com/v4/letter/t/82dd89/32.png) [@TimG](https://discourse.julialang.org/u/TimG)\
**Post date:** [July 27, 2025, 3:23pm UTC](https://discourse.julialang.org/t/plotting-in-geostats/130770/20 "2025-07-27T15:23:08Z")

</div>

> [@technocrat](#):
>
> I still haven’t found the geoscience _Hello, World_, a simple recipe to plot a choropleth map given one of the most widespread formats in existence, shapefiles.

Back home now. Here is a quick MWE based on the shapefile you linked. This file includes boundaries for several areas outside the contiguous USA (eg Guam), so I’ve excluded these. The shapefile uses the ‘NAD83’ reference system, and I’ve stuck to this in the MWE.

```julia-auto
using GeoIO
using GeoStats

import CairoMakie as Mke

gt = GeoIO.load("cb_2024_us_state_500k.shp") |> Filter(row -> row.NAME ∉ ["Alaska", "Hawaii", "Guam", "Puerto Rico", "American Samoa", "United States Virgin Islands", "Commonwealth of the Northern Mariana Islands"])

fig = Mke.Figure()
ax = Mke.Axis(fig[1, 1], title = "US 'Lower 48' States")
viz!(ax, gt.geometry, color=1:nrow(gt), segmentcolor="black", showsegments=true, segmentsize=0.3f0)
ax.aspect = Mke.DataAspect()
Mke.display(fig)

```

This produces a simple chloropleth map straight from the shapefile:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/a/da23d700a76a3ff85ac612d942d3cb4d05a52f7c.png)

This is obviously a much simpler task than the OP, which @juliohm solved, but it shows how easy it is to create a “hello world” map.

[Next page](https://discourse.julialang.org/t/plotting-in-geostats/130770.md?page=2)