Comparison between packages RnetCDF, ncdf4 and ncdf

Let use the same nc file obtained from NOAA (http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html) be used.

Here I show the 3 ways to get annual data depending on the package used.

In conclusion, RnetCDF does not scale data using scale_factor and add_offset automatically. Take care.
utcal.nc from package RnetCDF is very useful.

dest <- "sst.day.mean.2014.v2.nc"

For all the package, I will use the same coordinates: lat=35, lon=28

# Note that the function ind_long_lat can be used data from any of these packages
library(RNetCDF)
dat.RNetCDF <- open.nc(dest)
lat=35
lon=28
# it will be used for all data
library(phenology)
indice <- ind_long_lat(ncdf=dat.RNetCDF, long=lon, lat=lat)

Packages RnetCDF
# take care, it does not scale automatically the data
point_SST.RNetCDF_wrong <- var.get.nc(ncfile=dat.RNetCDF, variable="sst", 
                                start=c(indice ,1), count=c(1, 1, 1))
# correct method
point_SST.RNetCDF_corr <- var.get.nc(ncfile=dat.RNetCDF, variable="sst", 
                                     start=c(indice ,1), count=c(1, 1, 1))*
  att.get.nc(ncfile=dat.RNetCDF, variable="sst", attribute="scale_factor")+
  att.get.nc(ncfile=dat.RNetCDF, variable="sst", attribute="add_offset")
# Get the time
label.time <- dim.inq.nc(dat.RNetCDF, 0)$name
date.char <- utcal.nc(att.get.nc(dat.RNetCDF, label.time, "units"), 
                            var.get.nc(dat.RNetCDF, label.time, 
                                       start=1, 
                                       count=dim.inq.nc(dat.RNetCDF, 0)$length), 
                            type="s")
date.POSIXlt <- strptime(date.char, format="%Y-%m-%d %H:%M:%S", tz="UTC")

Package ncdf4
library(ncdf4)
dat.ncdf4 <- nc_open(dest)
point_SST.ncdf4 <- ncvar_get(dat.ncdf4, varid="sst", start=c(indice ,1),count=c(1, 1, 1))
# Get the time
# We must check what is the number of time element using names(dat.ncdf$dim); here it is the first one
label.time <- names(dat.ncdf4$dim)[1]
library(RNetCDF)
date.char <- utcal.nc(dat.ncdf4$dim[[label.time]]$units, dat.ncdf4$dim[[label.time]]$vals, type="s")
date.POSIXlt <- strptime(date.char, format="%Y-%m-%d %H:%M:%S", tz="UTC")

Package ncdf
library(ncdf)
dat.ncdf <- open.ncdf(dest)
point_SST.ncdf <- get.var.ncdf(dat.ncdf, varid="sst", start=c(indice ,1),count=c(1, 1, 1))
# Get the time
# We must check what is the number of time element using names(dat.ncdf$dim); here it is the first one
label.time <- names(dat.ncdf$dim)[1]
library(RNetCDF)
date.char <- utcal.nc(dat.ncdf$dim[[label.time]]$units, dat.ncdf$dim[[label.time]]$vals, type="s")
date.POSIXlt <- strptime(date.char, format="%Y-%m-%d %H:%M:%S", tz="UTC")

# Alternative without utcal() function from RNetCDF package:
ref.date <- strptime(gsub("days since ", "", dat.ncdf$dim[[label.time]]$units), format="%Y-%m-%d %H:%M:%S", tz="UTC")
if (substr(dat.ncdf$dim[[label.time]]$units, 1, 4)=="days") 
  {date.POSIXlt <- ref.date+dat.ncdf$dim[[label.time]]$vals*24*3600} else 
    {if (substr(dat.ncdf$dim[[label.time]]$units, 1, 5)=="hours") 
  {date.POSIXlt <- ref.date+dat.ncdf$dim[[label.time]]$vals*3600} else 
  {date.POSIXlt <- rep(NA, length(dat.ncdf$dim[[label.time]]$vals))}}


Commentaires

Enregistrer un commentaire

Posts les plus consultés de ce blog

Standard error from Hessian Matrix... what can be done when problem occurs

An Hessian matrix is a square matrix of partial second order derivative. From the square-root of the inverse of the diagonal, it is possible to estimate the standard error of parameters. Let take an example: val=rnorm(100, mean=20, sd=5) # Return -ln L of values in val in Gaussian distribution with mean and sd in par fitnorm<-function(par, val) {   -sum(dnorm(val, par["mean"], par["sd"], log = TRUE)) } # Initial values for search p<-c(mean=20, sd=5) # fit the model result <- optim(par=p, fn=fitnorm, val=val, method="BFGS", hessian=TRUE) # Inverse the hessian matrix to get SE for each parameters mathessian <- result$hessian inversemathessian <- solve(mathessian) res <- sqrt(diag(inversemathessian)) # results data.frame(Mean=result$par, SE=res)   library(MASS) (r<-fitdistr(val, "normal")) It works well. However, the inverse of the matrix cannot be calculated if the second order derivative for some paramete
Lire la suite

stepAIC from package MASS with AICc

Due to a bug in dropterm() and addterm() in package MASS, it is impossible to use AICc with stepAIC() . Here is a solution. Enter these commands and now you can use: stepAIC(g0, criteria="AICc") or stepAIC(g0, criteria="BIC") For example: datax <- data.frame(y=rnorm(100), x1=rnorm(100),                     x2=rnorm(100), x3=rnorm(100), x4=rnorm(100)) g0 <- glm(y ~ x1+x2+x3+x4, data=datax) gx_AICc <- stepAIC(g0, criteria="AICc") gx_AIC <- stepAIC(g0, criteria="AIC") gx_AIC <- stepAIC(g0, criteria="BIC")
Lire la suite

Install treemix in ubuntu 20.04

Informations are available here  https://bitbucket.org/nygcresearch/treemix/downloads/ In terminal, enter for  Ubuntu 20.04: wget https://bitbucket.org/nygcresearch/treemix/downloads/treemix-1.13.tar.gz tar -xvf treemix-1.13.tar.gz sudo apt-get install libgsl-dev sudo apt install libboost-dev sudo apt install libboost-all-dev cd treemix-1.13/ ./configure make sudo make install And it works: $ treemix TreeMix v. 1.13 $Revision: 231 $ Options: -h display this help -i [file name] input file -o [stem] output stem (will be [stem].treeout.gz, [stem].cov.gz, [stem].modelcov.gz) -k [int] number of SNPs per block for estimation of covariance matrix (1) -global Do a round of global rearrangements after adding all populations -tf [file name] Read the tree topology from a file, rather than estimating it -m [int] number of migration edges to add (0) -root [string] comma-delimited list of populations to set on one side of the root (for migration) -g [vertices file name] [edges file name] read the g
Lire la suite