Projecting a MODIS SST netcdf file in R

Question

I've downloaded monthly SST data (11u, nighttime) for July 2014.

https://oceancolor.gsfc.nasa.gov/cgi/l3?per=MO&prd=NSST_sst.nc&sen=A&res=4km&num=24&ctg=Standard&date=1Jul2014

The data is stored as a NetCDF (.nc) file. So, I've used the following code to import the file into R.

library(ncdf4)
library(raster)
sst<-raster(path.expand("~/Remote Sensing/SST/SST_July2014_L3_MO_4km.nc"))

Next, I checked the file.

sst 

This was the output:

class       : RasterLayer 
dimensions  : 4320, 8640, 37324800  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -90.00001, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : ~\Remote Sensing\SST\SST_25June2016_L3_8D_4km.nc 
names       : Sea.Surface.Temperature 
zvar        : sst 

Does the raster() function in ncdf4 automatically project the file? I'm confused because the MODIS SST metadata states that the projection is Equidistant Cylindrical. Sorry if this question is very elementary, but this is my first time working with rasters or MODIS data.

EDIT

I successfully specified the projection using the code provided by @Phil

sst<-raster(path.expand("~/Remote Sensing/SST/SST_July2014_L3_MO_4km.nc"))
sst@crs <- sp::CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs")

My other data files are in latlong WGS84, so I reprojected the SST raster.

sst<-projectRaster(sst, crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

I then tried to crop the bathymetry raster by the extent of my polygon layer called 'polygonsp.' (code for polygonsp below)

library(sp)
vertices<-c(-131.8246,-131.8246,-127.6479,-127.6479,50.14335,54.53591,54.53591,50.14335)
vertices<-matrix(vertices,4,2)
polygon<-Polygon(vertices)
polygono<-Polygons(list(polygon),1)
polygonsp<-SpatialPolygons(list(polygono))
proj4string(polygonsp)<-CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

sstcrop<-crop(sst,polygonsp)

However, I now get an error that says "Error in .local(x, y, ...) : extents do not overlap".

This code worked for me before I added the line specifying the projection of the SST data as equidistant cylindrical. Any idea what is happening?

Answer 1

It doesn't look like there's an associated projection with the file, and the default behaviour of raster::raster() when the CRS is missing is (taken from the function documentation, ?raster):

If this argument is missing, and the x coordinates are withing [sic] -360 .. 360 and the y coordinates are within -90 .. 90, "+proj=longlat +datum=WGS84" is used.

So it looks like WGS84 has been assigned for you. To change it, simply use sp::CRS() and assign it to the @crs slot:

sst@crs <- sp::CRS("+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371007 +b=6371007 +units=m +no_defs")

CRS definition taken from spatial reference.

And the result:

rasterVis::levelplot(sst)

Answer 2

raster() function use image metadata to set projection. By metadata, SST product should came projected on Equidistant Cylindrical:

library(ncdf4)

nc_open("~/Downloads/A20151822015212.L3m_MO_NSST_sst_4km.nc")

File ~/Downloads/A20151822015212.L3m_MO_NSST_sst_4km.nc (NC_FORMAT_NETCDF4):

     3 variables (excluding dimension variables):
        unsigned short sst[lon,lat]   (Chunking: [64,64])  (Compression: level 4)
(...)
    65 global attributes:
        product_name: A20151822015212.L3m_MO_NSST_sst_4km.nc
        instrument: MODIS
        title: MODIS Level-3 Standard Mapped Image
        (...)
        map_projection: Equidistant Cylindrical
        latitude_units: degrees_north
        longitude_units: degrees_east
        northernmost_latitude: 90
        southernmost_latitude: -90
        westernmost_longitude: -180
        easternmost_longitude: 180
        geospatial_lat_max: 90
        geospatial_lat_min: -90
        geospatial_lon_max: 180
        geospatial_lon_min: -180
        grid_mapping_name: latitude_longitude
        latitude_step: 0.0416666679084301
        longitude_step: 0.0416666679084301
        sw_point_latitude: -89.9791641235352
        sw_point_longitude: -179.97917175293
        geospatial_lon_resolution: 4.59999990463257
        geospatial_lat_resolution: 4.59999990463257
        geospatial_lat_units: km
        geospatial_lon_units: km
        spatialResolution: 4.60 km
(...)

The weird thing is raster() open the image with lat-lon projection:

library(raster)

r2 <- raster("~/Downloads/A20151822015212.L3m_MO_NSST_sst_4km.nc",var="sst")

proj4string(r2)
[1] "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

But, you can create raster from .nc file with correct projection:

library(ncdf4)

library(raster)

nc.r <- nc_open("~/Downloads/A20151822015212.L3m_MO_NSST_sst_4km.nc")

sst <- ncvar_get(nc = nc.r,varid = "sst")

proj <- "+proj=eqc +lat_ts=60 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

r <- raster(nrows=dim(sst)[2], ncols=dim(sst)[1], xmn=-180, xmx=180, ymn=-90, ymx=90, crs = proj)

r <- setValues(r, values=as.vector(sst))

r
class       : RasterLayer 
dimensions  : 4320, 8640, 37324800  (nrow, ncol, ncell)
resolution  : 0.04166667, 0.04166667  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=eqc +lat_ts=60 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
data source : in memory
names       : layer 
values      : -2, 45.00072  (min, max)

It has the same boundary box, resolution and projection than metadata information. Both rasters, are really different?

r2 <- raster("~/Downloads/A20151822015212.L3m_MO_NSST_sst_4km.nc", var="sst")

test <- as.matrix(r) - as.matrix(r2)
max(test)
[1] 0

min(test)
[1] 0

plot(r, alpha=0.5)
plot(r2,add=T, alpha=0.5)

And, checking raster in QGIS I got this. You should not have problems using raster() or creating the raster from values.