首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >R ggmap覆盖geom_contour

R ggmap覆盖geom_contour
EN

Stack Overflow用户
提问于 2021-03-27 23:18:31
回答 1查看 166关注 0票数 0

给定这个可重现的示例

代码语言:javascript
运行
复制
# set the origin of the grid 
# in cartesian coordinates (epsg 32632)
xmin<-742966
ymin<-5037923

# set x and y axis
x<-seq(xmin, xmin+25*39, by=25)
y<-seq(ymin, ymin+25*39, by =25)

# define a 40 x 40 grid
mygrid<-expand.grid(x = x, y = y)

# set the z value to be interpolated by the contour
set.seed(123)
mygrid$z<- rnorm(nrow(mygrid))

library(tidyverse)

# plot of contour is fine
ggplot(data=mygrid, aes(x=x,y=y,z=z))+
  geom_contour()

library(ggspatial)

# transform coordinates to wgs84 4326
# (one of the possible many other ways to do it)
mygrid_4326<-xy_transform(mygrid$x, mygrid$y, from = 32632, to = 4326)

# create new grid with lon and lat 
# (geographical coordinates espg 4326)
mygrid_4326<-mygrid_4326%>%
  mutate(z=mygrid$z)

# define the bounding box
my_bb<-c(min(mygrid_4326$x), min(mygrid_4326$y), 
         max(mygrid_4326$x), max(mygrid_4326$y))
names(my_bb)<-c('left', 'bottom', 'right', 'top')

library(ggmap)

# get the background map (by a free provider)
mymap<-get_stamenmap(bbox = c(left = my_bb[['left']], 
                              bottom = my_bb[['bottom']], 
                              right = my_bb[['right']], 
                              top = my_bb[['top']]), 
                     zoom = 15, 
                     maptype = 'toner-lite')

# plot of the map is fine
mymap%>%
  ggmap()

# overlay the contour of z is failing
mymap%>%
  ggmap()+
  #geom_contour(data=mygrid_4326, mapping=aes(x = x, y = y, z = z))
  stat_contour(data=mygrid_4326, mapping=aes(x = x, y = y, z = z))

Warning messages:
1: stat_contour(): Zero contours were generated 
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning -Inf

现在我正在寻找一种可行的解决方案来完成用ggplot制作的等高线图到用ggmap制作的底图的叠加;由于一些我不完全理解的原因,这似乎是不可能的;

这可能与坐标的变换有关,以某种方式影响了网格的形状(变得不完全规则)?

这篇文章:Plotting contours on map using ggmap似乎接近问题的核心,但还没有给出明确的解决方案(如果它存在的话)

如果可能,我希望留在ggplot (tidyverse)工具中,而不是求助于基础R的函数(例如contourLines)。

谢谢你的帮助

EN

回答 1

Stack Overflow用户

发布于 2021-06-28 21:41:44

在漫长而曲折的道路之后,我设法找到了我在这里张贴的解决方案

代码语言:javascript
运行
复制
# set the origin of the grid 
# in cartesian coordinates (epsg 32632)
xmin<-742966
ymin<-5037923

# set x and y axis
x<-seq(xmin, xmin+25*39, by=25)
y<-seq(ymin, ymin+25*39, by =25)

# define a 40 x 40 grid
mygrid<-expand.grid(x = x, y = y)

# set the z value to be interpolated by the contour
set.seed(123)
mygrid$z<- rnorm(nrow(mygrid))

library(sf)
# create simple feature in original crs
mygrid_sf <- st_as_sf(mygrid, coords = c("x", "y"), crs = 32632)

# transform to crs 4326
mygrid_sf_4326<-st_transform(mygrid_sf, 4326)    

# create contourlines
cL<-contourLines(x,
                 y,
                 matrix(mygrid$z, 40, 40, byrow = TRUE))

library(maptools)
# spatial line df
epsg_in<-'+init=epsg:32632'
cLdf<-ContourLines2SLDF(cL, proj4string=CRS(epsg_in))

# transform to sf
cLdf_sf <- st_as_sf(cLdf)

# reproject sf to 4326
cLdf_sf_4326<-st_transform(cLdf_sf, 4326)

# define the bounding box
min_lon<-min(st_coordinates(mygrid_sf_4326)[,1])
max_lon<-max(st_coordinates(mygrid_sf_4326)[,1])
min_lat<-min(st_coordinates(mygrid_sf_4326)[,2])
max_lat<-max(st_coordinates(mygrid_sf_4326)[,2])

library(ggplot2) 
library(ggmap)

# get the map
mymap<-get_stamenmap(bbox = c(left = min_lon, 
                              bottom = min_lat, 
                              right = max_lon, 
                              top = max_lat), 
                     zoom = 15, 
                     maptype = 'toner-lite')


# this aproach is failing
mymap%>%
  ggmap()+
  geom_sf(data=cLdf_sf_4326, aes(colour=level))

#Adding new coordinate system, which will replace the existing one.
#Error in FUN(X[[i]], ...) : object 'lon' not found

# but this is fine, because of this hack
# https://github.com/r-spatial/sf/issues/336

mymap%>%
  ggmap()+
  geom_sf(data=cLdf_sf_4326, aes(color=level), inherit.aes = FALSE)


# and finally this is another, cleaner, approach using ggspatial facilities

library(ggspatial)

ggplot() +
  annotation_map_tile(type = "osm") +
  geom_sf(data=cLdf_sf_4326, aes(color=level))

# rosm::osm.types()

ggplot() +
  annotation_map_tile(type = "stamenbw") +
  geom_sf(data=cLdf_sf_4326, aes(color=level))
票数 1
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/66832903

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档