0

In an effort to detect spatial autocorrelation using Moran's I in some variables, I was wondering how I could go about creating a polygonal shapefile from a ggplot-created map. My data lists multiple years (2010-2017) of variables for the 48 contiguous US states to which I've joined to the 'states' shapefile of the'urbnmapr' R package. However, whenever I try to create a shapefile from the joined data, it describes 'map_data' as a 'SpatialPointsDataFrame' as opposed to 'SpatialPolygon'. This gets in the way of the function 'poly2nb' where try to define the neighbor list (Reference I'm using to calculate Moran's I).

I'm just trying to figure out what I'm doing wrong and how I can go about fixing this. Thank you.

> head(df, n = 10)
# A tibble: 10 x 13
# Groups:   YEAR, CENSUS_DIVISION, STATE_FULL, STATE, FIPS, CT_LAT [10]
    YEAR CENSUS_DIVISION STATE_FULL STATE FIPS  CT_LAT CT_LON  GDPr    LI     SD     MM      MA      PA
   <int> <fct>           <fct>      <fct> <chr>  <dbl>  <dbl> <dbl> <dbl>  <dbl>  <dbl>   <dbl>   <dbl>
 1  2010 PACIFIC         California CA    06      37.1  -120. 55147    18  0.838  3.03   1.52    1.77  
 2  2010 PACIFIC         Oregon     OR    41      44.0  -121. 44043     5  0.206  0.219 -0.486   1.28  
 3  2010 PACIFIC         Washington WA    53      47.4  -121. 56541     5 -0.327  0.507 -0.0230  2.42  
 4  2010 MOUNTAIN        Arizona    AZ    04      34.2  -112. 40183     4  0.603  1.32   1.07    0.173 
 5  2010 MOUNTAIN        Colorado   CO    08      39.0  -106. 52452     3 -0.368  0.566  0.115  -0.535 
 6  2010 MOUNTAIN        Idaho      ID    16      44.3  -115. 36865     4 -0.152 -0.568  1.38    0.0780
 7  2010 MOUNTAIN        Montana    MT    30      47.1  -110. 40924     1 -0.823 -1.14   0.256  -0.215 
 8  2010 MOUNTAIN        Nevada     NV    32      39.3  -117. 47670     4  0.842  2.00  -0.118   0.513 
 9  2010 MOUNTAIN        New Mexico NM    35      34.4  -106. 42141     1  0.846  1.76   0.743   0.269 
10  2010 MOUNTAIN        Utah       UT    49      39.3  -112. 44438    -3 -0.571 -0.285  3.23   -0.223 

> states <- data(urbnmapr::states)
> head(states, n = 10)
# A tibble: 10 x 9
    long   lat order hole  piece group state_fips state_abbv state_name
   <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>      <chr>      <chr>     
 1 -88.5  31.9     1 FALSE 1     01.1  01         AL         Alabama   
 2 -88.5  31.9     2 FALSE 1     01.1  01         AL         Alabama   
 3 -88.5  31.9     3 FALSE 1     01.1  01         AL         Alabama   
 4 -88.5  32.0     4 FALSE 1     01.1  01         AL         Alabama   
 5 -88.5  32.0     5 FALSE 1     01.1  01         AL         Alabama   
 6 -88.5  32.1     6 FALSE 1     01.1  01         AL         Alabama   
 7 -88.4  32.2     7 FALSE 1     01.1  01         AL         Alabama   
 8 -88.4  32.2     8 FALSE 1     01.1  01         AL         Alabama   
 9 -88.4  32.2     9 FALSE 1     01.1  01         AL         Alabama   
10 -88.4  32.3    10 FALSE 1     01.1  01         AL         Alabama 

> map_data <- left_join(df, states, by = c("FIPS" = "state_fips"))
> map_data %>%
   ggplot(aes(long, lat, group = group, fill = GDPr)) +
   facet_wrap("YEAR") + 
   geom_polygon(color = NA) +
   coord_map(projection = "albers", lat0 = 39, lat1 = 45) +
   labs(fill = "")

enter image description here `

> WGScoor <- subset(map_data, YEAR == 2010)

> coordinates(WGScoor)=~long+lat

> proj4string(WGScoor)<- CRS("+proj=longlat +datum=WGS84")

> LLcoor<-spTransform(WGScoor,CRS("+proj=longlat"))

> raster::shapefile(LLcoor, filename = "Folder/Shapefile.shp", overwrite = T)

> SP2010 <- readOGR("Folder", layer = "Shapefile")
> SP2010
class       : SpatialPointsDataFrame 
features    : 72247 
extent      : -124.7332, -66.94989, 24.51496, 49.38436  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +ellps=WGS84 +no_defs 
variables   : 19
names       : YEAR, CENSUS_, STATE_F, STATE, FIPS,...

> poly2nb(SP2010, queen = T)
Error in poly2nb(SP2010, queen = T) : Not a polygon object

`

2

At this point:

 > WGScoor <- subset(map_data, YEAR == 2010)

map_data is a data frame with one point per vertex and a grouping variable, like piece, that defines the polygons. You then do:

 > coordinates(WGScoor)=~long+lat

which turns it into a spatial object, one spatial feature per row, and each row is a point. Its not polygons, because the sp package that defines these objects doesn't work with the polygons defined in ggplots "one row per point with a grouping variable" form.

This urbnmapr package doesn't appear on CRAN at this time. Does it have objects in sp format you can work with? You should be able to find a map of the USA states with matching FIPS codes and state abbreviations from a number of sources. Try the USABoundaries package (on CRAN).

Further, I'd suggest working with sf class objects instead of either sp or the needlessly bloated ggplot formats. sf objects behave a lot more like data frames when joining, and you can use packages such as tmap to make really good maps.

| improve this answer | |
1

The data in states <- data(urbnmapr::states) are not SpatialPolygons - they're just the SpatialPoints info to build the polygons if you wanted to (which apparently ggplot does automatically for you).

Look at the sample data you shared - there are 10 listings just for Alabama. If each row was a polygon, you should just see one row per state. The listings are of spatial points (lat/long) along with the order that they go in if you want use them to build a polygon, which you can do with your preferred spatial package.

| improve this answer | |
  • 1
    states isn't strictly even SpatialPoints - its just a dataframe with lat and long columns. And no coordinate system reference. WGScoor magically becomes a SpatialPoints object when it gets coordinates(WGScoor) assigned to it. – Spacedman May 13 at 22:08

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.