Stata update

[ericmelse]

Dear Ben, not an issue but a note (in case that you are not aware) that the current Stata update ( "04 Oct 2023") includes:

graph twoway now supports aspect ratios that provide common units in the x and y dimensions. Adding option aspectratio(1, units) creates a graph where one unit in the x dimension and one unit in the y dimension each take up the same distance on the plot. This is useful for plotting things such as latitude against longitude, which have common scaling of their units.

I suppose that this is relevant to geoplot.

Best, Eric

[benjann]

Oh wonderful, thanks for letting me know. I was asking for such a feature at the last UK Stata meeting, but I was not aware, that it has now already been introduced.

[ericmelse]

Yes, I thought so! By the way, I want to create geoplots for The Netherlands but I only need the boundary lines of the country as such and the (12) regions (provinces). So, I do not need rivers etc. Would you know where to get the most appropriate file for that?

[benjann]

I don't know. Eurostat? Statistics Netherlands?

[benjann]

Note that since the update from 09oct2023 geoplot now employs aspectratio(1, units) if used in Stata 18.

[ericmelse]

Dear Ben, after a substantial effort to get my hands on a set of shapefiles (shp & dbf) for The Netherlands, I finally succeeded. For this I had to resort to using a R package tmap that was created by Martijn Tennekes who is affiliated with Statistics Netherlands (CBS, see his Github here). Actually, for that package he worked with the author of a book on using R to create maps, see here for the second edition (which is work in progress, see his Github here). The little trick that I used was to export the shapefiles for The Netherlands from this R package, here I learned here about this trick (the map there is not as good as there are errors in it), using in R the following:

library(tmap)
data(NLD_prov)
new_file <- glue("NLD_prov.shp")
st_write(NLD_prov, new_file, append = FALSE)

which will put the set of shape related files in the current working folder of R. Using these I am now able to use them in Stata with your fine package, like:

* Setup
geoframe create regions NL_Data.dta, id(_ID) coord(xcoord ycoord) shpfile(NLD_prov_shp.dta)

* Create map
geoplot (area regions popultn, color(vlag, intensity(.5)) levels(17) ) ///
    (line regions, lw(*.3) lc(gs7)) , xsize(4.2) ysize(5) aspectratio(1, units) ///
    clegend(pos(5) ) /// clegend_suboption(color(gs5))) /// cuts(label) lc(gs5)) ///
    ztitle("Population", orientation(rvertical) s(*.6) c(gs5) m(l-0 r+0)) ///
    zscale(lc(gs8)) zlabel(, labc(gs5) labsize(*.8)) /// 
    ztick(, tlcolor(blue) tlwidth(*.5) tlength(*.5)) /// this is not used by geoplot but accepted
    compass(t(2) pos(2) size(3) xm(.2) ym(.2) c(gs9)) /// 
    sbar(units(km) ym(2) c(gs8) lab(size(*.5) c(gs5)))

which results in: This is a fine result that certainly meets my data visualization needs.

ztick issue But note that the line of code controlling the color, width and length of the clegend axis is accepted by geoplotbut not processed: ztick(, tlcolor(blue) tlwidth(*.5) tlength(*.5)) I have tried many alternatives but, reading the help file instructions, I assume that it should work. Which might possibly mean that there is something amiss in the ado file.

compass labels issue Another rather minor issue is the current control of the compass size. I used here type(2) but my remark applies to each type. What happens when the size is changed from the default 5, like decreasing or increasing it, that then the size of the compass labels are not changing accordingly. As a result, in my example above, the point of the 'pointer' is moving inward the capital letter N (for North). Looking for perfection, as always, I suggest to include an option to control the size as well as the font type used for compass labels, like: compass(t(2) pos(2) size(3) fontsize(medium) font("Merriweather")) Or something similar. Time permitting, maybe you can implement this.

I attach the files should you want to replicate. Geoplot_NL.zip

[benjann]

compass labels issue: see suboption label() in compass(). For example, you could type compass(label(size(medium))); don't know how to change font though. Mayne I will add an option that allows you to define custom labels (so that the font can be specified). Possibly I will then change syntax to

label(["text" ...][, suboptions])

ztick issue: the option is processed, but ztick() is only for ticks without label, not for labeled ticks. Since no unlabeled ticks are printed, the option has no effect. What you probably want is zlabel(, tlcolor(blue) tlwidth(*.5) tlength(*.5)).

[ericmelse]

Thank you Ben, for your helping hand. I have continued with adjusting and cleaning my code as well as somewhat improving it. I did the following:

