Creating examples of GMT based geographic data display

[dlakelan]

I'd like to create an example Tutorial for using GMT.jl in Julia to do displays of social and scientific data, targeted at people who know something about working with datasets but generally are very naive about terminology, file formats, and specifics of mapping jargon etc.

Some examples of the kinds of things I'm imagining people in the target audience might want to create:

The audience would be people who understand how to use CSV, DataFrames, DataFramesMeta, merging, joining, data munging, regression, maybe simulations etc, but would like to display their data in the context of a spatial display, a map, or a 3D model or something along those lines.

Suggestions are welcome here for specific projects that can be tackled, either using specific datasets, or reproducing specific kinds of visualizations, or just general suggestions of "it'd be nice to see how to..."

@joa-quim welcome and thanks for your help!

[joa-quim]

Quick thoughts (@Pwessel)

• 1 Looks doable but need to know what type of data was used to create the image
• 2 A challenge
• 3 & 4 Doable
• 5 Looks simple (aside from the text positioning)
• 6 Trivial ... if we know how to create such an image
• 7 May be trivial also if we know where to get that image that looks the output of a meteorological model forcast.

[dlakelan]

So I don't know if we have to follow those exactly but they were more to imagine the sorts of things that might be of interest.

I suggest we start with the Choropleth of some Census bureau data since the data is usually easy to acquire relatively.

I will find a suitable county dataset maybe income or age or something. I would like to start off with an explanation of the geo files that are needed, so for example the county shapefiles and the census info... Then describe for example what is in a shapefiles and how you can build a dataset that provides the measurements for each county, how GMT interprets that dataset, how to build a color palette, how to put it all together, what kinds of projections you can use etc etc

For something like the first map with one dot per X people we might be able to get Census PUMS/PUMA data (public use microdata area), each area represents close to 100k people in the US, so if we could find the centroid of each area we could plot something like that map.

[joa-quim]

Ok, if you have true data is better. Won't have much time to look at this today but tomorrow afternoon should have some.

[joa-quim]

My quick attempt with random data worked better than I thought so I'll leave it for you to play with when you have true data

countymap = gmtread("cb_2018_us_county_5m/cb_2018_us_county_5m.shp");
        This file has islands (holes in polygons).
        Use `gmtread(..., no_islands=true)` to ignore them.

# Extract the county names
ids = [countymap[k].attrib["NAME"] for k=1:length(countymap)];

# Invent the data to be represented
vals = rand(length(countymap));

# Create the variable that associates a value from `vals`  to each of the `countymap` polygon
zvals = GMT.make_zvals_vec(countymap, ids, vals, attrib="NAME");

# Create a colormap. Here between [0,1] becuase we used rand() above
C = makecpt(range=(0,1),);

# Quick display
imshow(countymap, level=zvals, cmap=C, close=true, fill="+z", pen=0.5, region=(-125,-65,24,50), proj=:guess, colorbar=true)

[dlakelan]

Great, I will start with that. Can you tell me what kind of object is "countymap"? what sorts of things can we do with it in Julia? Explaining this will help those who don't understand much about mapping! (like me!)

when we call make_zvals_vec it seems to be doing a kind of "join" between values in "vals" and something inside "countymap" is that correct? and the join is based on some aspect of countymap (in particular "NAME") having the same value as the corresponding value in "ids"?

[dlakelan]

Also, can we use a custom colormap created in Julia? Suppose for example I want something created in Colors.jl like: http://juliagraphics.github.io/Colors.jl/stable/colormapsandcolorscales/ or something premade from ColorSchemes.jl: https://juliagraphics.github.io/ColorSchemes.jl/stable/

[joa-quim]

This is going to be interesting between what we think the users should know and what they actually do. Don't have time now to expand (will do in 3-4 hours or tomorrow). Quick answers

See the manual https://www.generic-mapping-tools.org/GMT.jl/dev/types/ for the GMT.jl types, namely the GMTdataset and vectors of them (the countymap).

See the online help of make_zvals_vec (? make_zvals_vec)

[dlakelan]

Great, thanks. I really appreciate any time you can give I know it's quite limited!

[joa-quim]

Hmm, not sure where to start with this so let me sate the basics and go a bit beyond that.

We use x,y[,z] to plot lines, points (or any other symbols) or polygons. All plotting packages can do that by providing the data in files or Mx2[3] arrays. GMT is no exception and accepts input from text files with data as x y [z] or other field separators (e.g. x,y[,z]) in what is generically known as the comma separated values (or CSV).

