R script to extract Aqua MODIS SST statistics from a nc file based on a polygon

MODISA_sst_pol.R

# Antonio Olinto Avila-da-Silva, Instituto de Pesca, Brasil
# https://gist.github.com/aolinto/3a0872e0fb61b7ca3e69d35d286ae26c
# script to process Aqua MODIS Sea Surface Temperature
# files downloaded from https://oceancolor.gsfc.nasa.gov/cgi/l3
# Aqua MODIS Sea Surface temperature 11 µ daytime Monthly 4 km SMI images or
# Aqua MODIS Sea Surface temperature 11 µ nigthtime Monthly 4 km SMI images
# all .L3m_MO_SST_sst_4km.nc (daytime) or .L3m_MO_NST4_sst4_4km (nigthtime)
# files must be in the working directory
# the script will open each nc file to read date information
# the script will also transform nc file to raster, read sst data
# for a given area indicated by a shapefile, compute its statistics
# and write them into a single csv file named MODISA_sst.csv
# Some reference pages
# http://geog.uoregon.edu/GeogR/topics/netCDF-read-ncdf4.html
# https://scottishsnow.wordpress.com/2014/08/24/many-rastered-beast/
# version 2019/02/11

# run in the terminal
# Aqua MODIS Sea surface temperature 11 µ nighttime Monthly 4 km
# wget -q --post-data="sensor=aqua&sdate=2008-01-01&edate=2018-12-31&dtype=L3m&addurl=1&results_as_file=1&search=A*L3m_MO_NSST_sst_4km.nc" -O - https://oceandata.sci.gsfc.nasa.gov/api/file_search |wget -i -
# Aqua MODIS Sea surface temperature 11 µ daytime Monthly 4 km
# wget -q --post-data="sensor=aqua&sdate=2008-01-01&edate=2018-12-31&dtype=L3m&addurl=1&results_as_file=1&search=A*L3m_MO_SST_sst_4km.nc" -O - https://oceandata.sci.gsfc.nasa.gov/api/file_search |wget -i -

# load libraries
# ncdf4 needs libnetcdf-dev netcdf-bin in Linux
# install.packages(c("rgeos","maptools","ncdf4","raster","rgdal"))
library(maptools)
library(ncdf4)
library(raster)
library(rgdal)

# set working directory
setwd("/mnt/...")   # indicate the path to nc files

# verify the existence of MODISA_nsst.csv file
file.exists("MODISA_nsst.csv") # caution! new data will be appended to this file if it already exists
# if TRUE choose an option
#~ file.rename("MODISA_nsst.csv","MODISA_nsst.org")
#~ file.remove("MODISA_nsst.csv")

# list and remove objects from workspace
ls()
rm(list = ls())

# create a list of nc files and indicate its length
(f <- list.files(".", pattern="*.L3m_MO_NSST_sst_4km.nc",full.names=F)) # for daytime *.L3m_MO_SST_sst_4km.nc
(lf<-length(f))

# load shapefile
shp.area <- readOGR("/home/...","shapefile name") # indicate the location and name of the shapefile
shp.area <- spTransform(shp.area, CRS("+proj=longlat +datum=WGS84"))
proj4string(shp.area)
summary(shp.area)

# plot shapefile
plot(shp.area)

# get the spatial extent of the shapefile
ext.area <- extent(shp.area)

for (i in 1:lf) {
  # progress indicator
  print(sprintf("Processing file %i from %s: %s",i,length(f),f[i]))
  
  # open netCDF file
  nc.data <- nc_open(f[i])
  
  # extract date information from nc file
  dateini <- ncatt_get(nc.data,0,"time_coverage_start")$value
  dateend <- ncatt_get(nc.data,0,"time_coverage_end")$value
  datemean <- mean(c(as.Date(dateend,"%Y-%m-%dT%H:%M:%OSZ"),as.Date(dateini,"%Y-%m-%dT%H:%M:%OSZ")))
  year <- substring(datemean,0,4)
  month <- substring(datemean,6,7)
  
  # close netCDF file
  nc_close(nc.data)
  
  # create a raster from nc file
  rst.data <- raster(f[i],varname="sst")
  proj4string(rst.data) <- CRS("+proj=longlat +datum=WGS84")
    
  # crop the raster to area extent
  crp.data <- crop(rst.data,ext.area,snap="out")
  
  # set values higher than 45 to NA
  crp.data[crp.data>=45] <- NA
  
  # create a dummy raster with NAs
  crp.na <- setValues(crp.data,NA)

  # create a raster mask with the area boundaries
  rst.mask <- rasterize(shp.area,crp.na)

  # apply the mask to the raster with data
  msk.data <- mask(x=crp.data,mask=rst.mask)
  
  # get statistics
  sta.min <- cellStats(msk.data, stat='min',na.rm=TRUE)
  sta.mean <- cellStats(msk.data, stat='mean',na.rm=TRUE)
  sta.median <- median(na.omit(values(msk.data)))
  sta.max <- cellStats(msk.data, stat='max',na.rm=TRUE)
  
  # prepare final data set
  dat.output <- data.frame(year,month,sta.min,sta.mean,sta.median,sta.max)
  names(dat.output)<-c("year","month","SSTmin","SSTmean","SSTmedian","SSTmax")
 
  # save csv file
  fe <- file.exists("MODISA_nsst.csv")
  write.table(dat.output,"MODISA_nsst.csv",row.names=FALSE,col.names=!fe,sep=",",dec=".",append=fe) # change separator and decimal strings if necessary 
  
  # clean workspace
  rm(nc.data,dateini,dateend,datemean,year,month,rst.data,crp.data,crp.na,rst.mask,msk.data,sta.min,sta.mean,sta.median,sta.max,dat.output,fe)
}

rm(ext.area,f,i,lf,shp.area)