Add another Choropleths tutorial.

Author: joa-quim

documentation/utilities/getbyattrib.md

@@ -22,6 +22,11 @@ NOTE: Instead of ``getbyattrib`` one can use instead ``filter`` (..., `index=fal
   Add as many as wished. If using two `kwargs` the second works as a condition. ``(..., NAME=("Antioquia", "Caldas"), feature_id=0)``
   means select all elements from ``Antioquia`` and ``Caldas`` that have the attribute `feature_id` = 0.
 
+  A second set of attributes can be used to select elements by region, number of polygon vertices and area.
+  For that, name the keyword with a leading underscore, e.g. `_region`, `_nps`, `_area`. Their values are
+  passed respectively a 4 elements tuple and numbers. E.g. `_region=(-8.0, -7.0, 37.0, 38.0)`, `_nps=100`, `_area=10`.
+  Areas are in square km when input is in geographic coordinates, otherwise squre data unites.
+
 - `attrib` or `att`: (OLD SYNTAX) A NamedTuple with the attribname, attribvalue as in `att=(NAME_2="value",)`.
   Use more elements if wishing to do a composite match. E.g. `att=(NAME_1="val1", NAME_2="val2")` in which
   case only segments matching the two conditions are returned.
@@ -41,6 +46,10 @@ Examples
 D = filter(D, NAME_2="Porto");
 ```
 
+```julia
+D = filter(D, _region="(-8.0, -7.0, 37.0, 38.0)", _nps=100);
+```
+
 See Also
 --------
 

documentation/utilities/polygonlevels.md

@@ -1,11 +1,11 @@
 # polygonlevels
 
 ```julia
-zvals = polygonlevels(D::GDtype, ids::Vector{String}, vals::Vector{<:Real}; kw...) -> Vector{Float64}
+zvals = polygonlevels(D::GDtype, ids::VecOrMat{String}, vals::Vector{<:Real}; kw...) -> Vector{Float64}
 
 or
 
-zvals = polygonlevels(D::GDtype, ids::Matrix{String}, vals::Vector{<:Real}; kw...) -> Vector{Float64}
+zvals = polygonlevels(D::GDtype, idvals::GMTdataset; kw...) -> Vector{Float64}
 ```
 
 Creates a vector with `zvals` to use in `plot` and where length(zvals) == length(D)
@@ -18,6 +18,10 @@ The elements of `zvals` are made up from the `vals`.
 
 - `vals`: is a vector with the numbers to be used in plot `level` to color the polygons.
 
+- `idvals`: is a GMTdataset with the `text` field containing the ids to match against the `ids` strings.
+            The first column of `id_vals` must contain the values to be used in `vals`. This is a comodity
+            function when both the `ids` and `vals` are store in a GMTdataset.
+
 - `attrib` or `att`: keyword to select which attribute to use when matching with contents of the `ids` strings.
 
 - `nocase` or `insensitive`: a keyword from `kw`. Perform a case insensitive comparision between the contents of

tutorials.md

@@ -45,6 +45,17 @@
     ~~~</a>~~~
   @@
 
+  @@box
+    ~~~<a class="boxlink" href="choropleths/choropleths/">~~~
+    @@title US Choropleths@@
+    @@box-content
+      ~~~
+      <img src="/tutorials/choropleths/tilelogo.png">
+      ~~~
+    @@
+    ~~~</a>~~~
+  @@
+
   @@box
     ~~~<a class="boxlink" href="dlakelan/GMTMaps/">~~~
     @@title @dlakelan choropleths@@

tutorials/choropleths.md

@@ -0,0 +1,3 @@
+# Choropleths
+
+{{list_folder_with_images choropleths}}

tutorials/choropleths/choropleths.md

@@ -0,0 +1,84 @@
+# Choropleths
+
+### What color is your State?
+
+This example shows an adaptation of the [Average color of World](https://erdavis.com/2021/06/26/average-colors-of-the-world/) examples.
+It differs slightly for the also shown US case likely due to the different data used. Here we retrieve
+the image data (Sentinel 2) from the [EOX WMS server](https://eox.at/) and use the year 2020.
+
+Note how the sates polygons are read directly from source, without any previous download, uncompressing file format conversions, etc...
+We will also filter the states polygons to use only those inside the US continental zone and limit them by number
+of points and with that dropping those polygons of small islands (smaller in area than 10 km^2).
+
+\begin{examplefig}{}
+```julia
+using GMT
+
+# Fetch the state polygons from the US Census Bureau 
+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)")
+```
+\end{examplefig}
+
+### Represent colors by average altitude
+
+We can also use the same polygons to represent some other variable, such as the average altitude of each state.
+We will use the [Earth Relief 06m](https://www.generic-mapping-tools.org/remote-datasets/earth-relief.html) dataset
+
+\begin{examplefig}{}
+```julia
+using GMT   # Hide
+D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");  # Hide
+
+G = gmtread("@earth_relief_06m");
+
+# Calculate the mean elevation per State.
+Dh = zonal_statistics(G, D, mean);
+
+# Create a color table for the mean elevation values.
+C = makecpt(range=(0, ceil(Dh.ds_bbox[2])), cmap=:bamako);
+
+viz(D, region=(-125,-66,24,50), proj=:guess, levels=Dh, cmap=C,
+               plot=(data=D,lw=0), title="Mean Elevation", colorbar=true)
+```
+\end{examplefig}
+
+The case above was relatively easy because the data was already prepared in a form that GMT could use.
+That is, the `zvals` vector, which together with the colorscale, determines the color of the polygons was
+already in the same order as the polygons themselves (in the GMTdataset `D`). However, often that is not the case,
+and we have the variable that contains the information that we want to _colorize_ in a different source and
+with a different order than that of the polygons in ``D``. So, we need to do a kind of _join_ operation.
+That can be done with _join_ functions or use the internal `polygonlevels` function that links the polygon
+names provide by a stored attribute in ``D`` and the values in a vector or a table that must also have some
+text information (a name) associated with each value. It is that we will do in the next example.
+
+
+\begin{examplefig}{}
+```julia
+using GMT   # Hide
+D = gmtread("/vsizip//vsicurl/https://www2.census.gov/geo/tiger/GENZ2024/shp/cb_2024_us_state_500k.zip");  # Hide
+
+# Read a simple text file that has population and State name, one per row.
+pop = gmtread(TESTSDIR * "assets/uspop.csv");
+
+# Use the polygonlevels function to get the values in the same order as the polygons in D.
+zvals = polygonlevels(D, pop, att="NAME") / 1e6;
+
+# Create a color table for the values in zvals.
+C = makecpt(zvals, auto=:r, reverse =true, cmap=:bamako);
+
+viz(D, region=(-125,-66,24,50), level=zvals, cmap=C, proj=:guess,
+    plot=(data=D,lw=0), title="Population (Millions)", colorbar=true)
+```
+\end{examplefig}