However, the above is a pure Cartesian that suites many usage cases but carries no information about the metadata. Imagine, and that is a very common situation, that the x,y are points in a projection system, UTM for instances. We need then to specify the projection parameters so that the software knows how to plot the data and eventually deals with data in other projection systems. We need then to store that info in the metadata. Think also that data contains the polygons delineating US Counties and we want to keep a trace of their names, area, centroids, zip code, etc. The most traditional format that stores this kind of information is the ESRI shapefiles (but several database formats may also be used) and it stores it in what is known as the attributes table. Those attributes can be anything but I'll stop here I think it's the limits of my knowledge. When GMT reads one of these files via gmtread (which lowers down to use the most used open source code, GDAL) it returns a GMTdataset object that holds the data, it's geometry type (point, line, polygon) and the eventual attributes. "NAME" in the choropleth example above is one of those attributes.

I'll stop here for now and let you come with more questions (we still have to talk about grid and image formats that are the base map data).

[joa-quim]

Sorry, didn't answer to the colormaps question. Don't know how difficult is to use the colormaps from Colors.jl, probably not much but never looked at it because we have many of them and overall a quite long list in GMT. See
https://docs.generic-mapping-tools.org/dev/cookbook/cpts.html

[dlakelan]

Ok, so here are some more real-world consideration. As you know, when working with real data things get messy!

I'm using the dataset at https://www2.census.gov/programs-surveys/popest/datasets/2010-2020/counties/totals/co-est2020-alldata.csv which has county population estimates.

In the US county names are not unique, but we have FIPS standards (federal information processing standards) which assign a unique code to each state and to each county in the state... so there are two fields we need to match on.

I read in the counties shape file:

counties = gmtread("data/cb_2018_us_county_5m.shp")

and I see that counties[1].attrib is a dictionary and I can extract the STATEFP and COUNTYFP values to make a dataset:

dfc = DataFrame(STATE = map(x->x.attrib["STATEFP"],counties),COUNTY=map(x -> x.attrib["COUNTYFP"],counties))
joineddata = leftjoin(dfc,cpop,on= [:STATE => :STATEFIPS,:COUNTY=>:CFIPS],makeunique=true)

Now there are several non-matching lines towards the end of the dataset... with missing values...perhaps they correspond to some outlying territories like Guam or something I don't know. If I try to use them gmt doesn't like it:

cptval = makecpt(range=(0,11e6))

imshow(counties,level=joineddata.POPESTIMATE2020,cmap=cptval,close=true,fill="+z",pen=0.25,region=(-125,-65,24,50),proj=:guess,colorbar=true)

So I try eliminating the missingness:

countiesnom = counties[.! ismissing.(joineddata.POPESTIMATE2020)]
joineddatanom = joineddata[.! ismissing.(joineddata.POPESTIMATE2020),:]
imshow(countiesnom,level=Vector{Int64}(joineddatanom.POPESTIMATE2020),cmap=cptval,close=true,fill="+z",pen=0.25,region=(-125,-65,24,50),
    proj=:guess,colorbar=true)

This results in a plot... but it's clearly wrong. For example Los Angeles County is in the lower left of this map, and should have 9M people or so, and yet it is colored dark blue meaning basically less than 1M and some random counties in the middle of the country are colored crazy high populations.

[dlakelan]

and yes I do realize I didn't do the make_zvals_vec method, but partly that's because I don't understand what it does, and partly because it feels like maybe the manipulation with the DataFrames join, which feels way way more natural, did something similar?

[dlakelan]

Ok @joa-quim I've pushed a new version with GMTMaps.jmd which includes the code I'm using to make the choropleth. You can see what I'm up to more explicitly. If you want to try to run it, just:

using Weave
weave("GMTMaps.jmd")

[joa-quim]

Tried it but I'm getting an error at this step (and other errors in consequence of it)

