0

given this reproducible example

# 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

Now I'm searching for a viable solution to accomplish the overlay of a contour plot made with ggplot to a base map made with ggmap; it seems not possible for some reasons I do not fully understand;

might it be related to the transformation of coordinates affecting somehow the shape of the grid (becoming not perfectly regular)?

This post: Plotting contours on map using ggmap seems close to the core of the problem but not yet giving a definitive solution (if it exists)

If possibile I would like to stay within the ggplot (tidyverse) facilities without resorting to the base R's functions (e.g. contourLines).

Thank you for your help

maxbre
  • 161
  • 9

1 Answers1

1

after a long and winding road I managed to find the solution that I'm here posting

# 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))
maxbre
  • 161
  • 9