给定这个可重现的示例
# 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)。
谢谢你的帮助
发布于 2021-06-28 21:41:44
在漫长而曲折的道路之后,我设法找到了我在这里张贴的解决方案
# 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))https://stackoverflow.com/questions/66832903
复制相似问题