cpop = cpop[cpop.COUNTY != "000",:]
ERROR: ArgumentError: invalid row index of type Bool
Stacktrace:
 [1] DataFrameRow
   @ C:\Users\joaqu\.julia\packages\DataFrames\MA4YO\src\dataframerow\dataframerow.jl:98 [inlined]
 [2] getindex(df::DataFrame, rowind::Bool, #unused#::Colon)
...

Also, please remove the reference to the header in the description because I'm going to change a bit of what's get stored when reading an OGR format (shapefiles is only one of them).

The reason is that I started to store the aspatial attributes as a comma separated string in header at the time of the make_zvals_vec() first version, but later created the attrib dictionary field that is much easier to access. So, being repeated I'll remove the attributes from the header. This header is an heritage from the GMT command line that started to be a place where were stored ASCII file information but also evolved to a place holder where one can store plotting GMT options on a per segment basis ... but this is another story.

EDIT: needs to be cpop = cpop[cpop.COUNTY .!= "000",:]

[joa-quim]

So, the problem is the order of the values in the level variable. When data arrives to GMT it will paint each polygon with the color determined by the corresponding value in level (after consulting the colormap ofc). The issue here is that the shapefile with the counties and the csv with the population do not have the data by the same order. It starts to have it but then somewhere in the middle orders get different.

julia> joineddata[3450,:]
DataFrameRow
  Row │ STATE   COUNTY  STATE_1  COUNTY_1  STNAME     CTYNAME         POPESTIMATE2020
      │ String  String  Int64?   Int64?    String31?  String?         Int64?
──────┼───────────────────────────────────────────────────────────────────────────────
 3450 │ 13      307          13       307  Georgia    Webster County             2595

vs

counties[3450].attrib
Dict{String, String} with 9 entries:
  "LSAD"     => "06"
  "COUNTYFP" => "095"
  "GEOID"    => "13095"
  "COUNTYNS" => "00351259"
  "AWATER"   => "15352849"
  "ALAND"    => "851183464"
  "AFFGEOID" => "0500000US13095"
  "STATEFP"  => "13"
  "NAME"     => "Dougherty"

There is some issue when reading the shapefile that is loosing part of the name when it's a multi words name, but that's something to address latter.
One of the roles of the make_zvals_vec() function is to make sure that the order of the level vector (the output of make_zvals_vec) is the same as the polygons in the GMTdataset vector (the counties in this case). But you are also right that as is make_zvals_vec() is not useful because it only checks one condition and in this case we need to check for two, the county and state codes. I liked the DataFrame approach. It's more elegant and we can solve the other minor issues with the data types, but it must be able to solve the order issue.

[dlakelan]

Ok, so it's easy to reorder the dataframe, perhaps it loses order when I do the join. Let me check that out.

[dlakelan]

Great, so now I did this:

dfc = DataFrame(STATE = map(x->x.attrib["STATEFP"],counties),COUNTY=map(x -> x.attrib["COUNTYFP"],counties),ORDER=1:length(counties))


joineddata = @chain leftjoin(dfc,cpop,on= [:STATE => :STATEFIPS,:COUNTY=>:CFIPS],makeunique=true) begin
@orderby(:ORDER)
end

And then when I plot I get a more reasonable map (this is using logarithm of the population so the colors aren't so starkly different)

[dlakelan]

Without the logarithm transform... it's clear the LA county should be in charge of everything :-)

[dlakelan]

Ok, so this is fantastic because it shows a clear workflow that will be understandable to data analysts and uses general tools such as DataFrames for achieving the necessary data munging. Thank you so much for help with this first example!

I think an interesting next question perhaps is how to calculate population density (ie. take the population and divide by the area of the polygon). I assume GMT probably has a way to calculate the polygon area?

[joa-quim]

Areas (I knew this would come quickly)
One can compute them with geomarea that uses GDAL and reports
the area of the feature in square units of the spatial reference system in use.

geomarea(counties[1])
0.1507503541454951

or use gmtspatial

gmtspatial(counties[1], area=true)
...
1×3 Matrix{Float64}:
 -83.6008  39.1847  0.15075

area is third value and it's the same as the other value. However, according to the GMT manual (sorry, no extnded GMT.jl manual for this) the default units should be square km, but it's not (Paul @Pwessel). To get them in km^2

gmtspatial(counties[1], area="k")
...
1×3 Matrix{Float64}:
 -83.6009  39.1847  1446.7

or just

 gmtspatial(counties, area="k")

to compute area of all polygons at once

[dlakelan]

Ok, I've done an example where it calculates density, which is a very interesting plot. but when I do gmtspatial(counties,area="k") I get back a length-1 vector of GMTdataset. So then I just immediately take the 1st element. What is the reason it gives back a vector of datasets rather than just a dataset? Is there some other situation where the vector makes sense?

Also the dataset seems to have 3 columns as you say. It looks like the first two are lat/long values? what do they correspond to? Centroids of the polygons?

[joa-quim]

Yes, first two columns are the x/y (lon/lat in this case) coordinates of the centroid.
Regarding the fact that it returns a vector even when there is only one segment is something that annoys me a bit too but is generalized thing that happens at the translation from the data that comes from C (the GMT lib) and the Julia type. I have been hesitant to try to change that, which now would impact in lots of places in the GMT.jl code. Though it shouldn't be difficult to change that in this particular case.

[dlakelan]

Also, suppose I want to save the image to a particular filename. Right now imshow is outputting to /tmp/GMTjl_tmp.png and overwriting the previous thing when I do a new plot. Suppose I want to plot to particular output files say for inclusion into a document or whatever? I imagine there is a keyword argument but I'm not seeing anything obvious in the docs to tell me which one, or perhaps I should be calling a different function to output to a file rather than show it in a window?

Finally, is there a way to make the image show up in the output of the "weave" or in a Jupyter notebook etc?

[joa-quim]

Docs is an issue in GMT.jl, I know. A quite hard task as the GMT docs are big and one needs to adapt them but all plotting modules have this in the quick online docstrings

savefig | figname | name :: [Type => Str]

Save the figure with the figname=name.ext where ext chooses the figure format (e.g. figname="name.png")

Finally, is there a way to make the image show up in the output of the "weave" or in a Jupyter notebook etc?

Still didn't try weave but Jupyter and Pluto/Neptune display the images automatically.
This is a Jupyter nootebook

[joa-quim]

Ah, same fig as mine with a modified version of the make_zvals_vec function.
Note, I don't like this function name and that's why I have not yet exported it. It's just not a name a user would get any feeling of what it does. If you have suggestions for other names, I'm all ears.

Don't mistake me, I like the DataFrames solution and I think it will be a more natural for people used to it but I'm too biased into the internal workings and prefer something more straight to the goal.

Slightly different from yours because my colormap is: cptval = makecpt(range=(log(100),log(10e6)),C=:plasma); to avoid th backs that fall of the colormap minimum.

[joa-quim]

Maybe the simplest is to look at the code. The reason for this function is simply that I didn't know how to make what you showed me with DataFrames. And you are right that it intended to deal with the case of repeated polygons with the same identifiers (for examples a country with islands) but even that is not been taken care in latest versions of that function. I wonder if that is not the reason why you ended up with some missing in the population DataFrame, which ended up in deleting the corresponding polygons. Regarding this, the strategy should be to replace missing by NaN and not deleting the polygons. With NaNs the polygons are painted with the NaN color of the colormap (cpt) that is white by default but configurable.

BTW, do you know to identify variables with type Vector{Union{Missing, Int64}}? I could then do replacement under the hood and also accept that type in the legend option. The basic issue here is that I cannot let go those exotic types (plus missing) into the C lib that has no idea what that is.

[dlakelan]

BTW, do you know to identify variables with type Vector{Union{Missing, Int64}}?

julia> typeof(myvar)
Vector{Union{Missing, Int64}} (alias for Array{Union{Missing, Int64}, 1})

julia> if typeof(myvar) == Vector{Union{Int64,Missing}} println("yes"); end
yes

Although another way might be to use type dispatch so that the plot functions have multiple methods which handle the missing cases.

However, I think it's maybe best to simply document that vals should be a vector which assigns a value to each polygon and let people construct such a vector on their own. Perhaps some utility functions are worthwhile, but the beauty of Julia is that we have available generic machinery to work together with our specialized machinery. DataFrames in particular give us lots of data manipulation techniques to use.

Another option might be to allow to put a Dictionary which maps from identifier to value rather than a vector. The plot function could take this dictionary, and a function that takes a polygon and returns a key. This would be efficient when there are MANY polygons which all have the same identifier (for example if there is a particular "county" in Indonesia or something which consists of 30 small islands). Then inside the plotting function to construct the vector you iterate through the polygons, extract the key (which might be a tuple of values) and then lookup the key in the Dictionary to extract the numeric value, something like:

function plot(polygons, keyfun, valdict,...)
   valvec = [valdict[keyfun(polygon)] for polygon in polygons] 
   plot(polygons,valvec,...)
end

[joa-quim]

Thanks, the straight test isa(myvar, Vector{Union{Int64,Missing}}) I knew but then there are the other possible numeric types. Maybe isa(myvar, Vector{Union{Real,Missing}}) works.
We had an attempt to relax some types to facilitate multidispatch but it collided with the AbstractArray interface (see).

The idea of letting the level variable be a dictionary is interesting.

Also realized that the level option is not documented in plot. Quite likely because this was a latter addition to GMT itself and I forgot to update. Must fix that.

[dlakelan]

julia> isa(myvar,Vector{Union{Real,Missing}})
false

julia> typeof(myvar)
Vector{Union{Missing, Int64}} (alias for Array{Union{Missing, Int64}, 1})

Looks like no, seems like you'd have to do something more complicated.

[dlakelan]

Regarding this, the strategy should be to replace missing by NaN and not deleting the polygons. With NaNs the polygons are painted with the NaN color of the colormap (cpt) that is white by default but configurable.

Ok, I'm modifying my example, this is actually really easy, just

joineddata.POPESTIMATE2020 = replace(joineddata.POPESTIMATE2020,missing => NaN)

[dlakelan]

For including output into a weave / jmd document, i'm just using plot with figname="..." and then including the image into the document.

[joa-quim]

Will have to confirm with weave but in jupyter and pluto one only have to do plot(...) and no figname. The image is automatically inserted in the notebook.

[dlakelan]

when I just do plot as the last command in a code block it doesn't include the plot, so for now I'm using figname and explicitly including the file. I'm happy to look into what needs to happen to make that work, we should maybe ask on the Julia Discourse, or on the Weave repo discussion or something.

[joa-quim]

Want to finish some other stuff before and then will make a try with weaver. Probably we'll have to ask for some help. I had to do it for Pluto because it was not obvious at all what had to be done but maybe now it will be an adjustment of previous recipes.

[dlakelan]

Ok, I have an interesting additional plot type. We can calculate the whole of the population estimate of the US... what is the smallest area that contains 50% of the population? We can plot this as 0,1 variables with a :plasma based cpt.

So the question is how to calculate this set of counties... I'm thinking I calculate density, and then order by decreasing density, and then cumsum the populations, and include everything with cumsum less than 50% of total... Will give it a try!

[dlakelan]

Ok, solved... I pushed the code, and here's the 50% high density subset:

[dlakelan]

I suspect a subtle bug, when counties consist of more than one polygon. I'll fix it later.

[dlakelan]

I believe I went through and fixed a number of problems in processing the data, including removing the "000" county (ie. the whole state), and calculating the area of counties that have more than one polygon correctly (by grouping data frames and adding up the areas for given county). The current version should produce valid maps I think.

I think this shows choropleth maps pretty well so far. What do you think should be the second kind of map we look at? Maybe something that colors things continuously based on a field like temperature?

[joa-quim]

Hi,
I've made a couple of changes in master that now set under the hood the close=true,fill="+z" options so it means they can be dropped when register a new version. Also dropped the attributes from the header field of GMTdataset and learned how to recognize the Vector{Union{Missing, Real}} so now you can pass one of such vectors in level option and the missing get replaced by NaNs automatically.

Regarding an example with temperatures I know where o get temperatures for oceans since they are easily measured by satellites but the same is no true for atmosphere. I think the most reasonable source for this is from meteorological models but I have no experience to get that type of data. Anyway, I think you are now aiming for an example that uses grid or image data that I mentioned in one of my first posts.

You may want to have a look at RemoteS.jl Gallery in the docs for a some examples on using Satellite data.

[joa-quim]

Hi, some more changes.
I've just tagged a new version where all singletons of Vector{GMTdataset} are converted to GMTdataset before returning them. That is, al length-1 vector of GMTdataset are converted to GMTdataset . This will break slightly your code where you compute the population density.

[dlakelan]

Hi @joa-quim it's been a complicated couple weeks, I will look into some of these things next week after other stuff settles down. thanks for continuing to work on this project! Yes I think next we try some temperature or similar data. I'll see what I can find.

[joa-quim]

That's fine, no worries. Regarding the temperatures, as I told you sea surface temperatures are easy to find. Atmospheric I think we can't escape from model data.
Have a look too a https://www.generic-mapping-tools.org/RemoteS.jl/dev/gallery/Aqua_sst/remotes_L2_SST/