How to correctly add a map to a bitmap in R

Question

I am trying to build sea surface temperature data and add a color image of the earth so that the data is not confused with NAs. I tried several methods for this, but, as you will see in the figures below, the maps do not align correctly with the data.

To make this problem reproducible, here is the link to dropbox with the file I'm working with: https://www.dropbox.com/s/e8pwgmnhvw4s0nf/sst.nc4?dl=0 p>

Here is the code that I have developed so far;

library(ncdf4)
library(raster)
library(mapdata)
library(mapproj)
library(rgeos)
library(ggplot2)

Via ncdf4, rasterToPoints, map_data and ggplot2

eight = nc_open("Downloads/sst.nc4")
sst = ncvar_get(eight, "sst")
sst = raster(sst)
sst = t(flip(sst, 1)) # have to orient the data properly

# extract the dimensions and set the extent
lat.min = min(eight$dim$lat$vals)
lat.max = max(eight$dim$lat$vals)
lon.min = min(eight$dim$lon$vals)
lon.max = max(eight$dim$lon$vals)
sst = setExtent(sst, ext = c(lon.min, lon.max, lat.min, lat.max))

# provide proper projection
crs(sst) = "+init=epsg:4326"

# convert raster to points
sst.p <- rasterToPoints(sst)
df <- data.frame(sst.p)
colnames(df) <- c("Longitude", "Latitude", "sst")
usa = map_data("usa")
ggplot(data=df, aes(y=Latitude, x=Longitude)) +
  geom_raster(aes(fill=sst)) +
  theme_bw() +
  coord_equal() +
  scale_fill_gradient("SST (Celsius)", limits=c(0,35)) +
  geom_polygon(data = usa, aes(x=long, y = lat, group = group)) + 

  theme(axis.title.x = element_text(size=16),
        axis.title.y = element_text(size=16, angle=90),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        legend.position = "right",
        legend.key = element_blank()
  )

Via Raster, maptools data, SP Transform and basic build

#read in the data
sst = raster("Downloads/sst.nc4",  varname = "sst", stopIfNotEqualSpaced=FALSE)

# get world map data
data("wrld_simpl", package="maptools")

## Crop to the desired extent, then plot
newext <- c(lon.min, lon.max, lat.min, lat.max)
out <- crop(wrld_simpl, newext)

#transform to proper CRS
out = spTransform(out, "+init=epsg:4326")

#plot
plot(out, col="khaki", bg="azure2")
plot(sst, add = T)

- The projection that I use for this spatial data, EPSG:4326

-This is an XML fragment defining output output sst.nc4

<crs>PROJCS["Mercator_1SP / World Geodetic System 1984",
         GEOGCS["World Geodetic System 1984",
         DATUM["World Geodetic System 1984",
         SPHEROID["WGS 84", 6378135.0, 298.257223563, AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
         UNIT["degree", 0.017453292519943295],
         AXIS["Geodetic longitude", EAST],
         AXIS["Geodetic latitude", NORTH]],
         PROJECTION["Mercator_1SP"],
         PARAMETER["latitude_of_origin", 0.0],
         PARAMETER["central_meridian", 0.0],
         PARAMETER["scale_factor", 1.0],
         PARAMETER["false_easting", 0.0],
         PARAMETER["false_northing", 0.0],
         UNIT["m", 1.0],
         AXIS["Easting", EAST],
         AXIS["Northing", NORTH]]</crs>

map() mapproj projection, , , - .

Answer 1

. -

sst = raster("sst.nc4",  varname = "sst")

:

"cells are not equally spaced; you should extract values as points"

, :

library(ncdf4)
library(raster)
library(maptools)

d <- nc_open("sst.nc4")
sst <- ncvar_get(d, "sst")
lon <- ncvar_get(d, "lon")
lat <- ncvar_get(d, "lat")
nc_close(d)

xy <- cbind(rep(lon, length(lat)), rep(lat, each=length(lon)))

NA ( ...

xyv <- na.omit(cbind(xy, as.vector(sst)))

RasterLayer , ,

 r <- raster(extent(range(lon), range(lat)), res=1/6)
 r <- rasterize(xyv[, 1:2], r, xyv[,3], fun=mean) 

Plot

data(wrld_simpl)
w <- crop(wrld_simpl, r)

plot(r)
plot(w, col='gray', add=TRUE)