How to make a map showing temperature across an area

Question

I am hoping to get some advice/guidance on how to make a map as an R newbie.

I want to make a map showing a temperature gradient across a river. I have the longitudes, latitudes, temperature and other environmental variables such as salinity. I am envisioning something like this:

map

This is what my data looks like:

> dput(water)
structure(list(ID = c(482914L, 482555L, 482536L, 482928L, 482592L, 
482563L, 482520L, 482530L, 482942L, 482508L, 482573L, 482570L, 
482524L, 482505L, 482599L, 482544L, 482539L, 482533L, 482514L, 
482516L, 482933L, 482951L, 482909L, 482904L, 482578L, 482511L, 
482502L, 482549L, 482583L), Lat = c(49.1655, 48.2795, 48.057667, 
49.764, 48.954, 48.558667, 47.742333, 47.696, 49.616, 47.094667, 
48.575667, 48.473667, 47.696, 46.887167, 49.0365, 48.348333, 
48.072833, 48.027667, 47.403833, 47.3935, 49.749167, 49.412833, 
49.040333, 48.830667, 48.787833, 47.415667, 49.0365, 48.317333, 
48.666833), Long = c(-67.2685, -69.072167, -69.548833, -66.8895, 
-68.0945, -68.845333, -69.927833, -69.821667, -66.661167, -70.716167, 
-68.489833, -68.6655, -69.821667, -70.9235, -68.267667, -69.327667, 
-69.6245, -69.485667, -70.283833, -70.229833, -65.500333, -66.481, 
-67.104167, -67.857, -68.711167, -70.440833, -68.267667, -69.199667, 
-68.581167), Temp = c(10.48, 5, 6.77, 9.48, 10.72, 6.04, 8.33, 
10.05, 11.37, 19.26, 6.92, 6.4, 10, 23.7, 4.87, 5.95, 2.88, 5.49, 
17.55, 15.27, 9.54, 10.74, 6.32, 8.97, 9.84, 13.45, 23.4, 5.2, 
11.12), Salinity = c(29.16, 28.87, 26.78, 29.53, 24.77, 29.11, 
24.14, 22.81, 28.99, 6.82, 28.63, 28.78, 23.13, 0.12, 29.9, 29.42, 
30.77, 28.37, 9.96, 13.57, 29.97, 29.18, 28.99, 27.65, 27.2, 
16.31, 0.12, 29.41, 27.54)), class = "data.frame", row.names = c(NA, 
-29L))

I am not sure how to go about starting and which packages to use... I have gotten as far as to crop the area I am interested in using ggmap

library(ggmap)
qc <- c(left = -74, bottom = 45, right = -58, top = 52)
map <- get_stamenmap(qc, zoom = 7, maptype = "toner-lite")
ggmap(map)

Any help is appreciated, thank you.

Your observations in the example are relatively few points. Either you could add them as geom_points(), or you would need to think about some sort of spatial model to interpolate between them - maybe something like: https://stackoverflow.com/questions/41580004/how-to-plot-interpolating-data-on-a-projected-map-using-ggplot2-in-r — Miff, Jul 17 '23 at 14:44

Allan Cameron · Answer 1 · 2023-07-18T19:21:43.987

You only have a few points, so if you want something similar to the picture in your question, you need to interpolate these measurements over a 2D grid. You can do that with the interp package:

library(ggmap)
library(interp)

qc <- c(left = -74, bottom = 45, right = -58, top = 52)
map <- get_stamenmap(qc, zoom = 7, maptype = "toner-lite")

ggmap(map) +
  geom_tile(aes(x = x, y = y, fill = z) ,
             data = with(water, interp(Long, Lat, Temp, nx = 100, ny = 100,
                                       duplicate = "mean", extrap = TRUE)) |> 
               interp2xyz() |> as.data.frame()) +
  scale_fill_gradientn("Temperature", na.value = NA,
                       colours = c("blue", "green", "yellow", "orange"))

An easier, more honest, but less "meteorological" plot would be to draw your measurements as large colored points:

