映射状态边界ggplot

问题描述 投票:0回答:1

我试图将状态边界映射到ggplot中的栅格(显示高程)。我得到一个有趣的错误,我添加我的geom_polygon时无法克服

 dir=("YOUR DIR")
library (tidyverse)
library (raster)
library (rasterVis)
library (rgeos)
require ("maps")
usa2 <- getData('GADM', country='USA', level=2)
counties<- c("El Paso","Jeff Davis","Culberson", "Hudspeth", "Persidio", "Brewster", "Pecos", "Loving", "Winkler", "Ward", "Lea", "Eddy", "Chaves", "Lincoln","Otero", "Dona Ana", "Socorro", "Sierra", "Luna", "Catron", "Grant", "Hidalgo", "Cochise", "Greenlee", "Graham", "Gila", "Pinal", "Pima", "Santa Cruz")
usa2 <- subset(usa2,NAME_2 %in% counties)

usa2c <- gCentroid(usa2) %>% coordinates()

dem2 <- getData("SRTM",lat=usa2c[2],lon=usa2c[1],path=dir)
dem=dem2
dem <- crop(dem,usa2,filename=file.path(dir,"dem_usa2.tif"))

# reduce raster resolution
dem_lower_res<- aggregate(dem, fact=10)

dem.p  <-  rasterToPoints(dem_lower_res )
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")


state<-getData("GADM", country="USA", level=2)
st.names<- c("Texas", "New Mexico", "Arizona")
state2<-subset(state,NAME_1 %in% st.names) %>% subset (NAME_2 %in% counties)
st.gg<- fortify(state2)

usa2 <- getData('GADM', country='USA', level=2)
p2 <- ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt)) + scale_fill_distiller(palette="BrBG") + theme_dark() +    
  theme(legend.position="bottom") + geom_path (aes (st.gg$long, st.gg$lat, group=st.gg$group))
p2

shapefile of counties of interest raster of altitude错误:是Error: mapping must be created by aes()

在提出建议后,我认为我们更接近:e

r ggplot2 raster projection shapefile
1个回答
1
投票

我无法获得SRTM数据,所以我使用了更好的name = "alt"。我也使用了sf包和geom_sf,因为它使事情变得更容易

dir=(getwd())
library (tidyverse)
library (raster)
library (rasterVis)
library (rgeos)
library(sf)
require ("maps")
usa2 <- getData('GADM', country='USA', level=2)
counties<- c("El Paso","Jeff Davis","Culberson", "Hudspeth", "Persidio", "Brewster",
"Pecos", "Loving", "Winkler", "Ward", "Lea", "Eddy", "Chaves", "Lincoln","Otero",
"Dona Ana", "Socorro", "Sierra", "Luna", "Catron", "Grant", "Hidalgo", "Cochise", 
"Greenlee", "Graham", "Gila", "Pinal", "Pima", "Santa Cruz")

usa2 <- subset(usa2,NAME_2 %in% counties)

state<-getData("GADM", country="USA", level=2)
st.names<- c("Texas", "New Mexico", "Arizona")
state2<-subset(state,NAME_1 %in% st.names) %>% subset (NAME_2 %in% counties)

dem2 <- getData("alt", country = "USA")
dem=dem2
dem <- crop(dem[[1]], state2, filename=file.path(dir,"dem_usa2.tif"), overwrite = TRUE)

# reduce raster resolution
dem_lower_res<- aggregate(dem, fact=10)

dem.p  <-  rasterToPoints(dem_lower_res )
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

state3 <- sf::st_as_sf(state2)

p2 <- ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt)) + 
  scale_fill_distiller(palette="BrBG") + 
  theme_dark() +    
  theme(legend.position="bottom") + 
  geom_sf(data = state3, inherit.aes = FALSE, fill = NA)
p2
© www.soinside.com 2019 - 2024. All rights reserved.