# Let's Choroplethize the US
**URL:** https://discourse.julialang.org/t/lets-choroplethize-the-us/131254
**Category:** Geo
**Tags:** geo, gmt
**Created:** [July 31, 2025, 4:11pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254 "2025-07-31T16:11:08Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [July 31, 2025, 4:11pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/1 "2025-07-31T16:11:08Z")
Hi,
Making choropleth maps of US seems to a popular request/example here. Advanced solution often lack, IMO, the necessary map quality. In these last days before vacations I spent some time preparing these [examples](https://www.generic-mapping-tools.org/GMTjl_doc/tutorials/choropleths/choropleths/) that, I hope show the simplicity and mainly the map quality of the GMT.jl solutions
One of such examples
```julia-auto
using GMT
# Fetch the state polygons from the US Census Bureau (no need to download-uncompress-file_format_change. Just read it)
D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");
# Filter to keep only the continental US states (except Alaska), with an area larger than 10 km^2.
Df = filter(D, _region=(-125,-66,24,50), _area=10);
# Fetch the Sentinel 2 image that is used to calculate the average color. Restrain to pixel size of 2000 m.
wms = wmsinfo("http://tiles.maps.eox.at/wms?");
img = wmsread(wms, layer=4, region=(-125,-66,24,50), pixelsize=2000);
# Calculate the average color per State.
colorzones!(Df, median, img=img)
viz(Df, region=img, proj=:guess, plot=(data=Df, lw=0), title="Sate Color (median)")
```
---
### Author:  [@Jake](https://discourse.julialang.org/u/Jake)
#### Post date: [July 31, 2025, 4:48pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/2 "2025-07-31T16:48:58Z")
This is my first time using GMT. I received the error:
```julia-auto
julia> using GMT
julia> D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");
ERROR: Must select the input data type (grid, image, dataset, ogr, cmap or ps)
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] gmtread(_fname::String; kwargs::@Kwargs{})
@ GMT C:\Users\jakez\.julia\packages\GMT\W1WH8\src\gmtreadwrite.jl:176
[3] gmtread(_fname::String)
@ GMT C:\Users\jakez\.julia\packages\GMT\W1WH8\src\gmtreadwrite.jl:75
[4] top-level scope
@ REPL[6]:1
julia>
```
I am using 1.11.6 on Win11.
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [July 31, 2025, 4:50pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/3 "2025-07-31T16:50:58Z")
Sorry, forgot to warn that to replicate this example one needs the latest GMT.jl (version 1.31.1)
---
### Author:  [@Jake](https://discourse.julialang.org/u/Jake)
#### Post date: [July 31, 2025, 5:18pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/4 "2025-07-31T17:18:04Z")
That was interesting. It was my first time loading GMT and it loaded 1.30.0, I did an `]up` and it updated to 1.30.1.
I got past the first line and then the third line errored. When I add quotes it seems to work however.
```julia-auto
julia> wms = wmsinfo(http://tiles.maps.eox.at/wms?);
ERROR: ParseError:
# Error @ REPL[11]:1:20
wms = wmsinfo(http://tiles.maps.eox.at/wms?);
# └┘ ── not a unary operator
Stacktrace:
[1] top-level scope
@ none:1
julia> wms = wmsinfo("http://tiles.maps.eox.at/wms?");
julia>
```
Everything then seems to run except the last line which errors:
```julia-auto
julia> viz(Df, region=img, proj=:guess, plot=(data=Df, lw=0), title="Sate Color (median)")
Access is denied.
psconvert [ERROR]: System call [@"C:\Users\jakez\.julia\artifacts\cca2b0693e0ae8df2daca37d0779bb7209f32c5f\bin\gs.exe" -q -dNOSAFER -dNOPAUSE -dBATCH -sDEVICE=bbox -DPSL_no_pagefill -dMaxBitmap=2147483647 -dUseFastColor=true "C:\Users\jakez\AppData\Local\Temp/GMTjl_jakez.ps" 2> "./psconvert_27032c.bb"] returned error 1.
ERROR: Something went wrong when calling the module. GMT error number = 79
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] gmt(::String)
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\gmt_main.jl:166
[3] showfig(d::Dict{Symbol, Any}, fname_ps::String, fname_ext::String, opt_T::String, K::Bool, fname::String)
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\common_options.jl:4194
[4] finish_PS_module_barr_last(d::Dict{…}, cmd::Vector{…}, fname::String, fname_ext::String, opt_extra::String, output::String, opt_T::String, K::Bool, P::Nothing)
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\common_options.jl:4558
[5] finish_PS_module(d::Dict{…}, cmd::Vector{…}, opt_extra::String, K::Bool, O::Bool, finish::Bool, args::Vector{…})
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\common_options.jl:4496
[6] prep_and_call_finish_PS_module
@ C:\Users\jakez\.julia\packages\GMT\acUHa\src\common_options.jl:4407 [inlined]
[7] prep_and_call_finish_PS_module(d::Dict{…}, cmd::Vector{…}, opt_extra::String, K::Bool, O::Bool, finish::Bool, arg1::Vector{…}, arg2::Nothing, arg3::Nothing, arg4::Nothing)
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\common_options.jl:4405
[8] _common_plot_xyz(cmd0::String, arg1::Vector{…}, caller::String, O::Bool, K::Bool, is3D::Bool, d::Dict{…})
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\psxy.jl:253
[9] #common_plot_xyz#924
@ C:\Users\jakez\.julia\packages\GMT\acUHa\src\psxy.jl:13 [inlined]
[10] common_plot_xyz
@ C:\Users\jakez\.julia\packages\GMT\acUHa\src\psxy.jl:11 [inlined]
[11] plot(arg1::Vector{…}; first::Bool, kw::@Kwargs{…})
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\plot.jl:122
[12] plot
@ C:\Users\jakez\.julia\packages\GMT\acUHa\src\plot.jl:101 [inlined]
[13] imshow(arg1::Vector{…}, x::Vector{…}, y::Vector{…}; kw::@Kwargs{…})
@ GMT C:\Users\jakez\.julia\packages\GMT\acUHa\src\imshow.jl:86
[14] imshow (repeats 2 times)
@ C:\Users\jakez\.julia\packages\GMT\acUHa\src\imshow.jl:36 [inlined]
[15] top-level scope
@ REPL[15]:1
Some type information was truncated. Use `show(err)` to see complete types.
```
---
### Author:  [@TimG](https://discourse.julialang.org/u/TimG)
#### Post date: [July 31, 2025, 5:24pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/5 "2025-07-31T17:24:48Z")
Just added GMT into a clean environment:
```julia-auto
julia> using GMT
julia> D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");
ERROR: Must select the input data type (grid, image, dataset, ogr, cmap or ps)
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] gmtread(_fname::String; kwargs::@Kwargs{})
@ GMT C:\Users\tim\.julia\packages\GMT\hx804\src\gmtreadwrite.jl:176
[3] gmtread(_fname::String)
@ GMT C:\Users\tim\.julia\packages\GMT\hx804\src\gmtreadwrite.jl:75
[4] top-level scope
@ REPL[2]:1
julia>
pkg> st
Status `C:\Users\tim\OneDrive\Documents\Julia\GMT\Project.toml`
[5752ebe1] GMT v1.31.0
```
It seems the latest version isn’t the one that gets downloaded [Just noticed previous posts on this, too!]
---
### Author:  [@TimG](https://discourse.julialang.org/u/TimG)
#### Post date: [July 31, 2025, 5:35pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/6 "2025-07-31T17:35:05Z")
Now updated to v1.31.1
Then:
```julia-auto
julia> wms = wmsinfo("http://tiles.maps.eox.at/wms?")
serverURL: http://tiles.maps.eox.at/wms?
OnlineResource: http://tiles.maps.eox.at/?
version: 1.1.1
request: GetMap
layernames: ["hydrography", "coastline_3857", "s2cloudless-2021_3857", "overlay_base", "s2cloudless-2023_3857", "coastline_black", "s2cloudless-2017_3857", "s2cloudless-2019_3857", "overlay", "s2cloudless-2020", "s2cloudless-2021", "osm_3857", "s2cloudless-2022", "s2cloudless-2023", "s2cloudless-2024", "s2cloudless_3857", "streets", "terrain_3857", "osm", "blackmarble", "overlay_base_bright_3857", "bluemarble_3857", "overlay_base_3857", "s2cloudless", "graticules", "terrain", "overlay_bright_3857", "overlay_3857", "s2cloudless-2020_3857", "magnetic_graticules", "terrain-light", "overlay_base_bright", "bluemarble", "s2cloudless-2022_3857", "s2cloudless-2024_3857", "streets_3857", "s2cloudless-2018_3857", "overlay_bright", "coastline", "blackmarble_3857", "hydrography_3857", "s2cloudless-2017", "s2cloudless-2018", "s2cloudless-2019", "terrain-light_3857"]
layer: 45 Layers. Use layer[k] to see the contents of layer k
julia> img = wmsread(wms, layer=4, region=(-125,-66,24,50), pixelsize=2000);
┌ Warning: This request returned an empty image.
└ @ GMT C:\Users\tim\.julia\packages\GMT\acUHa\src\extras\webmapserver.jl:156
```
Edit: If I cut and paste into a file and “run without debugging” in VS-code, this works fine. Only a problem line by line in a VS-Code REPL.
Edit2: Went away and came back and now it works as it should. Don’t know what happened. Sorry for the noise!
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [July 31, 2025, 5:53pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/7 "2025-07-31T17:53:31Z")
Hmm,
- The installed version. No idea why. I just created the 1.31.1 version before the post. Maybe registry is not updating with the required frequency.
- The error `wms = wmsinfo(http://tiles.maps.eox.at/wms?);└┘ ── not a unary operator`.
My bad, the copy past must have removed the quotes. This same [example](https://www.generic-mapping-tools.org/GMTjl_doc/tutorials/choropleths/choropleths/#what_color_is_your_state) in the Docs have them.
- The `Warning: This request returned an empty image.` error. Honestly don’t know what to do with this to make it robust. Apparently the WMS server does not always treturn the same response. For me it shows:
layernames: [“terrain”, “overlay\_bright\_3857”, “coastline\_3857”, “s2cloudless-2020\_3857”,
but for you “s2cloudless-2020\_3857” shows up as layer 3 (but 4 for me and Github where docs are built). The solution is to change `layer=4` to `layer=3`.
- The `ghostscript` error (@Jake). Also `dono`, but I notice the first error message is `Access is denied.`, which continues to be strange as the images output from GMT are written in `/tmp` where everybody is supposed to have read/right access. See if you have a file called /tmp/GMTjl\_(your\_user\_name).ps
---
### Author:  [@Jake](https://discourse.julialang.org/u/Jake)
#### Post date: [July 31, 2025, 8:08pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/8 "2025-07-31T20:08:15Z")
> [@joa-quim](#):
>
> See if you have a file called /tmp/GMTjl\_(your\_user\_name).ps
Since it is a windows machine there is no folder called c:\tmp. That is probably why the error.
Edit: Upon further investigation, it is trying to put the file in “C:\Users\jakez\AppData\Local\Temp/GMTjl\_jakez.ps” according to the error message. The folder exists, has full permissions, but the file does not exist.
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [July 31, 2025, 8:29pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/9 "2025-07-31T20:29:54Z")
Yes, the destiny is the content of `GMT.PSname[1]`. Check the `img` is not nothing. You might have suffered the same thing as TimG (see above, and my reply). What does this prints?
```julia-auto
wms = wmsinfo("http://tiles.maps.eox.at/wms?")
```
---
### Author:  [@Jake](https://discourse.julialang.org/u/Jake)
#### Post date: [July 31, 2025, 8:37pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/10 "2025-07-31T20:37:45Z")
> [@joa-quim](#):
>
> What does this prints?
Nothing
```julia-auto
julia> wms = wmsinfo("http://tiles.maps.eox.at/wms?");
julia>
```
Without the semicolon at the end of the line it prints:
```julia-auto
julia> wms = wmsinfo("http://tiles.maps.eox.at/wms?")
serverURL: http://tiles.maps.eox.at/wms?
OnlineResource: http://tiles.maps.eox.at/?
version: 1.1.1
request: GetMap
layernames: ["terrain", "overlay_bright_3857", "coastline_3857", "s2cloudless-2020_3857", "magnetic_graticules", "terrain-light", "overlay_base_bright", "overlay_base", "s2cloudless-2022_3857", "s2cloudless-2024_3857", "s2cloudless-2018_3857", "coastline", "s2cloudless-2020", "s2cloudless-2021", "osm_3857", "s2cloudless-2022", "s2cloudless-2023", "s2cloudless-2024", "s2cloudless_3857", "terrain_3857", "overlay_base_bright_3857", "hydrography", "graticules", "overlay_3857", "bluemarble", "s2cloudless-2021_3857", "s2cloudless-2023_3857", "coastline_black", "s2cloudless-2017_3857", "streets_3857", "overlay_bright", "s2cloudless-2019_3857", "overlay", "blackmarble_3857", "hydrography_3857", "s2cloudless-2017", "streets", "s2cloudless-2018", "s2cloudless-2019", "osm", "blackmarble", "terrain-light_3857", "bluemarble_3857", "overlay_base_3857", "s2cloudless"]
layer: 45 Layers. Use layer[k] to see the contents of layer k
julia>
```
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [July 31, 2025, 8:49pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/12 "2025-07-31T20:49:22Z")
You have to do the commands one by one to see where it fails. Does this get you an image that you can vizualize with `viz(img)`?
`img = wmsread(wms, layer=4, region=(-125,-66,24,50), pixelsize=2000)`
But if this continues it’s better that you open an issue and we’ll continue there to not bother every body with this debug attempts.
---
### Author:  [@Jake](https://discourse.julialang.org/u/Jake)
#### Post date: [July 31, 2025, 8:59pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/13 "2025-07-31T20:59:35Z")
No errors until I get to the viz line. I will open an issue.
---
### Author:  [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)
#### Post date: [August 1, 2025, 12:24am UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/14 "2025-08-01T00:24:51Z")
I just tried this with 1.31.1 and got the same error
```julia-auto
julia> D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");
ERROR: Must select the input data type (grid, image, dataset, ogr, cmap or ps)
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] gmtread(_fname::String; kwargs::@Kwargs{})
@ GMT ~/.julia/packages/GMT/acUHa/src/gmtreadwrite.jl:177
[3] gmtread(_fname::String)
@ GMT ~/.julia/packages/GMT/acUHa/src/gmtreadwrite.jl:75
[4] top-level scope
@ REPL[11]:1
```
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [August 1, 2025, 1:07am UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/15 "2025-08-01T01:07:33Z")
> [@pdeffebach](#):
>
> `D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");`
OK, make it as it wishes. Change to
`D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip", ogr=true);`
But next lines will error if version is not 1.31.1
---
### Author:  [@technocrat](https://discourse.julialang.org/u/technocrat)
#### Post date: [August 1, 2025, 3:15am UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/16 "2025-08-01T03:15:05Z")
@joa-quim and I have exchanged friendly correspondence about my own struggles with a “hello, world” under macOS (I’ve had trouble with the GeoStats suite as well). Both are probably due to deficiencies in my understanding. What I have found doable was Makie/GeoMakie, to do the attached.
That said, for state-level work, I’m coming to the belief that choropleths are a flawed approach because of the perceptual over-emphasize on area. Tiny Washington, DC and vasty Wyoming have roughly the same population and wield the same Electoral College, is the most obvious example.
That’s why I have StateBins.jl pending registration, shamelessly ripped off (well, inspired by) the `R` package of the same name by the estimable Bob Rudis.
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [August 1, 2025, 2:03pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/17 "2025-08-01T14:03:32Z")
Shouldn’t be too hard to modify the [bubblechart](https://www.generic-mapping-tools.org/GMTjl_doc/documentation/modules/bubblechart/index.html#bubblechart) function to do that as well (and projected if wished). Maybe after vacations.
---
### Author:  [@technocrat](https://discourse.julialang.org/u/technocrat)
#### Post date: [August 1, 2025, 9:49pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/18 "2025-08-01T21:49:55Z")
Cool, thanks. I think I’d use centroids without projections, since we’re not concerned with preserving areal dimensions in a cartogram.
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [August 1, 2025, 10:45pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/19 "2025-08-01T22:45:59Z")
Yes, I get that but for GMT.jl, projected or not is for free and I was (vaguely) thinking in a tool not particularly tailored for a specific Earth region.
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [August 31, 2025, 10:34pm UTC](https://discourse.julialang.org/t/lets-choroplethize-the-us/131254/20 "2025-08-31T22:34:53Z")
I ended up adding [one version](https://www.generic-mapping-tools.org/GMTjl_doc/tutorials/choropleths/choropleths/#choropleths_by_symbol_color) based on centroids. Not specific to any country, which may lead to some overlapings.
