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:

Map of WA State Precincts

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

enter image description here

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
 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?

  • $\begingroup$ 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? $\endgroup$
    – Carl Lange
    Commented Jan 3, 2021 at 14:39
  • $\begingroup$ 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. $\endgroup$ Commented Jan 3, 2021 at 18:20
  • $\begingroup$ It should be noted that election results are not available for all precincts, just for those in "participating" counties. $\endgroup$ Commented Jan 3, 2021 at 19:07
  • $\begingroup$ 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. $\endgroup$
    – Carl Lange
    Commented Jan 3, 2021 at 21:52
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Jan 8, 2021 at 21:32

2 Answers 2


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 = 

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

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

labeled = extractRuleData[fullData, "LabeledData"];
counties = extractRuleData[labeled, "CountyName"]
(* {{{Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams,Adams, ...12917... ,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima,Yakima}}}  *)

(* Precints in Adams County *)
precincts = 
 labeledDataForRuleValue[labeled, "CountyName", "Adams", "PrecName"]

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

(* PPartName associated with Adams County *)
partName = 
 labeledDataForRuleValue[labeled, "CountyName", "Adams", "PPartName"]
(*{"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;
  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}]]

enter image description here

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

namecolor = <||>;
getColor[name_] := (
   AssociateTo[namecolor, name -> RandomColor[]]; 
   Key[name][namecolor], Key[name][namecolor]]
i = 1;
  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}]]

enter image description here

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]] === 
         precName) &], {"Votes"}];
  rep = precData[
    Select[(#Candidate === 
         "Donald J. Trump / Michael R. Pence" && #CountyCode === 
         countyCode && 
        ToUpperCase[ToString[#PrecinctName]] === 
         precName) &], {"Votes"}];
  demVotes = First@First@Normal@dem;
  repVotes = First@First@Normal@rep;
  totalVotes = demVotes + repVotes;
  ratio = 
   Which[totalVotes == 0, Return[RGBColor[0.5, 0.5, 0.5]], 
    totalVotes > 0, demVotes/totalVotes];
  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.]

(* the new dataset *)
precData = 
   "Dataset", HeaderLines -> 1];

(* the code - note that we are using the symbol "labeled" defined previously with our original data file *)
partName = 
  labeledDataForRuleValue[labeled, "CountyName", "Chelan", 
precincts = 
  labeledDataForRuleValue[labeled, "CountyName", "Chelan", "PrecName"];
geodata = geometryDataForRuleValue[fullData, "CountyName", "Chelan"];
namecolor = <||>;
getColor[countyCode_, precName_] := (
    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;
  EdgeForm[{Black}], {getColor["CH", ToString[precincts[[i]]]], 
      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}]]

enter image description here

  • $\begingroup$ 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.) $\endgroup$ Commented Jan 4, 2021 at 7:59
  • $\begingroup$ 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. $\endgroup$ Commented Jan 4, 2021 at 8:03
  • $\begingroup$ The link to the data works for me (it's the first entry on that page). $\endgroup$ Commented Jan 4, 2021 at 8:06
  • $\begingroup$ 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. $\endgroup$ Commented Jan 4, 2021 at 10:45
  • 1
    $\begingroup$ 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]. $\endgroup$ Commented Jan 7, 2021 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[
    Import["~/Downloads/statewide_splits_2020general.zip", \

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}]];

enter image description here

Now we can group by CountyName quite simply:


enter image description here

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}}] &]

enter image description here

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

  • $\begingroup$ 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? $\endgroup$ Commented Jan 4, 2021 at 7:41
  • $\begingroup$ Yes, my mistake, I have corrected my code! Sorry about that. $\endgroup$
    – Carl Lange
    Commented Jan 4, 2021 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.