4

I have a data set with global observations with lat/lon information. There are 288 longitudes from -179.375 to 179.375 with spacing 1.25, 180 latitudes from -89.5 to 89.5 with spacing 1. There is one observation at each lat/lon combination. The total number of observations is 288*180=51840.

I want to make a heat map with world map embedded as well. Here is the code to make a rectangular map.

library(fields)
library(maps)
library(maptools)

set.seed(1)
unique.lat <- seq(-89.5, 89.5, by = 1)
unique.lon <- seq(-179.375, 179.375, by = 1.25)
lat <- rep(unique.lat, each = length(unique.lon))
lon <- rep(unique.lon, length(unique.lat))
y <- lat + abs(lon) + rnorm(length(lat))  # generate some data for plotting

dat <- data.frame(lat = lat, lon = lon, y = y)
data.mat <- matrix(dat$y, nrow = length(unique.lon), ncol = length(unique.lat), byrow = F)
surface <- list(x = unique.lon, y = unique.lat, z = data.mat)
plot.surface(surface, type = "I", axes = FALSE, xlab = "Longitude", ylab = "Latitude", legend.cex = 0.8, legend.mar = 7)
map("world", add = T)
map.axes()

enter image description here

What I really want is a heat map using mollweide projection. I learned how to make a simple world map with

data(wrld_simpl)
plot(spTransform(wrld_simpl, CRS = ("+proj=moll +lon_0=0 +x_0=0 +y_0=0")))

enter image description here

How do I somehow combine these two plots and make a heat map on the mollweide-projected map?

2

Here's a way, use the raster package. (If you really have data for this it's probably not a great idea to grid in longlat and then resample to Mollweide - but depends on your purpose).

R library(raster) library(maptools) r <- raster(surface) data(wrld_simpl) prj <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84" w <- spTransform(wrld_simpl, CRS = prj) p <- projectRaster(r, crs = prj) plot(p); plot(w, add = TRUE)

Please note that projectRaster gives warnings about infinite values, if you project the centre points from the original grid you don't see any out of range, so it may be raster doing extra sampling at the margins - outside of these centre points, I don't know exactly. sp/rgdal never lets you explore this outside of using project(), but clearly raster swallows the out of bounds values and works to return a result.

R a <- project(as.matrix(expand.grid(surface$x, surface$y)), "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84") sum(!is.finite(a)) [1] 0

  • Thanks. When I used p <- projectRaster(r, crs = prj), I got the following warning message: Warning message: In rgdal::rawTransform(projto_int, projfrom, nrow(xy), xy[, 1], : 13369 projected point(s) not finite. What does it mean? Should I worry about it? – Patrick Li May 30 '16 at 10:13
  • I added a note, whether you should worry depends on your data and your purpose, (I'm more worried about warping from longlat to Mollweide generally since you may be losing information unnecessarily). – mdsumner May 30 '16 at 14:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.