0

The polygons are plotted using the US census county file:

https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip

Read into R with :

UScounties <- rgdal::readOGR(dsn = paste0(getwd(),"/mapdata/Counties"),
                      layer = "cb_2018_us_county_500k")

> head(UScounties@data)
  STATEFP COUNTYFP COUNTYNS       AFFGEOID  **GEOID**    NAME LSAD      ALAND   AWATER
0      21      007 00516850 0500000US21007 21007 Ballard   06  639387454 69473325
1      21      017 00516855 0500000US21017 21017 Bourbon   06  750439351  4829777
2      21      031 00516862 0500000US21031 21031  Butler   06 1103571974 13943044
3      21      065 00516879 0500000US21065 21065  Estill   06  655509930  6516335
4      21      069 00516881 0500000US21069 21069 Fleming   06  902727151  7182793
5      21      093 00516893 0500000US21093 21093  Hardin   06 1614569777 17463238

spTransform...

UScounties <- spTransform(UScounties,CRS("+init=epsg:4326"))

And then the correct map is drawn when plotting with leaflet.

Now, I want to join other county data by GEOID, e.g., Unemployment data from USDA.gov:

https://www.ers.usda.gov/webdocs/DataFiles/48747/Unemployment.xlsx?v=307.4

UScounties@data <- merge(UScounties@data,Unemployment, by = c("GEOID"= "FIPS_Code") )

Check if the merge is correct:

> UScounties@data[UScounties@data$FIPS == 36013,c(1:7,11:12)]

      FIPS STATEFP COUNTYFP COUNTYNS       AFFGEOID GEOID       NAME State             Area_name
1834 36013      36      013 00974105 0500000US36013 36013 Chautauqua    NY Chautauqua County, NY

Use Leaflet to plot counties with popups:

UScounties@data$AWATER <- as.numeric(UScounties@data$AWATER) 

pal <- colorNumeric(
  palette = "Blues",
  domain = UScounties$AWATER)

leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addPolygons(data = UScounties,
              color = "white",
              fillColor = pal(UScounties$AWATER),
              weight = 1,
              smoothFactor = 0.2, 
              opacity = .77,
              fillOpacity = 0.9,
              highlightOptions = highlightOptions(color = "white",
                                                  weight = 2,
                                                  bringToFront = TRUE),
              labelOptions = labelOptions(style = list("font-weight" = "normal", padding= "3px 8px"),
                                          textsize = "15 px",
                                          direction = "auto"),
              popup = paste("County: ", UScounties@data$NAME, "<br>",
                            "Land: ", UScounties$ALAND, "<br>",
                            "Water: ", UScounties$AWATER, "<br>")) %>%
  addLegend("bottomleft",
            title = "Legend",
            pal = pal,
            values = UScounties$AWATER,
            opacity = .4)

Check map for correct plot:

Popup for Manhattan

Data for this polygon is displayed incorrectly.

Anyone come across this before?

1
  • The problem with this is that there's no Wilson County hereabouts - this is the New York area.
    – Spacedman
    Aug 15 at 10:31
1

I can't get that merge to work as you've written it. It looks wrong.

Without knowing exactly what you did, with code we can re-run, I don't think we can totally get to the bottom of this. But I suspect what's happened is that you've replaced the @data element of the spatial polygons data frame with one that is the wrong size because your match didn't match for everything. If you did a "left join" then non-matching rows will mean fewer rows in the output table, and if you change the @data part of a spatial polygons data frame without also adjusting the @polygons part you will break the link between geometry and attributes, which is what this looks like (now that I know there's no "Wilson County" in New York).

I'd do it like this. First read the data:

UScounties <- rgdal::readOGR(dsn = ".", layer = "cb_2018_us_county_500k")
UScounties <- spTransform(UScounties,CRS("+init=epsg:4326"))

Then read the unemployment table. I'm reading this directly from the XLSX using the openxlsx package:

Unemployment = openxlsx::read.xlsx("./Unemployment.xlsx",startRow=5)

Then do the merge by matching the GEOID from the first item with FIPS_Code on the second item:

USMerge = merge(UScounties, Unemployment, by.x="GEOID", by.y="FIPS_Code")

Note this keeps rows in UScounties that don't match identifiers. So the number of rows is the same:

> nrow(USMerge) == nrow(UScounties)
[1] TRUE

but the missing ones have NA values, eg:

> which(is.na(USMerge$Metro_2013))
 [1]  327  328  502  503  504  505  506 1068 1069 1210 1475 1476 1477 2725

If you are expecting all items to match, you can test this beforehand with:

any(is.na(match(UScounties$GEOID, Unemployment$FIPS_Code)))

(read as "are there any missing matches for the GEOID in the FIPS_code table?").

Plotting this now with leaflet shows this is New York county...

enter image description here

Note 1: First I saved the XLSX to a CSV and read it, but that meant lots of leading zeroes in codes like "06032" were dropped, leading to more mismatches. Instead of trying to fix that, I discovered reading with openxlsx::read.xlsx preserved leading zeroes.

Note 2: this kind of problem is less frequent if you use the sf spatial data classes instead of the sp ones. Amongst other advantages its also much much harder to break the polygons-attributes link since the polygons are just another column in the data frame.

6
  • Ahh interesting about the merge(). I was previously using dplyr::left_Join, and switched to merge as I was troubleshooting (tried both- same wrong result). I still don't get any errors when I enter merge as is. I just checked the merge library documentation and the call as I have it is ok. I tried merging directly to the spatialpolygon layer and got the same incorrect results.
    – T2029
    Aug 13 at 14:12
  • I'll have another go with the real data in a bit.
    – Spacedman
    Aug 13 at 17:54
  • I cant see anything in help(merge) that implies by=c("a"="b") is a valid specification. You can do by=c("a","b") to merge on a and b in both data frames but that's not what c("a"="b") does in dplyr_*join functions.
    – Spacedman
    Aug 13 at 18:12
  • The merge doesn't fail. I showed the code and output proving the data is correctly linked. The polygon has the wrong data.
    – T2029
    Aug 13 at 21:16
  • The merge line as you wrote it cannot work. Its not valid syntax. You must have done it some other way, in which case you should edit your question to show how you did it so we can replicate the error you see in the data. It also means I can't run your "check if the merge is correct" line because when I do fix the merge I don't get a FIPS column.
    – Spacedman
    Aug 14 at 8:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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