Manipulate ncdf files




library("HelpersMG")


# Guatemala map
# get the ncdf file here for 2010, 2011, 2012
# <a href="http://www.esrl.noaa.gov/psd/cgi-bin/db_search/DBListFiles.pl?did=132&tid=33865&vid=2423">http://www.esrl.noaa.gov/psd/cgi-bin/db_search/DBListFiles.pl?did=132&tid=33865&vid=2423</a>
# ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.2012.v2.nc
# ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.2011.v2.nc
# ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.2010.v2.nc

url <- "ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.2012.v2.nc"
dest <- paste(Sys.getenv("HOME"), "/sst.day.mean.2012.v2.nc", sep="")
download.file(url, dest)

url <- "ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.2011.v2.nc"
dest <- paste(Sys.getenv("HOME"), "/sst.day.mean.2011.v2.nc", sep="")
download.file(url, dest)

url <- "ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.2010.v2.nc"
dest <- paste(Sys.getenv("HOME"), "/sst.day.mean.2010.v2.nc", sep="")
download.file(url, dest)

library("ncdf4")
dta2010<-nc_open(paste(Sys.getenv("HOME"), "/sst.day.mean.2010.v2.nc", sep=""))
dta2011<-nc_open(paste(Sys.getenv("HOME"), "/sst.day.mean.2011.v2.nc", sep=""))
dta2012<-nc_open(paste(Sys.getenv("HOME"), "/sst.day.mean.2012.v2.nc", sep=""))

# get the daily sst temperature for Pacific Guatemala
indice <- ind_long_lat(ncdf=dta2010, long=-90.25, lat=13.5)
Monterrico2010<-ncvar_get(dta2010, varid="sst", start=c(indice,1), count=c(1, 1, dta2010$dim$time$len))
Monterrico2011<-ncvar_get(dta2011, varid="sst", start=c(indice,1), count=c(1, 1, dta2011$dim$time$len))
Monterrico2012<-ncvar_get(dta2012, varid="sst", start=c(indice,1), count=c(1, 1, dta2012$dim$time$len))

# 366 because 2012 is bisextille year
Monterrico2012<-c(Monterrico2012, rep(NA, 366-length(Monterrico2012)))
Monterrico<-c(Monterrico2010, Monterrico2011, Monterrico2012)
names(Monterrico)<-seq(from=as.Date("2010-01-01"), by="1 days", length.out=length(Monterrico))

# plot the SST temperatures
plot(as.Date(names(Monterrico)), Monterrico, xlab="Days", ylab="SST C", type="l", bty="n", axes=FALSE)

axis(side = 1, at = seq(from=as.Date("2010-01-01"), by="1 month",
                        length.out=length(Monterrico)),
     labels=c(labels=rep(c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"),
                         round(length(Monterrico)/365.25)+1), ((round(length(Monterrico)/365.25)+1)*12+1):(length(Monterrico))))
axis(side=2)

library("fields")

# Map centered to Pacific
day <- as.Date("2010-02-01")
carte<-ncvar_get(dta2010, varid="sst", start=c(1,1,as.numeric(day-as.Date("2010-01-01")+1)), count=c(dta2010$dim$lon$len, dta2010$dim$lat$len, 1))
image(dta2010$dim$lon$vals, dta2010$dim$lat$vals, carte, col=rainbow(128), xlab="Longitude", ylab="Latitude");segments(0,0,360,0, lwd=1, lty=2)

# Arrow to Guatemala
HelpersMG:::.Arrows(x1=-110%%360, 5, x0=-90.25%%360, y0=13.5, y1=13.5)


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