• I removed the part to draw line regions because the area regions have the option to be drawn with a boundary line. This means 'less computing' which reduces the time to draw the plot. But, more importantly, I observe a more fine drawing of the boundary as such. So the plot of the map has improved.
• Instead of using levels(17) to control the color coding, I used the recommended alternative option cuts(#min(#step)#max) that requires some upfront calculations to get the parameters right.
• Using cuts provides another benefit: it is more easy to also use the option zmlabel but without labels (using nolab). Then I adjust their length and position: tl(*4.35) tp(i)) so that these are drawn inside the scale, as to separate the color coded steps. I do this because the colors are used as such in the map within the boundary of a region. Also, this is helpful to mitigate the visual effect of simultaneous contrast, typically close to the border where two color steps meet.
• Regarding the font that I now use (Sabon Small Caps ), I used graph set window fontface (before using geoplot) that applies to all text and labels. Because the legend labels are populated with many digits (now formatted) using Sabon Small Caps is really effective and aesthetically pleasing. Also the compass labels now have a nice 'classic look'. But, their distance to the compass is still a bit too close (in my opinion). I see this with all the fonts that I used. Maybe you can increase the distance a little bit between the compass and its labels.
• I also use colors for the scale that meet the requirement of proper color perception by people who suffer from CVD as per Jann, B. (2023).

So, my revised code is:

local textclr ""58 144 254"*1.4"
geoplot (area regions popultn, color("58 144 254*.6" gs14 "168 144 8", ipolate(100, positions(0 .35 1)) ///
        intensity(.5)) lw(*.3) lc(gs7) cuts(350000(196875)3696875) ) ///
    , xsize(4.2) ysize(5) aspectratio(1, units) graphreg(fc(white) lc(white) ilc(white)) ///
    clegend(pos(5)) ztitle("Population", orientation(rvertical) s(*.7) c(`textclr')) ///
    zscale(lc(gs8) lw(*.6)) ///
    zlabel(, format(%9,0fc) labc(gs5) labsi(*.9) tlc(gs5) tlw(*.3) tl(*.5)) ///
    zmlabel(350000(196875)3696875, nolab tlc(gs5) tlw(*.3) tl(*4.35) tp(i)) ///
    compass(t(1) pos(2) si(2.3) xm(.2) ym(.2) c(gs9) lab(si(*.5) c(`textclr'))) /// 
    sbar(units(km) ym(2) c(gs8) lab(si(*.55) c(gs5)))

which results in:

[benjann]

Hi Eric,

some comments:

• whether you use line regions or area regions with lines turned on doesn't really matter; the result will look the same as long as you use the same line options (i.e. your comment "I observe a more fine drawing of the boundary as such" indicates that you used different line options)
• ipolate(100) generates 100 colors, but you only define 17 bins; so I would type ipolate(17) (it doesn't really matter much because geoplot will automatically reduce the 100 colors to 17 by interpolation, but why not give it the right number of colors in the first place)
• aspectratio() is not needed; geoplot takes care of that automatically (it has always done so, even before the units suboption became available)
• graphreg(fc(white) lc(white) ilc(white)) is not needed; geoplot sets the style of the graph region to none internally
• you can use zmtick() rather than zmlabel(, nolab)
• use compass(label(gap())) to set gap between the compass and its labels

Here's my reproduction of your graph (ignoring the font settings):

geoframe create regions NLD_prov
local textclr `""58 144 254"*1.4"'
geoplot ///
    (area regions popultn, lw(*.3) lc(gs7) ///
        cuts(350000(196875)3696875) /// => 17 bins
        color("58 144 254*.6" gs14 "168 144 8", ///
            ipolate(17, positions(0 .35 1)) intensity(.5))) ///
    , xsize(4.2) ysize(5) ///
      clegend(pos(5)) ///
      ztitle("Population", orientation(rvertical) s(*.7) c(`textclr')) ///
      zscale(lc(gs8) lw(*.6)) ///
      zlabel(, format(%9,0fc) labc(gs5) labsi(*.9) tlc(gs5) tlw(*.3) tl(*.5)) ///
      zmtick(350000(196875)3696875, tlc(gs5) tlw(*.3) tl(*4.35) tp(i)) ///
      compass(t(1) pos(2) si(2.3) xm(.2) ym(.2) c(gs9) ///
          lab(si(*.5) c(`textclr') gap(50))) /// 
      sbar(units(km) ym(2) c(gs8) lab(si(*.55) c(gs5)))

ben

[benjann]

By the way, I think you should apply

geoframe generate plevel

after loading the data because there seem to be some enclaves (e.g. in the Tilburg region). This makes sure that the enclaves will be treated as such. That is, type

geoframe create regions NLD_prov
geoframe generate plevel
local textclr `""58 144 254"*1.4"'
geoplot ...

ben

[ericmelse]

Right, Ben, many thanks for you suggestions to clean up my code. And, yes, indeed I also noticed the enclaves and worried what to do about them. Beside the Belgium enclaves, near Tilburg, there are 30 of them, see here), another 'enclosure' that got colored like the regional area of Zeeland in my previous examples is the so-called Veerse Meer. So, indeed, your suggestion to use:

geoframe create regions NLD_prov
geoframe generate plevel

worked perfectly (21 nested polygons were found): My next step will be if I can add rivers to this map (I am looking for a shape file of those).

[benjann]

re compass labels: With the latest update it is now possible to define custom labels. The syntax is

compass(label(text(<labels>)))

Where <labels> are the labels to be used for north, south, east, and west. For example, you could type

compass(label(text("{bf:N}" "S" "E" "W")))

to print N in boldface. The different compass types have different default settings, but all of them now support labels in the same way. For example, the following would label north and south in compass type 3:

compass(type(3) label(text("{bf:N}" "S")))

ben