How to generate high resolution temperature map using unevenly spaced coordinates in R?

Question

I want to generate 12 maps using ggplot and facet_wrap which illustrate ocean temperature across the Scotian shelf from Jan-Dec. I was given a ".csv" file with 25,000 observations, which has monthly temperature data from 2016-2018 with the associated latitude and longitude values. I have tried using geom_raster to map this data but I get the following error message

Error: cannot allocate vector of size 39.7 Gb
In addition: Warning messages:
1: In f(...) :
 Raster pixels are placed at uneven horizontal intervals and will be shifted. Consider using geom_tile() instead.
2: In f(...) :
 Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead

So of course, I have also tried using geom_tile but when I run the code I get a blank map with no color even though I have specified both fill and color. Here is some sample data similar to my data frame, please note that my real data for latitude and longitude values are not evenly spaced by seq(...,...,0.1) (I'm not sure how to make sample data without the sequencing, sorry) so you won't get the same error as me, and the code will work

#In my data frame the lat and long values are not equally spaced!
mapoc_temp = expand.grid(data.frame(Longitude= seq(-64.5,-62.4,0.1), 
                                    Latitude= seq(42.7,44.8,0.1),
                                    year = sample(c(2016,2017,2018), 22, replace = T),
                                    month = sample(month.abb, 22, replace = T)))
mapoc_temp$Temp = runif(nrow(mapoc_temp))

Here is my code for geom_raster which gives me the error I mentioned

library(mapdata)
library(ggplot2)

#map I'm using for ggplot
Canada = map_data("wildfires", "Canada")

ggplot(mapoc_temp, aes(x=Longitude, y=Latitude)) + 
#try to map my temp on the ocean
  geom_raster(aes(fill = Temp, x = Longitude), interpolate = TRUE) +
  geom_polygon(data = Canada, aes(x=long, y=lat, group=group), colour="grey50", fill = 'grey55')+
  coord_sf(xlim=c(-64.5,-62.8), ylim=c(42.7,45)) +
#get months
facet_wrap(vars(month))

here is the code when I try geom_tile, again I noticed with my sample data it works, but with my real data it doesn't and I believe it has something to do with my coordinates not being equal distances apart.

ggplot(mapoc_temp, aes(x=Longitude, y=Latitude)) + 
  geom_tile(aes(fill = Temp, x = Longitude, y = Latitude), colour = mapoc_temp$Temp) +
  geom_polygon(data = canada, aes(x=long, y=lat, group=group), colour="grey50", fill = 'grey55')+
  coord_sf(xlim=c(-64.5,-62.8), ylim=c(42.7,45)) +
  facet_wrap(vars(month))

here is the picture I get with geom_tile

and here is roughly what I am trying to produce, but obviously, 12 of these maps because I want a map for each month, with better resolution (if possible, but at this point, I'll take any map colored with my temperature data).

I have a feeling I'll have to do something more along the lines with this type of code (a snippet of something I found a while ago), but I have tried manipulating this with no success. Any suggestions?

#my CRS value!
latlong = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

#for transforming
planar ="+proj=utm +zone=20 +datum=NAD83 +units=km +no_defs +ellps=GRS80 +towgs84=0,0,0"


    out <- x%>%
      st_as_sf(coords=c("Longitude","Latitude"),crs=latlong)%>%
      st_transform(planar)%>%
      st_coordinates()%>%
      raster::rasterize(.,grid,field=xyz[,3],fun=rasterFun)%>%
      raster::projectRaster(.,crs=latlong)%>%
      raster::rasterToPolygons(.)%>% # this part is slow
      st_as_sf()%>%
      rename(MAP = layer)

Answer 1

You have sort of answered your own question with the example data you offered. The geom_raster and geom_tile functions have limits to their use. A big one is that because they require the data to be on an even grid, with very high res station data the machine will try to create an even grid before plotting, hence the error of needing 39.7 GB of memory. By rounding your lon/lat values to an even 0.1x0.1 grid the machine is now able to plot the data in the other answers in this post. Try running the following code chunk on your data.

mapoc_temp <- mapoc_temp %>% 
  ungroup() %>% 
  mutate(Longitude = plyr::round_any(Longitude, 0.1), 
         Latitude = plyr::round_any(Latitude, 0.1))

You may change the level of rounding as you prefer, but avoid going too fine scale unless it is absolutely necessary. Once your data have been rounded to an even grid you can plot them with the following chunk.

ggplot(mapoc_temp, aes(x = Longitude, y = Latitude)) +
  borders(fill = "grey80") +
  geom_raster(aes(fill = Temp)) +
  coord_quickmap(xlim = c(-64.5,-62.8), ylim = c(42.7,45)) +
  facet_wrap(~month)

There is no need to create any sf shape polygons or add any CRS information to the plot. The tidyverse is not designed to utilise this information and in my opinion it only complicates things unnecessarily to bring in the sf package. I hope this helps!

Answer 2

It looks like the geom_polygon(mapping layer) was overplotting the geom_tile layer. Below is an example using your example data with geom_tile called last, and with alpha set to 0.4 to show the map beneath.

#In my dataframe the lat and long values are not equally spaced!
  mapoc_temp = expand.grid(data.frame(Longitude= seq(-64.5,-62.4,0.1), 
                                      Latitude= seq(42.7,44.8,0.1),
                                      year = sample(c(2016,2017,2018), 22, replace = T),
                                      month = sample(month.abb, 22, replace = T)))
  mapoc_temp$Temp = runif(nrow(mapoc_temp))

  library(mapdata)
#> Loading required package: maps
  library(ggplot2)

  #map I'm using for ggplot
  canada = map_data("worldHires", "Canada")

  ## Use geom_raster as the last layer in the ggplot2 call,
  ##  otherwise the polygons plot over the tiles.
  ## Below alpha is set on the raster layer to show underlying map.
  ggplot(mapoc_temp, aes(x=Longitude, y=Latitude)) + 
    #try to map my temp on ocean
    geom_polygon(data = canada, aes(x=long, y=lat, group=group), colour="grey50", fill = 'grey55')+
    geom_raster(aes(fill = Temp, x = Longitude),alpha = .4, interpolate = TRUE) +
    coord_sf(xlim=c(-64.5,-62.8), ylim=c(42.7,45)) +
    #get months
    facet_wrap(vars(month))

^{Created on 2020-02-22 by the reprex package (v0.3.0)}