3

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?

  • 1
    The extents do not overlap error is because you're assigning a projection with units in meters to a dataset that clearly has units in degrees. So your sst object just represents a tiny little area somewhere in the gulf of Guinea :) – Loïc Dutrieux Jan 6 '17 at 0:02
  • Equidistant Cylindrical is like longlat, with units in meters. Because this dataset is in degrees, the assumption of raster to assign +proj=longlat to it was correct. If you open it in qgis, and define its projection to longlat, you'll see that it aligns perfectly with coastlines. – Loïc Dutrieux Jan 6 '17 at 11:32
  • @LoïcDutrieux thanks so much! is raster() also correct in assigning a datum of WGS84? – Splash1199 Jan 6 '17 at 16:24
  • Yes, when not specified it is usually correct to assume a WGS84 datum – Loïc Dutrieux Jan 9 '17 at 10:08
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)

enter image description here

  • Thanks! I'm now encountering a different problem, I've edited my question to explain – Splash1199 Jan 5 '17 at 19:38
  • @Splash1199 it may be that the image isn't the projection you thought. Either way, I'd ask a new question and link back to this. – Phil Jan 6 '17 at 12:34
1

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)

plot

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

plot2

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.