Skip to content

Instantly share code, notes, and snippets.

@aolinto

aolinto/COPERNICUS_Multiobs_3d.R

Created August 3, 2019 19:12
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save aolinto/9d9514eb453109425f5b35c42fa0eb95 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save aolinto/9d9514eb453109425f5b35c42fa0eb95 to your computer and use it in GitHub Desktop.
Download ZIP
R script to extract COPERNICUS Temperature and Salinity statistics from a nc file based on a polygon
Raw
COPERNICUS_Multiobs_3d.R
# Antonio Olinto Avila-da-Silva, Instituto de Pesca, Brasil
# https://gist.github.com/aolinto/9d9514eb453109425f5b35c42fa0eb95
# 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-armor-3d-rep-monthly
# the script will transform nc file to brick raster, read temperature and salinity data
# for a given area and depth, compute its statistics and write them into
# a single csv file named COPERNICUS_to.csv or COPERNICUS_so.csv
# version 2019/08/03
# nc help http://geog.uoregon.edu/GeogR/topics/netCDF-read-ncdf4.html
# load libraries
# --------------
# install.packages(c("ncdf4","rgdal","maptools","raster"))
# ncdf4 needs libnetcdf-dev netcdf-bin in Linux
library(rgdal)
library(maptools)
library(raster)
library(ncdf4)
# set workspace
# --------------
# list and remove objects from workspace
ls()
rm(list = ls())
# set working directory
#~ setwd("/mnt/dados/Sensoriamento_Remoto/Multiobs") # indicate the path to the nc file
setwd("/mnt/Guilherme/Antonio01/Sensoriamento_Remoto/Multiobs")
# verify the existence of COPERNICUS_to-so.csv file
file.exists("COPERNICUS_to-so.csv") # caution! new data will be appended to this file if it already exists
# if TRUE choose an option
#~ file.rename("COPERNICUS_to-so.csv","COPERNICUS_to-so.old")
#~ file.remove("COPERNICUS_to-so.csv")
# load shapefiles
# ---------------
# study area
#~ shp.area <- readOGR("/home/antonio/Documentos/Estudos/Polvo/SIG","AreaPolvoPPW") # indicate the location and name of the shapefile
#~ shp.area <- readOGR("/mnt/Guilherme/Antonio01/Documentos/Estudos/Polvo/SIG","AreaPolvoPPW_50-75") # indicate the location and name of the shapefile
shp.area <- readOGR("/mnt/Guilherme/Antonio01/Documentos/Estudos/Polvo/SIG","AreaPolvoPPW_75-100") # indicate the location and name 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)
# check ncfile information
# ------------------------
nc.data <- nc_open("armor-3d-monthly_2003-2018_0-100.nc")
nc.data
#. list of the nc variable names
var.names <- attributes(nc.data$var)$names
var.names
#. variables settings
ncatt_get(nc.data,var.names[1])
ncatt_get(nc.data,var.names[2])
#. list of the nc dimensions names
attributes(nc.data$dim)$names
#. dimensions info
ncatt_get(nc.data,"time")
ncatt_get(nc.data,"time","units")
(time_min <- ncatt_get(nc.data,"time","valid_min")$value)
as.Date(time_min/24, origin="1950-01-01")
(time_max <- ncatt_get(nc.data,"time","valid_max")$value)
as.Date(time_max/24, origin="1950-01-01")
ncvar_get(nc.data,"time")
as.Date(ncvar_get(nc.data,"time")/24, origin="1950-01-01")
ncatt_get(nc.data,"depth")
ncvar_get(nc.data,"depth") # [1] 0 10 20 30 50 75 100
ncatt_get(nc.data,"latitude")
ncvar_get(nc.data,"latitude")
ncatt_get(nc.data,"longitude")
ncvar_get(nc.data,"longitude")
#. close
nc_close(nc.data)
# ===========================
# loopings for data extraction
# ===========================
dl <- 6 # select depth level number: 5 for 50-75, 6 for 75-100 in this data set
dls <- ifelse(dl==5,"50-75","75-100")
for (br in 1:length(var.names)) { # looping variables
# brickraster
# -----------
# create a brickraster and transform long axis
rst.data <- brick("armor-3d-monthly_2003-2018_0-100.nc",varname=var.names[br],level=dl) # indicate the nc file, varname (to, so) and level (depth layer 1 to 7: 0, 10, ..., 100)
rst.data <- rotate(rst.data) # 0,360 to -180,180 longitude
# exploratory summaries and plots
# -------------------------------
# explore stack
#~ rst.data
#~ names(rst.data)[1:48]
#~ getZ(rst.data)[1:48] # dates in original format
#~ as.Date(getZ(rst.data)[1:48]/24, origin="1950-01-01") # converted dates
#. summary by month
#~ summary(rst.data[[1:18]])
#~ dat.min <- summary(rst.data[[1:179]])[1,]
#~ dat.median <- summary(rst.data[[1:179]])[3,]
#~ dat.max <- summary(rst.data[[1:179]])[5,]
#~ plot(rst.data,zlim=c(min(dat.min),max(dat.max)))
#~ plot(dat.median,type="l",col="orange",ylim=c(min(dat.min),max(dat.max)))
#~ points(dat.min,type="l",col="darkgreen")
#~ points(dat.max,type="l",col="blue")
#~ abline(v=seq(0,179,12),col="gray",lty=2)
#. plot a sample layer, study area and land
#~ plot(rst.data[[1]],main=names(rst.data)[1],zlim=c(min(dat.min),max(dat.max)))
#~ plot(shp.area,border="red",add=T)
#~ plot(shp.land0,border="darkgreen",col="gray",add=T)
#. plot variable by month/year
#~ dev.new(width=15, height=10, unit="in")
#~ par(mfrow=c(3,4),mar=c(3,3,3,5))
#~ for (i in seq(2,36,3)) {
#~ plot(rst.data[[i]],main=as.Date(getZ(rst.data)[i]/24, origin="1950-01-01"),zlim=c(min(dat.min),max(dat.max)))
#~ plot(shp.area,border="red",add=T)
#~ plot(shp.land0,border="darkgreen",col="gray",add=T)
#~ }
# crop data from the area
# -----------------------
# 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)
# view sample data
#~ dev.new(width=15, height=10, unit="in")
#~ par(mfrow=c(3,4),mar=c(3,3,3,5))
#~ for (i in seq(2,36,3)) {
#~ plot(msk.data[[i]],main=as.Date(getZ(rst.data)[i]/24, origin="1950-01-01"),zlim=c(min(dat.min),max(dat.max)))
#~ plot(shp.area,border="red",add=T)
#~ plot(shp.land0,border="darkgreen",col="gray",add=T)
#~ }
# looping for data statistical calculation
for (i in seq(1,nlayers(rst.data))) { # looping year/month
# get data
year <- as.numeric(substr(as.Date(getZ(rst.data)[i]/24, origin="1950-01-01"),1,4))
month <- as.numeric(substr(as.Date(getZ(rst.data)[i]/24, origin="1950-01-01"),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(var.names[br],dls,year,month,sta.min,sta.mean,sta.median,sta.max)
names(dat.output)<-c("variable","depth","year","month","min","mean","median","max")
# save csv file
fe <- file.exists("COPERNICUS_to-so.csv")
write.table(dat.output,"COPERNICUS_to-so.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)
} # finish looping year/month
rm(rst.data,ext.area,crp.data,crp.na,rst.mask,msk.data)
} # finish looping variables
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment