Let's Choroplethize the US

,

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 that, I hope show the simplicity and mainly the map quality of the GMT.jl solutions

One of such examples

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)")

2 Likes

This is my first time using GMT. I received the error:

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.

Sorry, forgot to warn that to replicate this example one needs the latest GMT.jl (version 1.31.1)

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

Just added GMT into a clean environment:

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!]

Now updated to v1.31.1

Then:

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!

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 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

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.

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?

wms = wmsinfo("http://tiles.maps.eox.at/wms?")

Nothing

julia> wms = wmsinfo("http://tiles.maps.eox.at/wms?");

julia>

Without the semicolon at the end of the line it prints:

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>

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.

No errors until I get to the viz line. I will open an issue.

I just tried this with 1.31.1 and got the same error

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

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

@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.

5 Likes

Shouldn’t be too hard to modify the bubblechart function to do that as well (and projected if wished). Maybe after vacations.

1 Like

Cool, thanks. I think I’d use centroids without projections, since we’re not concerned with preserving areal dimensions in a cartogram.

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.

1 Like