Add territories to map

[alicefeng]

It'd be great if we could map territories in addition to the states. I just got asked to include Puerto Rico in a few state-level maps.

[benjaminrobinson]

Hi! I really like the project as I need to have the combination of being able to map states as well as territories in my work. I'm having some trouble mapping the territories and laying a geom_polygon layer of data on top of it and I think I know why its happening. I looked at your data-raw/utils.R script and I'm not sure if you're getting the territory data to be in the same projection as the states.

Here's the straight mapping of Alaska:

Here's the straight mapping of the geom_polygon:

Here's the combination:

I think this arises from a failed projection in the get_shapefile function, specifically in lines 109-114 where you accidentally rbind the alaska object twice:

mapdata <-
    shapes[!shapes$STATEFP %in% exclude, ] %>%
    rbind(alaska, hawaii, alaska, guam, puerto_rico, virgin_islands,
          mariana_islands, american_samoa) %>%
    # convert from EPSG2163 to (US National Atlas Equal Area) WGS84
    spTransform(CRS("+init=epsg:4326"))

Either way, I believe the projection for the territories is off from what it should be (i.e. in alignment with the state shapes like you intended to do). Thanks for the project! I'm looking forward to using it a lot in my spatial R work!

• Ben.

[khueyama]

Hi @benjaminrobinson! Glad to hear that you like the project, and sorry that you're having issues. Would you be able to pass along the code you used to create these plots so that we can take a deeper dive into things?

[benjaminrobinson]

Sure! Happy to oblige:

library(RJSONIO)
library(dplyr)
library(ggplot2)
library(urbnmapr)

url <- "https://services.arcgis.com/YnOQrIGdN9JGtBh4/arcgis/rest/services/Public_Study_Areas/FeatureServer/0/query?where=SAC=613023&f=json&outSR=4326"
len <- url %>%
        RJSONIO::fromJSON() %>%
        .$features %>%
        first %>%
        .$geometry %>%
        .$rings %>%
        lapply(., length) %>%
        unlist

sac <- url %>%
          RJSONIO::fromJSON() %>%
          .$features %>%
          first %>%
          .$geometry %>%
          .$rings %>%
          unlist(recursive = FALSE) %>%
          do.call(rbind, .) %>%
          as.data.frame %>%
          setNames(c("long", "lat")) %>%
          mutate(sac = '613023',
                 id = rep(1:length(len),  times = len))

ggplot() +
            geom_polygon(data = urbnmapr::territories_counties %>% filter(state_abbv == 'AK'), 
            aes(x = long, y = lat, group = group), fill = "#E0E0E0", color = "black") +
            geom_polygon(data = sac , 
            aes(x = long, y = lat, group = id), fill = "#FF7500") + 
            coord_map("mercator") + theme_void() + guides(fill = FALSE) + 
            ggtitle("ALASKA | SAC:  613023 - UNITED UTILITIES INC")

[benjaminrobinson]

And if you do it for a state in the states portion of urbnmapr:

library(RJSONIO)
library(dplyr)
library(ggplot2)
library(urbnmapr)

url <- "https://services.arcgis.com/YnOQrIGdN9JGtBh4/arcgis/rest/services/Public_Study_Areas/FeatureServer/0/query?where=SAC=305062&f=json&outSR=4326"
len <- url %>%
        RJSONIO::fromJSON() %>%
        .$features %>%
        first %>%
        .$geometry %>%
        .$rings %>%
        lapply(., length) %>%
        unlist

sac <- url %>%
          RJSONIO::fromJSON() %>%
          .$features %>%
          first %>%
          .$geometry %>%
          .$rings %>%
          unlist(recursive = FALSE) %>%
          do.call(rbind, .) %>%
          as.data.frame %>%
          setNames(c("long", "lat")) %>%
          mutate(sac = '305062',
                 id = rep(1:length(len),  times = len))

ggplot() +
            geom_polygon(data = urbnmapr::states %>% filter(state_abbv == 'OH'), 
            aes(x = long, y = lat, group = group), fill = "#E0E0E0", color = "black") +
            geom_polygon(data = sac , 
            aes(x = long, y = lat, group = id), fill = "#FF7500") + 
            coord_map("mercator") + theme_void() + guides(fill = FALSE) + 
            ggtitle("OHIO | SAC: 305062 - CINCINNATI BELL-OH")

[benjaminrobinson]

Just wondering if you're seeing what I'm seeing. Thanks!

[khueyama]

@benjaminrobinson I believe the issue is that we've elided Alaska, Hawaii, and the territories so that they display as insets within the map of the continental US. So your data is plotting to where Alaska truly is, while the map itself has been shifted down and to the right. I think there should be a way to elide your data as well, but I haven't had much luck so far.

[benjaminrobinson]

Hmm...Interesting. Thanks for the explanation! I'll look into it on my end as well because I think this is a really interesting and appealing use case. Thanks again for your work on this!

[benjaminrobinson]

I also want to be able to do this type of mapping on demand in a shiny application so I'd ideally like to have some functionality that could be scaled depending on which of Hawaii, Alaska, or the territories I'm looking at.

[awunderground]

@benjaminrobinson You could try transform_state() from urbnmapr/data-raw/utils.R. Plug in your shapefiles with the same arguments we use for Alaska.

transform_state <- function(object, rot, scale, shift) {
  object %>% elide(rotate = rot) %>%
    elide(scale = max(apply(bbox(object), 1, diff)) / scale) %>%
    elide(shift = shift)
}

transform_state(-35, 2, c(-2600000, -2300000))

proj4string(alaska) <- proj4string(shapes)

I can't guarantee that this will work. The package is designed for making simple choropleths. I've mapped points on the maps, but I haven't mapped shapes.

@alicefeng "territories" and "territories_counties" can now be used with get_urbn_map() and get_urbn_labels(). Both functions are explained in the README.

Closed with 1135665312ee679881c4fd7eea628f13691e9df7

[alicefeng]

Thanks for updating the readme @awunderground !