How to correctly add map to raster image in R

Question

I'm trying to plot sea surface temperature data and add a colored image of land so that data isn't confused with NAs. I've tried multiple methods to do so, but as you'll see in the images below, the maps do not line up properly relative to the data.

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

Here is the code I've developed thus 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 base plotting

#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 I'm using for this spatial data is EPSG:4326

-Here is the XML snippet dictating the sst.nc4 output projection

<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>

I've also attempted to use the map() function with mapproj's projection argument but it doesn't seem to have a pseudo-mercator projection as an option.

Patriciate

Can you get it properly in arcmap or Qgis? It seems to me that you have projection problem. Are you getting the data from AQUA/MODIS? – Romanticize

@Masoud I haven't tried arcmap or Qgis. Do you know of a tutorial for those packages? This is AQUA/MODIS data. – Patriciate

They are not R packages. They are gis software. NASA has a software called panoply (or a similar name). Install that and load your netcdf file into it. If that worked, then you need to change projection of your basemap in R to get the desired output. – Romanticize

Answer 1

This one is a bit confusing. The generally easiest approach would be

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

but, for this file, that gives this error:

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

So let's do that:

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)))

Combine and remove NA values (about half the cells...

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

Set up a RasterLayer with a resolution that is sufficient for your purposes, and rasterize the points

 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)

Trouvaille

this works! Will this method still work without omitting the NAs? I'd like to keep the NAs in if possible. My initial efforts to edit this were unsuccessful – Patriciate

You can keep the NAs but that won't change the results. – Trouvaille

is that due to the rasterizing that seems to have some sort of mean gapfilling functionality? – Patriciate