IN WHICH TEXAS COUNTIES ARE WOMEN MORE AFFECTED BY GONORRHEA? PART II
To accomplish the first task of creating boundaries around Public Health Regions (PHRs), I first located a list of all counties within each public health region. I read through the given table and noted down a list for my reference. Then, I dove right in!
- I introduced a new numeric variable ‘phr‘ to the original ‘texas_state’ map data table using a series of IF–THEN; ELSE IF–THEN; ELSE; statements that would read the mapsgfk.us_counties variable ‘County’ and assign each value (a single row with distinct LAT/LONG) a phr between 1-11 based on the information within the linked list. I ran a PROC CONTENTSand PROC SORT on data = texas_state; by phr; to prep for the next step: the removal of county borders, leaving me with a map of outlined Public Health Region borders.
- Here I also added a PROC GPROJECT step for my input map ‘texas_state’ in order to specify the type of projection (Albers) and set parameters via PARMOUT into work.projparm (project parameters). This will allow smooth implementation of additional annotations in the near future (See Step X)
- Because I needed to retain the borders around counties for the sake of my ultimate goal of visualizing counties where a particular sex experienced higher incidence rates of Gonorrhea than the other, I had to create a new map data table using PROC GREMOVE that I called ‘texas_phrmap’. I learned that the only way to “overlay” the two maps is to implement the PHR border map into the PROC GMAP procedure as an annotation (shout out to Mr. Eberhart for his incredibly helpful guide to the Annotate facility).
- Here is an excerpt from page 10 of his guide: “To create annotate dataset to draw polygon boundaries, first remove observations from the traditional map dataset that do not refer to the polygon boundaries using PROC GREMOVE. The GREMOVE procedure combines unit areas defined in a map data set into larger unit areas by removing shared borders between the original unit areas. Then use the data step to create an annotate dataset that has functions to start polygons and continue polygons as needed.“
- Using the example found in this guide, I was able to create an annotation I named ‘anno_phr’ that created a black polygon border around each phr segment from my texas_phrmap. I updated my PROC GMAP statement to coutline=beige and cempty=beige which would help distinguish between county borders and PHR boundaries.
- SAS GMAP is limited in how much you can increase the line size in the default resolution. In order to make the PHR boundaries much thicker compared to county borders, I used GOPTION XPIXELS = 4000 YPIXELS = 3500 to artificially increase the scale of the image. By doing this, I could increase the size of the PHR boundary annotation function.
- Note: if you opt for this method, you will also have to scale up the other components of the image accordingly, for example the font size of the title.
Now I was ready to tackle my next challenge: adding labels at major cities and the Public Health Regions! I figured I would need the LAT and LONG or X and Y coordinates for each of the cities I wished to have on my map. It was at this point I realized I was in over my head! I could not find any online resources similar to Mr. Eberhart’s guide for this task. I couldn’t find any archived SAS learner’s query similar enough to mine on StackExchange. This was when I stumbled upon The Graph Guy aka ex-SAS employee Robert Allison’s blog. Following along a few of his ‘How-to’ style data visualization posts, I finally encountered the instruction I so desperately needed. With newfound understanding of parameters and map data, I created a new output table:
data city_names; set mapsgfk.uscity (where=(statecode='TX'
and city in ('Austin' 'Dallas' 'Houston' 'Tyler' 'Temple'
'El Paso' 'Lubbock' 'Amarillo' 'Laredo' 'San Antonio'
'Harlingen' 'Odessa' 'Abilene' 'San Angelo' 'Arlington' )));
run;
proc gproject data=city_names out=anno_cities latlong eastlong degrees dupok
parmin=projparm parmentry=texas_state;
id ;
run;
I then created my second annotation dataset called ‘anno_cities’ from the gproject output that would place a large black dot [function=‘pie’ ; style =‘psolid’ ; rotate =360;] and the city name [function= ‘label’; color=‘black’; text=city; style=”Albany AMT/Bold”;] at the LAT, LONG values of each of the specified cities. I used IF–THEN statements here to set desired positions for the text about each dot. Here is a great guide to the POSITION statement.
Note: I fiddled around with color patterns numerous times throughout my graph making journey, so please don’t be alarmed by the change in appearance between the previous step and this upcoming one. Here is how I updated the color palette:

pattern1 v=s c=cx4dac26;
pattern2 v=s c=cxb8e186;
pattern3 v=s c=cxfffff0;
pattern4 v=s c=cxf1b6da;
pattern5 v=s c=cxd01c8b;
to

pattern1 v=s c=cx5aaefd;
pattern2 v=s c=cxa5d4fd;
pattern3 v=s c=cxf4eeff;
pattern4 v=s c=cxffb0ce;
pattern5 v=s c=cxf67eaa;
With the effects of PROC GPROJECT on the ‘texas_state’ map data and the addition of a second annotation ‘anno_cities’ , which I combined with the first annotation like so:
data anno_all; set anno_phr anno_cities; run;
My map now looks like this:

ADDING PUBLIC HEALTH REGION LABELS
I used Google Maps to right-click and identify LAT/LONG coordinates where PHR labels could be placed (low interference with city labels and PHR boundaries). I entered these values into a table called I ‘phr_numbers’ as datalines. I then ran a PROC GPROJECT for phr_numbers and formatted the output data anno_numbers as desired. Just as the previous step, I introduced the new annotation ‘anno_numbers’ into ‘anno_all’. As a bonus, I wanted to clean up the appearance of the map by adding white circles behind each PHR number to increase contrast thus visibility of the label. I created an fourth and final annotation called ‘phr_dots’ for which the input is identical to the phr_numbers. This ensures they are projected at the same coordinates within the map.


data anno_numbers; set anno_numbers;
xsys='2'; ysys='2'; hsys='3'; when='a';
length function $8 color $12 style $35 text $100;
function='label'; position='+';
style='albany amt/bold';
rotate=.; size=2; color='black'; text=trim(left(region));
run;
...
data anno_dots; set anno_dots;
xsys='2'; ysys='2'; hsys='3'; when='a';
function='pie'; style='Psolid'; rotate=360; size=1.5; color='white'; output;
function='pie'; style='Pempty'; rotate=360; size=1.5; width=2; color='black';
output;
run;
...
data anno_all; set anno_phr anno_cities anno_dots anno_numbers ;
run;
All done! Isn’t it so exciting!?
My other choropleth map projects were made very similarly. The primary difference was in the original data introduced to the map creation program as datalines (for incidence rate over time, 2018-2019, total_rate was calculated for each STD for an individual year and then each output table was converted into tab delimited file that was used to create a year-specific map)