# Creating Maps from Shape Files for “Mixed” Labeled Data

Following suggestions from an extremely useful and informative solution to a previous question regarding shape files (Extracting a single layer from a SHP file) by ssch, I have been able to create two maps of Washington State Counties and its precincts.

The map of the precinct data for the entire state looks like this:

The map for the county data for the entire state looks like this:

Data can be found at WA Secretary of State Website. The shape file in question can be downloaded at Statewide_Splits_2020General.zip. The file used to generate the graphics is given in the code below.

However, because there are a great many precincts and these are soon to change with redistricting, I would like to:

1. be able to create precinct maps separately for each county, rather than for the entire state. Because the scaling of the map for the entire state creates minute maps for areas with high population density, they are not easily distinguished.

2. create GeoRegionValuePlots for each county/precinct for different measured properties (ie. population, demographic composition, economic data, election results, etc.)

Data that I would like to see displayed graphically are the 2020 WA general election results for each precinct. This data can be downloaded at: Excel file for All Counties by Precinct.

Following, but slightly modifying ssch's code creates a nice map of precincts for the entire state.

 data =   First@Import["Statewide_Splits_2020General.shp", "Data"];
data /. Rule[a_, b_] :> a
labels = "PrecName" /. ("LabeledData" /. data);
geometry = "Geometry" /. data;
labeled = {labels, geometry}\[Transpose];
gathered = GatherBy[labeled, First@# &];
gathered // Length
Union@labels;
g =  Graphics[MapIndexed[{ColorData[1] @@ #2,
Tooltip[#[[All, 2]], #[[1, 1]]]} &, gathered],
ImageSize -> {1000, 1000}]


substituing "CountyName" for "PrecName" in defining the labels produces a nice map of the counties.

However, I haven't been able to figure out how to create maps for each individual county showing the precincts within them that can then be used by GeoRegionValuePlot. Various attempts with joins and selections when creating the variables labels and labeled create problems in transposition that I haven't been able to overcome.

1. How can this be done?

2. How can external data for each precinct be indexed to the precinct/county shape files and to a ColorFunction so that a GeoRegionValuePlot can be associated with the appropriate polygon characterizing its shape?

• Could you link directly to the file you're talking about and perhaps include one or two images to illustrate what you're looking for? – Carl Lange Jan 3 at 14:39
• Now edited in. Thanks for tie interest in helping. What I am looking for is to be able to display maps of each county and their precincts and then ultimately using GeoRegionValuePlots displaying various results for each county and it's included precincts using a color function, one in which the percentage of Democratic voters is shown in shades of blue and Republican voters in shades of red. This exercise is primarily to help me understand how to manipulate shape files to extract subsections of the full file. – Stuart Poss Jan 3 at 18:20
• It should be noted that election results are not available for all precincts, just for those in "participating" counties. – Stuart Poss Jan 3 at 19:07
• Would you mind joining me in this chat room when you have time to help me try and nail down what you need? chat.stackexchange.com/rooms/117993/… I'm not from the US so I may need clarification on some of the terms you're using. – Carl Lange Jan 3 at 21:52
• The algorithms provided by Carl and Jean-Pierre both solve the problem posed in part 1 of my question, with part 2 perhaps best asked as a separate question. They are difficult to compare in that the first creates a Dataset and generates code for several county maps without reference to the name of a specific file and less control for image size important for counties with more precincts, but it offers the advantage of quickly organizing maps for each county. Accounting for size, (Windows 10, 64 bit machine) first requires about 10.9 sec and the second about 10.35 sec. Thanks to you both. – Stuart Poss Jan 8 at 21:32

Now that I can access the data file, I can provide a more complete answer. First, here are a few functions:

(* FUNCTION: input a shape file and returns a data file *)

shapeFileToData[path_] := (
First@Import[path, "Data"]
)

(* FUNCTION: Extract the top rules left side *)

getRulesLeft[data_] := (data /. Rule[a_, b_] :> a)

(* FUNCTION: Extract data corresponding to a specific rule *)

extractRuleData[data_, rule_] := (
data /. {Rule[Except[rule], _] :> Nothing, Rule[a_, b_] :> b}
)

(* FUNCTION: Given LabeledData, a rule name and its value as a \
filter, get the values matching another rule *)
labeledDataForRuleValue[data_, ruleName_, ruleValue_,
otherRuleName_] := (
rd = extractRuleData[data, ruleName];
pos = Flatten@Position[Flatten[rd], ruleValue];
other = Flatten@extractRuleData[data, otherRuleName];
Map[other[[#]] &, pos]
)

(* FUNCTION: starting with the full data, returns the geometry \
associated with a LabeledData rule name and value *)

geometryDataForRuleValue[fullData_, ruleName_, ruleValue_] := (
ldata = extractRuleData[fullData, "LabeledData"];
rd = extractRuleData[ldata, ruleName];
pos = Flatten@Position[Flatten[rd], ruleValue];
geometry = First@extractRuleData[fullData, "Geometry"];
Map[geometry[[#]] &, pos]
)


Now let us use these functions:

fullData =
shapeFileToData[

topRules = getRulesLeft[fullData]
(*{"LayerName", "Geometry", "Labels", "LabeledData"}*)

labels = extractRuleData[fullData, "Labels"]
(*{{"PPartName", "PrecCode", "FIPS", "CountyName", "PrecName"}}*)

labeled = extractRuleData[fullData, "LabeledData"];
counties = extractRuleData[labeled, "CountyName"]

(* Precints in Adams County *)
precincts =
(* {"RITZVILLE RURAL SE", "RITZVILLE RURAL SE", "HATTON RURAL", \
"FAIRVIEW", "RALSTON", "RITZVILLE RURAL SE", "OTHELLO RURAL #5", \
"HATTON RURAL", "HATTON RURAL", "WASHTUCNA CITY #1", "OTHELLO RURAL \
#4", "RITZVILLE WARD #3", "WASHTUCNA RURAL", "OTHELLO RURAL #6", \
"RALSTON", "BENGE", "RITZVILLE RURAL NW", "OTHELLO RURAL #2", \
"FAIRVIEW", "LIND RURAL", "CUNNINGHAM", "RITZVILLE WARD #2", \
"RALSTON", "RITZVILLE RURAL NW", "RITZVILLE RURAL SE", "SCHRAG", \
"HATTON CITY", "OTHELLO RURAL #1", "LIND RURAL", "WASHTUCNA RURAL", \
"BENGE", "OTHELLO CITY #5", "RITZVILLE RURAL SE", "OTHELLO CITY #4", \
"OTHELLO CITY #3", "CUNNINGHAM", "LIND RURAL", "OTHELLO RURAL #6", \
"WASHTUCNA RURAL", "OTHELLO CITY #1", "OTHELLO CITY #2", "FAIRVIEW", \
"SCHRAG", "SCHRAG", "RITZVILLE WARD #1", "LIND #1", "OTHELLO RURAL \
#3", "RITZVILLE WARD #5", "WASHTUCNA RURAL", "BENGE", "OTHELLO CITY \
#6", "OTHELLO RURAL #6", "CUNNINGHAM", "RITZVILLE WARD #4", \
"RITZVILLE RURAL SE", "RITZVILLE RURAL NW", "RITZVILLE RURAL NW", \
"BATUM", "RITZVILLE RURAL SE", "BENGE", "BENGE", "OTHELLO RURAL #1", \
"LIND RURAL", "BATUM", "BATUM", "BATUM", "BATUM"} *)

(* County where the precints of "OTHELLO CITY #2" is located *)
county = labeledDataForRuleValue[labeled, "PrecName",
"OTHELLO CITY #2", "CountyName"]

(* PPartName associated with Adams County *)
partName =
(*{"0122.1", "0122.3", "0324.3", "0322.1", "0222.2", "0122.2", \
"0515.1", "0324.2", "0324.1", "0212.1", "0514.1", "0113.1", "0221.4", \
"0516.1", "0222.3", "0211.3", "0121.3", "0512.1", "0322.2", "0321.1", \
"0325.3", "0112.1", "0222.1", "0121.1", "0122.5", "0323.1", "0313.1", \
"0511.1", "0321.2", "0221.3", "0211.4", "0415.1", "0122.6", "0414.1", \
"0413.1", "0325.1", "0321.3", "0516.3", "0221.2", "0411.1", "0412.1", \
"0322.3", "0323.3", "0323.2", "0111.1", "0311.1", "0513.1", "0115.1", \
"0221.1", "0211.5", "0416.1", "0516.2", "0325.2", "0114.1", "0122.4", \
"0121.999", "0121.2", "0123.2", "0122.7", "0211.1", "0211.2", \
"0511.2", "0321.999", "0123.4", "0123.3", "0123.5", "0123.1"}*)


Note that precincts name appear more than once. These entries are differentiated by the PPartName. We will now get the geometries associated with Adams county. The map will show different colors for different PPartName. The tool tip shows the name of the precincts followed by the PPartName in bracket.

geodata = geometryDataForRuleValue[fullData, "CountyName", "Adams"];
i = 1;
GeoGraphics[{Red,
EdgeForm[{Black}], {RandomColor[],
Tooltip[#, precincts[[i]] <> " (" <> partName[[i++]] <> ")"]} & /@
geodata}, GeoBackground -> "VectorMinimal",
Frame -> {{None, None}, {None, True}},
FrameLabel -> {{None, None}, {None, Style["Adams", Red, Bold, 20]}},
GeoScaleBar -> Placed["Metric", {Left, Top}]]


The following uses a function to keep the same color for all PPartName associated with a specific precinct.

namecolor = <||>;
getColor[name_] := (
If[
MissingQ[namecolor[[name]]],
AssociateTo[namecolor, name -> RandomColor[]];
Key[name][namecolor], Key[name][namecolor]]
)
i = 1;
GeoGraphics[{Red,
EdgeForm[{Black}], {getColor[precincts[[i]]],
Tooltip[#, precincts[[i]] <> " (" <> partName[[i++]] <> ")"]} & /@
geodata}, GeoBackground -> "VectorMinimal",
Frame -> {{None, None}, {None, True}},
FrameLabel -> {{None, None}, {None, Style["Adams", Red, Bold, 20]}},
GeoScaleBar -> Placed["Metric", {Left, Top}]]


UPDATE: Part 2 of the question required coloring the precincts in shades of blue or red acccording to the percentage of vote in the presidential election. We will import the new file as a dataset. We will be using a very similar code as before, with some modifications to the function getColor. When it is time to introduce a color, instead of random color, we will call the function getPoliticalColor. The new file uses a two letter CountyCode instead of CountyName. This code will now be given to getColor with the name of the precinct (in case we have precincts with the same name in two different counties). The example uses the county of Chelan.

(* FUNCTION: given the county code and precinct name, returns a color \
(shades of red or blue) corresponding to the percentage of votes for \
either Trump or Biden *)

getPoliticalColor[precData_, countyCode_, precName_] := (
dem = precData[
Select[(#Candidate ===
"Joseph R. Biden / Kamala D. Harris" && #CountyCode ===
countyCode &&
ToUpperCase[ToString[#PrecinctName]] ===
rep = precData[
Select[(#Candidate ===
"Donald J. Trump / Michael R. Pence" && #CountyCode ===
countyCode &&
ToUpperCase[ToString[#PrecinctName]] ===

ratio =
Which[totalVotes == 0, Return[RGBColor[0.5, 0.5, 0.5]],

finalColor = Which[
Between[{0, 0.1}][ratio], RGBColor[0.8, 0, 0.],
Between[{0.1, 0.2}][ratio], RGBColor[1., 0., 0.],
Between[{0.2, 0.3}][ratio], RGBColor[1.0, 0.1, 0.1],
Between[{0.3, 0.4}][ratio], RGBColor[1.0, 0.3, 0.3],
Between[{0.4, 0.5}][ratio], RGBColor[1, 0.5, 0.5],

Between[{0.5, 0.6}][ratio], RGBColor[0.7, 0.8, 1.],
Between[{0.6, 0.7}][ratio], RGBColor[0.5, 0.67, 1.],
Between[{0.7, 0.8}][ratio], RGBColor[0.2, 0.47, 1.],
Between[{0.8, 0.9}][ratio], RGBColor[0.902, 0.749, 0],
Between[{0.9, 1.0}][ratio], RGBColor[0, 0, 1.]
];

finalColor
)

(* the new dataset *)
precData =

(* the code - note that we are using the symbol "labeled" defined previously with our original data file *)
partName =
labeledDataForRuleValue[labeled, "CountyName", "Chelan",
"PPartName"];
precincts =
labeledDataForRuleValue[labeled, "CountyName", "Chelan", "PrecName"];
geodata = geometryDataForRuleValue[fullData, "CountyName", "Chelan"];
namecolor = <||>;
getColor[countyCode_, precName_] := (
If[MissingQ[namecolor[[precName]]],
AssociateTo[namecolor,
precName -> getPoliticalColor[precData, countyCode, precName]];
Key[precName][namecolor], Key[precName][namecolor]]
)
(* Note below that we are using CH as the two letter code for Chelan *)
i = 1;
GeoGraphics[{GeoStyling[Directive[Opacity[0.8]]],
EdgeForm[{Black}], {getColor["CH", ToString[precincts[[i]]]],
Tooltip[#,
ToString[precincts[[i]]] <> " (" <> partName[[i++]] <> ")"]} & /@
geodata}, GeoBackground -> "StreetMap",
Frame -> {{None, None}, {None, True}},
FrameLabel -> {{None, None}, {None, Style["Chelan", Red, Bold, 20]}},
GeoScaleBar -> Placed["Metric", {Left, Top}]]


• This does work, but all the precincts are the same color and the Tooltips no longer point to the individual precincts, but rather solely to the county name. I hadn't anticipated this outcome, but hoped that the I could get more control over the boundary colors and thickness and not lose the tooltips. I am assuming that all I will need will be the boundaries themselves as I wish to use the maps in GeoRegionValuePlots, which would color the individual precincts based on computed values for external data (ie election results, etc.) – Stuart Poss Jan 4 at 7:59
• Perhaps I am mixing two separate questions, as I see that what you have provided can be used to map g over all counties, with modification of the SomeCountyName labeling and positioning. – Stuart Poss Jan 4 at 8:03
• The link to the data works for me (it's the first entry on that page). – Stuart Poss Jan 4 at 8:06
• I see that the link works now. Previously, I was always redirected to an error page when clicking on it. I'll have another look. – Jean-Pierre Jan 4 at 10:45
• The function replaceUnwantedCharacters[Flatten[rd]] appears to be undefined. Was something left out of the code provided? Perhaps this may explain why the assignment of precincts = labeledDataForRuleValue[labeled, "CountyName","Adams","PrecName"] reutnrs an error Flatten: Nonatomic expression expected at position 1 in Flatten[data]. – Stuart Poss Jan 7 at 4:31

For now, here is another way to extract the data and group it by county. I typically prefer to use Associations rather than lists of rules in Mathematica, mostly for ease of use (I've never really gotten my head around the "Geometry" /. data semantics).

Let's start by reading the shapefile and converting it into an association (in this case, we also want only the First layer, as there is only one).

data = First[
ResourceFunction["ToAssociations"][
{"Data"}]]];


One of the new additions in version 12.2 are GeoPolygons, which are slightly nicer to use (in my opinion) than just Polygon. You don't really need to convert them.

geoms = data["Geometry"] /. Polygon -> GeoPolygon


We also need LabeledData:

ld = d["LabeledData"]


Now we can create a Dataset containing the "de-interleaved" data:

ds = Dataset[Table[<|
"CountyName" -> ld["CountyName"][[n]],
"PrecName" -> ld["PrecName"][[n]],
"FIPS" -> ld["FIPS"][[n]],
"PrecCode" -> ld["PrecCode"][[n]],
"PPartName" -> ld["PPartName"][[n]],
"geom" -> geoms[[n]]|>, {n, Length@geoms}]];


Now we can group by CountyName quite simply:

ds[GroupBy["CountyName"]]


And for example now we can create GeoGraphics for each of the ten smallest counties, showing the precincts:

ds[GroupBy["CountyName"] /* TakeSmallestBy[Length@# &, 3],
GeoGraphics[{EdgeForm[Black], #[[All, "geom"]],
ImageSize -> {300, 300}}] &]


I recognise I may have made things harder by involving the Dataset construct but it's an alternative approach!

• Conceptually, you have made some things simpler. I assume you mean ld=data["LabeledData"], However, all works until I try your last line upon which I get a Part specification geom is not applicable Part error. I have substituted ds[TakeSmallestBy[Length@# &, 5], GeoGraphics[{EdgeForm[Black], #[[All, "geom"]]}, ImageSize -> {300, 300}] &] to get larger county maps and a black edgeform, but either this or your line balks with the Part geom error. Am I missing something here? – Stuart Poss Jan 4 at 7:41
• Yes, my mistake, I have corrected my code! Sorry about that. – Carl Lange Jan 4 at 11:12