R Leaflet SpatialPolygonsDataFrame @data join w/ dataframe; data on wrong polygon

Question

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:

Data for this polygon is displayed incorrectly.

Anyone come across this before?

Answer 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...

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.