R script to extract COPERNICUS sea surface salinity data from a polygon

COPERNICUS_sss_pol.R

# Antonio Olinto Avila-da-Silva, Instituto de Pesca, Brasil
# https://gist.github.com/aolinto/59e43cdbac3ff3e0cd5eb2eab5efcbc4
# script to process COPERNICUS Sea Surface Salinity
# product MULTIOBS_GLO_PHY_REP_015_002
# nc file downloaded from
# http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=MULTIOBS_GLO_PHY_REP_015_002
# dataset-sss-ssd-rep-monthly sos variable
# the script will transform nc file to stack raster, read sss data
# for a given area, compute its statistics and write them into
# a single csv file named COPERNICUS_sss.csv
# version 2018/12/28

# load libraries
# install.packages(c("rgdal","maptools","raster"))
library(rgdal)
library(maptools)
library(raster)

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

# set working directory
setwd("...")   # indicate the path to the nc file

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

# land area
shp.land0 <- getData('GADM', country='BRA', level=0)
shp.land0 <- spTransform(shp.land0, CRS("+proj=longlat +datum=WGS84"))
proj4string(shp.land0)

shp.land2 <- getData('GADM', country='BRA', level=2)
shp.land2 <- spTransform(shp.land2, CRS("+proj=longlat +datum=WGS84"))
proj4string(shp.land2)

# plot study area and land
plot(shp.area,col="lightblue",axes=T)
plot(shp.land2,border="darkgreen",col="gray",add=T)

# create a stackraster and transform long axis
rst.data <- stack("file.nc",varname="sos") # indicate the nc file
rst.data <- rotate(rst.data) # 0,360 to -180,180 longitude

# explore stack
rst.data
summary(rst.data[[1:48]])
names(rst.data)[1:48]
getZ(rst.data)[1:48]
plot(rst.data, zlim=c(34,38))

# plot a sample layer, study area and land
plot(rst.data[[1]],main=names(rst.data)[1],xlim=c(-47.5,-45),ylim=c(-25.5,-23.5),zlim=c(34,38))
plot(shp.area,border="red",add=T)
plot(shp.land2,border="darkgreen",col="gray",add=T)

# plot salinity by month/year
dev.new(width=10, height=10, unit="in")
par(mfrow=c(4,3))
for (i in seq(25,48,2)) {
plot(rst.data[[i]],main=substr(getZ(rst.data)[i],1,7),xlim=c(-47.5,-45),ylim=c(-25.5,-23.5),zlim=c(34,38))
plot(shp.area,border="red",add=T)
plot(shp.land0,border="darkgreen",col="gray",add=T)
}

# get the spatial extent of the shapefile
ext.area <- extent(shp.area)
    
# crop the raster to area extent
crp.data <- crop(rst.data,ext.area,snap="out")
  
# create a dummy raster with NAs
crp.na <- setValues(crp.data[[1]],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)

dev.new(width=10, height=10, unit="in")
par(mfrow=c(4,3),mar=c(3,3,3,5))
for (i in seq(25,48,2)) {
plot(msk.data[[i]],main=substr(getZ(rst.data)[i],1,7),xlim=c(-47.5,-45),ylim=c(-25.5,-23.5),zlim=c(34,38))
plot(shp.area,border="red",add=T)
plot(shp.land0,border="darkgreen",col="gray",add=T)
}

for (i in seq(1,nlayers(rst.data),2)) {

	# get data
	year <- as.numeric(substr(getZ(msk.data)[i],1,4))
	month <- as.numeric(substr(getZ(msk.data)[i],6,7))
	  
	# get statistics
	sta.min <- cellStats(msk.data[[i]], stat='min',na.rm=TRUE)
	sta.mean <- cellStats(msk.data[[i]], stat='mean',na.rm=TRUE)
	sta.median <- median(na.omit(values(msk.data[[i]])))
	sta.max <- cellStats(msk.data[[i]], 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","SSSmin","SSSmean","SSSmedian","SSSmax")

	# save csv file
	fe <- file.exists("COPERNICUS_sss.csv")
	write.table(dat.output,"COPERNICUS_sss.csv",row.names=FALSE,col.names=!fe,sep=",",dec=".",append=fe) # change separator and decimal strings if necessary   

	# clean workspace
	rm(year,month,sta.min,sta.mean,sta.median,sta.max,dat.output,fe) 
}