ggmap(map) +
  geom_point(data = water, aes(Long, Lat, colour = Temp), size = 5) +
  scale_color_gradientn("Temperature", na.value = NA,
                       colours = c("blue", "green", "yellow", "orange"))

Edit

If you want a terrain-type map, you can change the maptype

map <- get_stamenmap(qc, zoom = 7, maptype = "terrain")

Which gives you

If you want to plot the interpolation only on the water, you can strip out the water from the map object and make it transparent. You would then plot it as a geom_raster.

qc <- c(left = -74, bottom = 45, right = -58, top = 52)
map <- get_stamenmap(qc, zoom = 7, maptype = "terrain")

attrib <- attributes(map)
map <- c(map)
map[map %in% c('#99B3CC')] <- 'NA'

xvals <- seq(attrib$bb$ll.lon, attrib$bb$ur.lon, len = attrib$dim[2])
yvals <- seq(attrib$bb$ur.lat, attrib$bb$ll.lat, len = attrib$dim[1])
df <- cbind(expand.grid(x = xvals, y = yvals), z = map)

ggplot() +
  geom_raster(data = with(water, interp(Long, Lat, Temp, 
                                      nx = 100, ny = 100,
                                      duplicate = "mean")) |> 
              interp2xyz() |> as.data.frame(),
            aes(x, y, fill = z), alpha = 0.8) +
  scale_fill_gradientn("Temperature", na.value = NA,
                       colours = c("blue", "green", "yellow", "orange")) +
  ggnewscale::new_scale_fill() +
  geom_raster(data = df, aes(x, y, fill = z)) +
  scale_fill_identity() +
  coord_equal(ratio = 1.9, xlim = c(-74, -58), expand = FALSE) +
  theme(panel.background = element_rect(fill = '#99B3CC', color = NA),
        panel.grid = element_blank())

Thank you this helps a lot! I'm wondering if you have suggestions on where to acquire a better base map with topography. I am using ggmap but maybe a different package? Apologies if this is not the right forum for this question... — no_frills_30, Jul 17 '23 at 19:02
@no_frills_30 why not just change `maptype = "toner-lite"` to `maptype = "terrain"` in `get_stamenmap`? — Allan Cameron, Jul 17 '23 at 19:58
I was able to find a larger dataset with more coordinates. Would it be possible to have the temperature gradient following the true shape of the river? Instead of a polygon shape that currently overlaps with the land? And is it possible to keep the labels on the original map as the temperature gradient now overlaps. — no_frills_30, Jul 18 '23 at 11:25
@no_frills_30 it is possible to make the temperature map follow the river, but it would be considerably more complex, requiring a shapefile or raster to mask out the output of the interpolation (though it would be straightforward to obtain these). The labels that are blocked out by the temperature map would need to be drawn in again over the top. If you want to share your data I can have a look. — Allan Cameron, Jul 18 '23 at 12:28
Hi thank you. Here is a link to my data: https://docs.google.com/spreadsheets/d/e/2PACX-1vSSGCvqtwAuDK6kvwACJ1CACcTJO4ni_YnT1Uw7xx_iEpNW6jfiY5jHIbpFvIrdf7YwBJz9uRsfq9Dm/pubhtml?gid=850473751&single=true I have found a 'shapefile' of Canada here: https://www.diva-gis.org/datadown but want to focus in on a specific area, as outlined above. Is this all I need? — no_frills_30, Jul 18 '23 at 12:58
@no_frills_30 that data set isn't very useful for this purpose, since it follows a narrow band of measurements way out to sea, so doesn't fill the whole ocean area. I have just used your original data for illustration, but you can use the same code on other data you have to get a satisfactory result — Allan Cameron, Jul 18 '23 at 17:06
Thank you, I have tried the code. However, in geom_raster, what is 'df'? I changed 'df' to the original dataset of 'water' and I get this error: Error in `geom_raster()`: ! Problem while computing aesthetics. ℹ Error occurred in the 2nd layer. Caused by error in `FUN()`: ! object 'x' not found — no_frills_30, Jul 18 '23 at 18:44

How to make a map showing temperature across an area

1 Answers1