如何将城市点映射到具有偏移坐标的美国地图以允许区域之间 space ？
how to map cities points to US map with shifted coordinates to allow for space between regions?
我想使用 ggplot2 绘制美国地图，其中地图分为 4 个区域中的 1 个，每个区域都有空格 b/w。此外，我有一组城市坐标，我想映射到每个区域。我的问题如下。我可以很好地创建地图，但我无法让城市坐标点落在地图上。我知道在区域之间添加空间需要更改地图的坐标，但我也相应地更改了城市的坐标，因此我认为它们会在另一个区域移动，但整个过程一团糟...
library(maps)
library(ggplot2)
us.map <- map_data('state')
# add map regions
us.map$PADD[us.map$region %in%
c("connecticut", "maine", "massachusetts", "new hampshire", "rhode island", "vermont", "new jersey", "new york", "pennsylvania")] <- "PADD 1: East Coast"
us.map$PADD[us.map$region %in%
c("illinois", "indiana", "michigan", "ohio", "wisconsin", "iowa", "kansas", "minnesota", "missouri", "nebraska", "north dakota", "south dakota")] <- "PADD 2: Midwest"
us.map$PADD[us.map$region %in%
c("delaware", "florida", "georgia", "maryland", "north carolina", "south carolina", "virginia", "district of columbia", "west virginia", "alabama", "kentucky", "mississippi", "tennessee", "arkansas", "louisiana", "oklahoma", "texas")] <- "PADD 3: Gulf Coast"
us.map$PADD[us.map$region %in%
c("alaska", "california", "hawaii", "oregon", "washington", "arizona", "colorado", "idaho", "montana", "nevada", "new mexico", "utah", "wyoming")] <- "PADD 4: West Coast"
# subset the dataframe by region (PADD) and move lat/lon accordingly
us.map$lat.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$lat[us.map$PADD == "PADD 1: East Coast"]
us.map$long.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$long[us.map$PADD == "PADD 1: East Coast"] + 5
us.map$lat.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$lat[us.map$PADD == "PADD 2: Midwest"]
us.map$long.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$long[us.map$PADD == "PADD 2: Midwest"]
us.map$lat.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$lat[us.map$PADD == "PADD 3: Gulf Coast"] - 3
us.map$long.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$long[us.map$PADD == "PADD 3: Gulf Coast"]
us.map$lat.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$lat[us.map$PADD == "PADD 4: West Coast"] - 2
us.map$long.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$long[us.map$PADD == "PADD 4: West Coast"] - 10
# plot
ggplot(us.map, aes(x=long.transp, y=lat.transp), colour="white") +
geom_polygon(aes(group = group, fill="red")) +
theme(panel.background = element_blank(), # remove background
panel.grid = element_blank(),
axis.line = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()) +
coord_equal()+ scale_fill_manual(values="lightgrey", guide=FALSE)
结果如下：
这很好（一些代码来自：https://gis.stackexchange.com/questions/141181/how-to-create-a-us-map-in-r-with-separation-between-states-and-clear-labels）但我想将一组坐标映射到它。
Link 用于两个压缩数据集，下面使用 cities2.csv 和 PADDS.csv：
https://www.dropbox.com/s/zh9xyiakeuhgmdy/Archive.zip?dl=0（抱歉，数据太大无法使用dput输入）
#Two datasets found on dropbox link in zip
cities<-read.csv("cities2.csv")
padds<-read.csv("PADDS.csv")
padds$State<-NULL
colnames(padds)<-c("state","PADD")
points<-merge(cities, padds, by="state",all.x=TRUE)
#Shift city coordinates according to padd region
points$Long2<-ifelse(points$PADD =="PADD 1: East Coast", points$Long+5, points$Long)
points$Long2<-ifelse(points$PADD =="PADD 4: West Coast", points$Long-10, points$Long2)
points$Lat2<-ifelse(points$PADD =="PADD 3: Gulf Coast", points$Lat-3, points$Lat)
points$Lat2<-ifelse(points$PADD =="PADD 4: West Coast", points$Lat-2, points$Lat2)
结果如下：
显然这里出了点问题...非常感谢任何帮助。
我认为您的 cities CSV 文件中的坐标有误。下面是我如何检查坐标。我首先下载了您的 CSV 文件，将文件读取为 cities，然后我创建了一个 sf 对象并使用 mapview 包将其可视化。
colnames(cities) <- c("state", "Lat", "Long")
library(sf)
library(mapview)
cities_sf <- cities %>%
st_as_sf(coords = c("Long", "Lat"), crs = 4326)
mapview(cities_sf)
如你所见，纬度似乎是对的，但经度都错了。但是，您似乎只是弄错了经度符号，因为我仍然可以根据这些点看到美国的形状。
所以，这是一个快速修复。
library(dplyr)
cities2 <- cities %>% mutate(Long = -Long)
cities_sf2 <- cities2 %>%
st_as_sf(coords = c("Long", "Lat"), crs = 4326)
mapview(cities_sf2)
现在cities2中的坐标是正确的。所以我们可以运行你的代码来绘制城市地图。
colnames(padds)<-c("state","PADD")
points<-merge(cities2, padds, by="state",all.x=TRUE)
points$Long2<-ifelse(points$PADD %in% "PADD 1: East Coast", points$Long+5, points$Long)
points$Long2<-ifelse(points$PADD %in% "PADD 4: West Coast", points$Long-10, points$Long2)
points$Lat2<-ifelse(points$PADD %in% "PADD 3: Gulf Coast", points$Lat-3, points$Lat)
points$Lat2<-ifelse(points$PADD %in% "PADD 4: West Coast", points$Lat-2, points$Lat2)
# P is the ggplot object you created earlier
P + geom_point(data = points, aes(x = Long2, y = Lat2))
更新
这是 OP 要求的完整代码。
library(maps)
library(ggplot2)
library(dplyr)
#Two datasets found on dropbox link in zip
cities<-read.csv("cities.csv")
padds<-read.csv("PADDS.csv")
padds$State<-NULL
colnames(cities) <- c("state", "Lat", "Long")
colnames(padds)<-c("state","PADD")
cities2 <- cities %>% mutate(Long = -Long)
points<-merge(cities2, padds, by="state",all.x=TRUE)
#Shift city coordinates according to padd region
points$Long2<-ifelse(points$PADD =="PADD 1: East Coast", points$Long+5, points$Long)
points$Long2<-ifelse(points$PADD =="PADD 4: West Coast", points$Long-10, points$Long2)
points$Lat2<-ifelse(points$PADD =="PADD 3: Gulf Coast", points$Lat-3, points$Lat)
points$Lat2<-ifelse(points$PADD =="PADD 4: West Coast", points$Lat-2, points$Lat2)
us.map <- map_data('state')
# add map regions
us.map$PADD[us.map$region %in%
c("connecticut", "maine", "massachusetts", "new hampshire", "rhode island", "vermont", "new jersey", "new york", "pennsylvania")] <- "PADD 1: East Coast"
us.map$PADD[us.map$region %in%
c("illinois", "indiana", "michigan", "ohio", "wisconsin", "iowa", "kansas", "minnesota", "missouri", "nebraska", "north dakota", "south dakota")] <- "PADD 2: Midwest"
us.map$PADD[us.map$region %in%
c("delaware", "florida", "georgia", "maryland", "north carolina", "south carolina", "virginia", "district of columbia", "west virginia", "alabama", "kentucky", "mississippi", "tennessee", "arkansas", "louisiana", "oklahoma", "texas")] <- "PADD 3: Gulf Coast"
us.map$PADD[us.map$region %in%
c("alaska", "california", "hawaii", "oregon", "washington", "arizona", "colorado", "idaho", "montana", "nevada", "new mexico", "utah", "wyoming")] <- "PADD 4: West Coast"
# subset the dataframe by region (PADD) and move lat/lon accordingly
us.map$lat.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$lat[us.map$PADD == "PADD 1: East Coast"]
us.map$long.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$long[us.map$PADD == "PADD 1: East Coast"] + 5
us.map$lat.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$lat[us.map$PADD == "PADD 2: Midwest"]
us.map$long.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$long[us.map$PADD == "PADD 2: Midwest"]
us.map$lat.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$lat[us.map$PADD == "PADD 3: Gulf Coast"] - 3
us.map$long.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$long[us.map$PADD == "PADD 3: Gulf Coast"]
us.map$lat.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$lat[us.map$PADD == "PADD 4: West Coast"] - 2
us.map$long.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$long[us.map$PADD == "PADD 4: West Coast"] - 10
# plot
P <- ggplot(us.map, aes(x=long.transp, y=lat.transp), colour="white") +
geom_polygon(aes(group = group, fill="red")) +
theme(panel.background = element_blank(), # remove background
panel.grid = element_blank(),
axis.line = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()) +
coord_equal()+ scale_fill_manual(values="lightgrey", guide=FALSE)
# P is the ggplot object you created earlier
P + geom_point(data = points, aes(x = Long2, y = Lat2))
有趣的问题！这是一个基于 sf 包的解决方案，它还可能使将这样的图与其他空间数据结合起来变得更加容易。方法是：
- 从不同的包
USAboundaries::us_states() 而不是 ggplot2::map_data 获取州边界，因为将单个点转换为多边形会浪费时间
- 使用
st_as_sf将城市坐标转换为点，设置坐标系，并修复错误的坐标符号，如其他答案所述。 （N.B。手动从 colname 中删除了非标准字符）
- 合并文件，添加区域标签，并删除阿拉斯加、夏威夷和波多黎各（因为它们会使情节看起来很奇怪，而你没有使用它们）
- 使用
st_set_geometry 将形状替换为所需的翻译形状（只需执行 + c(x, y)）。请注意，sf 使用 (x, y) 进行仿射变换，即 (long, lat).
- 使用
geom_sf绘制点和形状。
我认为这种方法的主要优点是您可以随后使用您喜欢的来自 sf 的任何空间工具，并且代码可能更具可读性。如果您要制作类似的图，这可能是值得的。主要缺点可能是需要额外的包，包括 ggplot2 的开发版本以获得 geom_sf（使用 devtools::install_github("tidyverse/ggplot2" 来安装它）。这也比仅仅将经度更改为负数并使用现有代码要多得多...
library(tidyverse)
library(sf)
library(USAboundaries)
# Define regions
padd1 <- c("CT", "ME", "MA", "NH", "RI", "VT", "NJ", "NY", "PA")
padd2 <- c("IL", "IN", "MI", "OH", "WI", "IA", "KS", "MN", "MO", "NE", "ND",
"SD")
padd3 <- c("DE", "FL", "GA", "MD", "NC", "SC", "VA", "DC", "WV", "AL", "KY",
"MS", "TN", "AR", "LA", "OK", "TX")
padd4 <- c("AK", "CA", "HI", "OR", "WA", "AZ", "CO", "ID", "MT", "NV", "NM",
"UT", "WY")
us_map <- us_states() %>%
select(state_abbr) # keep only state abbreviation column
cities <- read_csv(here::here("data", "cities.csv")) %>%
mutate(Long = -Long) %>% # make longitudes negative
st_as_sf(coords = 3:2) %>% # turn into sf object
st_set_crs(4326) %>% # add coordinate system
rename(state_abbr = StateAbbr)
combined <- rbind(us_map, cities) %>%
filter(!(state_abbr %in% c("AK", "HI", "PR"))) %>% # remove non-contiguous cities and states
mutate( # add region identifier based on state
region = case_when(
state_abbr %in% padd1 ~ "PADD 1: East Coast",
state_abbr %in% padd2 ~ "PADD 2: Midwest",
state_abbr %in% padd3 ~ "PADD 3: Gulf Coast",
state_abbr %in% padd4 ~ "PADD 4: West Coast"
)
)
eastc <- combined %>%
filter(region == "PADD 1: East Coast") %>%
st_set_geometry(., .$geometry + c(5, 0)) # replace geometries with 5 degrees east
mwest <- combined %>%
filter(region == "PADD 2: Midwest") %>%
st_set_geometry(., .$geometry + c(0, 0))
gulfc <- combined %>%
filter(region == "PADD 3: Gulf Coast") %>%
st_set_geometry(., .$geometry + c(0, -3))
westc <- combined %>%
filter(region == "PADD 4: West Coast") %>%
st_set_geometry(., .$geometry + c(-10, -2))
ggplot(data = rbind(eastc, mwest, gulfc, westc)) + # bind regions together
theme_bw() +
geom_sf(aes(fill = region))
这是输出图
我想使用 ggplot2 绘制美国地图，其中地图分为 4 个区域中的 1 个，每个区域都有空格 b/w。此外，我有一组城市坐标，我想映射到每个区域。我的问题如下。我可以很好地创建地图，但我无法让城市坐标点落在地图上。我知道在区域之间添加空间需要更改地图的坐标，但我也相应地更改了城市的坐标，因此我认为它们会在另一个区域移动，但整个过程一团糟...
library(maps)
library(ggplot2)
us.map <- map_data('state')
# add map regions
us.map$PADD[us.map$region %in%
c("connecticut", "maine", "massachusetts", "new hampshire", "rhode island", "vermont", "new jersey", "new york", "pennsylvania")] <- "PADD 1: East Coast"
us.map$PADD[us.map$region %in%
c("illinois", "indiana", "michigan", "ohio", "wisconsin", "iowa", "kansas", "minnesota", "missouri", "nebraska", "north dakota", "south dakota")] <- "PADD 2: Midwest"
us.map$PADD[us.map$region %in%
c("delaware", "florida", "georgia", "maryland", "north carolina", "south carolina", "virginia", "district of columbia", "west virginia", "alabama", "kentucky", "mississippi", "tennessee", "arkansas", "louisiana", "oklahoma", "texas")] <- "PADD 3: Gulf Coast"
us.map$PADD[us.map$region %in%
c("alaska", "california", "hawaii", "oregon", "washington", "arizona", "colorado", "idaho", "montana", "nevada", "new mexico", "utah", "wyoming")] <- "PADD 4: West Coast"
# subset the dataframe by region (PADD) and move lat/lon accordingly
us.map$lat.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$lat[us.map$PADD == "PADD 1: East Coast"]
us.map$long.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$long[us.map$PADD == "PADD 1: East Coast"] + 5
us.map$lat.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$lat[us.map$PADD == "PADD 2: Midwest"]
us.map$long.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$long[us.map$PADD == "PADD 2: Midwest"]
us.map$lat.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$lat[us.map$PADD == "PADD 3: Gulf Coast"] - 3
us.map$long.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$long[us.map$PADD == "PADD 3: Gulf Coast"]
us.map$lat.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$lat[us.map$PADD == "PADD 4: West Coast"] - 2
us.map$long.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$long[us.map$PADD == "PADD 4: West Coast"] - 10
# plot
ggplot(us.map, aes(x=long.transp, y=lat.transp), colour="white") +
geom_polygon(aes(group = group, fill="red")) +
theme(panel.background = element_blank(), # remove background
panel.grid = element_blank(),
axis.line = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()) +
coord_equal()+ scale_fill_manual(values="lightgrey", guide=FALSE)
结果如下：
这很好（一些代码来自：https://gis.stackexchange.com/questions/141181/how-to-create-a-us-map-in-r-with-separation-between-states-and-clear-labels）但我想将一组坐标映射到它。
Link 用于两个压缩数据集，下面使用 cities2.csv 和 PADDS.csv：
https://www.dropbox.com/s/zh9xyiakeuhgmdy/Archive.zip?dl=0（抱歉，数据太大无法使用dput输入）
#Two datasets found on dropbox link in zip
cities<-read.csv("cities2.csv")
padds<-read.csv("PADDS.csv")
padds$State<-NULL
colnames(padds)<-c("state","PADD")
points<-merge(cities, padds, by="state",all.x=TRUE)
#Shift city coordinates according to padd region
points$Long2<-ifelse(points$PADD =="PADD 1: East Coast", points$Long+5, points$Long)
points$Long2<-ifelse(points$PADD =="PADD 4: West Coast", points$Long-10, points$Long2)
points$Lat2<-ifelse(points$PADD =="PADD 3: Gulf Coast", points$Lat-3, points$Lat)
points$Lat2<-ifelse(points$PADD =="PADD 4: West Coast", points$Lat-2, points$Lat2)
结果如下：
显然这里出了点问题...非常感谢任何帮助。
我认为您的 cities CSV 文件中的坐标有误。下面是我如何检查坐标。我首先下载了您的 CSV 文件，将文件读取为 cities，然后我创建了一个 sf 对象并使用 mapview 包将其可视化。
colnames(cities) <- c("state", "Lat", "Long")
library(sf)
library(mapview)
cities_sf <- cities %>%
st_as_sf(coords = c("Long", "Lat"), crs = 4326)
mapview(cities_sf)
如你所见，纬度似乎是对的，但经度都错了。但是，您似乎只是弄错了经度符号，因为我仍然可以根据这些点看到美国的形状。
所以，这是一个快速修复。
library(dplyr)
cities2 <- cities %>% mutate(Long = -Long)
cities_sf2 <- cities2 %>%
st_as_sf(coords = c("Long", "Lat"), crs = 4326)
mapview(cities_sf2)
现在cities2中的坐标是正确的。所以我们可以运行你的代码来绘制城市地图。
colnames(padds)<-c("state","PADD")
points<-merge(cities2, padds, by="state",all.x=TRUE)
points$Long2<-ifelse(points$PADD %in% "PADD 1: East Coast", points$Long+5, points$Long)
points$Long2<-ifelse(points$PADD %in% "PADD 4: West Coast", points$Long-10, points$Long2)
points$Lat2<-ifelse(points$PADD %in% "PADD 3: Gulf Coast", points$Lat-3, points$Lat)
points$Lat2<-ifelse(points$PADD %in% "PADD 4: West Coast", points$Lat-2, points$Lat2)
# P is the ggplot object you created earlier
P + geom_point(data = points, aes(x = Long2, y = Lat2))
更新
这是 OP 要求的完整代码。
library(maps)
library(ggplot2)
library(dplyr)
#Two datasets found on dropbox link in zip
cities<-read.csv("cities.csv")
padds<-read.csv("PADDS.csv")
padds$State<-NULL
colnames(cities) <- c("state", "Lat", "Long")
colnames(padds)<-c("state","PADD")
cities2 <- cities %>% mutate(Long = -Long)
points<-merge(cities2, padds, by="state",all.x=TRUE)
#Shift city coordinates according to padd region
points$Long2<-ifelse(points$PADD =="PADD 1: East Coast", points$Long+5, points$Long)
points$Long2<-ifelse(points$PADD =="PADD 4: West Coast", points$Long-10, points$Long2)
points$Lat2<-ifelse(points$PADD =="PADD 3: Gulf Coast", points$Lat-3, points$Lat)
points$Lat2<-ifelse(points$PADD =="PADD 4: West Coast", points$Lat-2, points$Lat2)
us.map <- map_data('state')
# add map regions
us.map$PADD[us.map$region %in%
c("connecticut", "maine", "massachusetts", "new hampshire", "rhode island", "vermont", "new jersey", "new york", "pennsylvania")] <- "PADD 1: East Coast"
us.map$PADD[us.map$region %in%
c("illinois", "indiana", "michigan", "ohio", "wisconsin", "iowa", "kansas", "minnesota", "missouri", "nebraska", "north dakota", "south dakota")] <- "PADD 2: Midwest"
us.map$PADD[us.map$region %in%
c("delaware", "florida", "georgia", "maryland", "north carolina", "south carolina", "virginia", "district of columbia", "west virginia", "alabama", "kentucky", "mississippi", "tennessee", "arkansas", "louisiana", "oklahoma", "texas")] <- "PADD 3: Gulf Coast"
us.map$PADD[us.map$region %in%
c("alaska", "california", "hawaii", "oregon", "washington", "arizona", "colorado", "idaho", "montana", "nevada", "new mexico", "utah", "wyoming")] <- "PADD 4: West Coast"
# subset the dataframe by region (PADD) and move lat/lon accordingly
us.map$lat.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$lat[us.map$PADD == "PADD 1: East Coast"]
us.map$long.transp[us.map$PADD == "PADD 1: East Coast"] <- us.map$long[us.map$PADD == "PADD 1: East Coast"] + 5
us.map$lat.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$lat[us.map$PADD == "PADD 2: Midwest"]
us.map$long.transp[us.map$PADD == "PADD 2: Midwest"] <- us.map$long[us.map$PADD == "PADD 2: Midwest"]
us.map$lat.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$lat[us.map$PADD == "PADD 3: Gulf Coast"] - 3
us.map$long.transp[us.map$PADD == "PADD 3: Gulf Coast"] <- us.map$long[us.map$PADD == "PADD 3: Gulf Coast"]
us.map$lat.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$lat[us.map$PADD == "PADD 4: West Coast"] - 2
us.map$long.transp[us.map$PADD == "PADD 4: West Coast"] <- us.map$long[us.map$PADD == "PADD 4: West Coast"] - 10
# plot
P <- ggplot(us.map, aes(x=long.transp, y=lat.transp), colour="white") +
geom_polygon(aes(group = group, fill="red")) +
theme(panel.background = element_blank(), # remove background
panel.grid = element_blank(),
axis.line = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
axis.text = element_blank()) +
coord_equal()+ scale_fill_manual(values="lightgrey", guide=FALSE)
# P is the ggplot object you created earlier
P + geom_point(data = points, aes(x = Long2, y = Lat2))
有趣的问题！这是一个基于 sf 包的解决方案，它还可能使将这样的图与其他空间数据结合起来变得更加容易。方法是：
- 从不同的包
USAboundaries::us_states()而不是ggplot2::map_data获取州边界，因为将单个点转换为多边形会浪费时间 - 使用
st_as_sf将城市坐标转换为点，设置坐标系，并修复错误的坐标符号，如其他答案所述。 （N.B。手动从 colname 中删除了非标准字符） - 合并文件，添加区域标签，并删除阿拉斯加、夏威夷和波多黎各（因为它们会使情节看起来很奇怪，而你没有使用它们）
- 使用
st_set_geometry将形状替换为所需的翻译形状（只需执行+ c(x, y)）。请注意，sf使用(x, y)进行仿射变换，即(long, lat). - 使用
geom_sf绘制点和形状。
我认为这种方法的主要优点是您可以随后使用您喜欢的来自 sf 的任何空间工具，并且代码可能更具可读性。如果您要制作类似的图，这可能是值得的。主要缺点可能是需要额外的包，包括 ggplot2 的开发版本以获得 geom_sf（使用 devtools::install_github("tidyverse/ggplot2" 来安装它）。这也比仅仅将经度更改为负数并使用现有代码要多得多...
library(tidyverse)
library(sf)
library(USAboundaries)
# Define regions
padd1 <- c("CT", "ME", "MA", "NH", "RI", "VT", "NJ", "NY", "PA")
padd2 <- c("IL", "IN", "MI", "OH", "WI", "IA", "KS", "MN", "MO", "NE", "ND",
"SD")
padd3 <- c("DE", "FL", "GA", "MD", "NC", "SC", "VA", "DC", "WV", "AL", "KY",
"MS", "TN", "AR", "LA", "OK", "TX")
padd4 <- c("AK", "CA", "HI", "OR", "WA", "AZ", "CO", "ID", "MT", "NV", "NM",
"UT", "WY")
us_map <- us_states() %>%
select(state_abbr) # keep only state abbreviation column
cities <- read_csv(here::here("data", "cities.csv")) %>%
mutate(Long = -Long) %>% # make longitudes negative
st_as_sf(coords = 3:2) %>% # turn into sf object
st_set_crs(4326) %>% # add coordinate system
rename(state_abbr = StateAbbr)
combined <- rbind(us_map, cities) %>%
filter(!(state_abbr %in% c("AK", "HI", "PR"))) %>% # remove non-contiguous cities and states
mutate( # add region identifier based on state
region = case_when(
state_abbr %in% padd1 ~ "PADD 1: East Coast",
state_abbr %in% padd2 ~ "PADD 2: Midwest",
state_abbr %in% padd3 ~ "PADD 3: Gulf Coast",
state_abbr %in% padd4 ~ "PADD 4: West Coast"
)
)
eastc <- combined %>%
filter(region == "PADD 1: East Coast") %>%
st_set_geometry(., .$geometry + c(5, 0)) # replace geometries with 5 degrees east
mwest <- combined %>%
filter(region == "PADD 2: Midwest") %>%
st_set_geometry(., .$geometry + c(0, 0))
gulfc <- combined %>%
filter(region == "PADD 3: Gulf Coast") %>%
st_set_geometry(., .$geometry + c(0, -3))
westc <- combined %>%
filter(region == "PADD 4: West Coast") %>%
st_set_geometry(., .$geometry + c(-10, -2))
ggplot(data = rbind(eastc, mwest, gulfc, westc)) + # bind regions together
theme_bw() +
geom_sf(aes(fill = region))
这是输出图