R/read.R

In jobonaf/pesco: Data fusion for air quality data

## functions to read:
##  - observed data
##  - CTM output
##  - stations' metadata
##  - proxies

## read NetCDF as provided by Arpa Emilia-Romagna
read.ncdf.arpaer <- function(con=NULL, pollutant="pm10", lev=1,
                             tz.in="UTC", tz.out="Africa/Algiers") {
  if(class(con)!="ncdf") {   ## existing file or connection
    nc   <- open.ncdf(con)
  } else {                    ## ncdf object
    nc <- con
  }
  lon  <- get.var.ncdf(nc,varid="lon")
  lat  <- get.var.ncdf(nc,varid="lat")
  Time <- as.POSIXct(get.var.ncdf(nc,varid="Times"), 
                     format="%Y-%m-%d_%H:%M:%S", tz=tz.in)
  Time <- tz.change(x=Time,tz.in=tz.in,tz.out=tz.out)
  opts <- c(pollutant, toupper(pollutant), tolower(pollutant))
  vars <- names(nc$var)
  pp <- intersect(opts, vars)
  if(length(pp)==0) stop(paste(pollutant,"not found in",con))
  value  <- get.var.ncdf(nc,varid=pp)
  if(length(dim(value))==4) {
    value <- value[,,lev,]
  } else if (length(dim(value))!=3) {
    stop(paste(pollutant,"must have dimensions X-Y-Z-T or X-Y-T"))
  }
  if (class(con)!="ncdf") {
    close.ncdf(nc)
  }
  coords <- ll2utm.grid(lat,lon)
  Coords <- list(x=matrix(coords$x,nrow=nrow(lon),ncol=ncol(lon)),
                 y=matrix(coords$y,nrow=nrow(lon),ncol=ncol(lon)))
  out <- list(coords=Coords, time=Time, data=value)
}


read.grads <- function(filectl) {
  
  # manages file name
  dum  <- strsplit(filectl,"/",fixed=TRUE)[[1]]
  file <- dum[length(dum)]
  path <- NULL
  if (length(dum)> 1) {
    for (i in 1:(length(dum)-1)) {
      path <- paste(path,dum[i],sep="/")
    } 
  }       
  nlines <-  length(scan(filectl,what="character", flush = TRUE))
  ctl <- scan(filectl,what="character",sep="\n")
  
  # reads ctl file
  ok <- 0
  for (iline in 1:length(ctl)) {
    line <- ctl[iline]
    if (xgrep("dset",small(line))) {
      dum <- subwrd(line,2)
      if ( xgrep("^",dum) ) { 
        filebin <- strsplit(dum,"^",fixed=TRUE)[[1]][2]
        if(length(path)>0) filebin <- paste(path,filebin,sep="/")
      } else {
        filebin <- dum
      }
      ok <- ok + 1 
    }
    if (xgrep("undef",small(line))) {
      rmis <- as.numeric(subwrd(line,2))
      ok <- ok + 1 
    }
    if (xgrep("xdef",small(line))) {
      nx   <- as.integer(subwrd(line,2))
      ok   <- ok + as.integer(subwrd(small(line),3)=="linear")
      xmin <- as.numeric(subwrd(line,4))
      dx   <- as.numeric(subwrd(line,5))
    }
    if (xgrep("ydef",small(line))) {
      ny   <- as.integer(subwrd(line,2))
      ok   <- ok + as.integer(subwrd(small(line),3)=="linear")
      ymin <- as.numeric(subwrd(line,4))
      dy   <- as.numeric(subwrd(line,5))
    }
    if (xgrep("zdef",small(line))) {
      nz   <- as.integer(subwrd(line,2))
      ztyp <- subwrd(line,3)
      if(ztyp=="linear") {
        zmin <- as.numeric(subwrd(line,4))
        dz   <- as.numeric(subwrd(line,5))
        ok <- ok + 1 
      } 
      if(small(ztyp)=="levels") {
        zdim <- as.numeric(subwrd(line,4:(3+nz)))
        ok <- ok + 1 
      } 
    }
    if (xgrep("tdef",small(line))) {
      nt   <- as.integer(subwrd(line,2))
      ok   <- ok + as.integer(subwrd(small(line),3)=="linear")
      tmin <- small(subwrd(line,4))
      dt   <- small(subwrd(line,5))
    }
    if (xgrep("vars",small(line))) {
      ok <- ok + 1 
      if(ok==7) {
        nv  <- as.integer(subwrd(line,2))
        var <- vector(mode="character",length=nv)
        for (ivar in 1:nv) {
          dum <- ctl[(iline+ivar)]
          var[ivar] <- subwrd(dum,1)
        }
      } else if(ok>8) {
        stop(paste("Problems with",filectl))
      }   
    }
  }
  
  # manages dimensions names
  xdim <- as.character(seq(xmin,length=nx,by=dx))
  ydim <- as.character(seq(ymin,length=ny,by=dy))
  if(ztyp=="linear") zdim <- as.character(seq(zmin,length=nz,by=dz))
  vdim <- as.character(var)
  
  # manages dates and time 
  date.fmt=c(dates="ymd", times="hms")
  tl <- nchar(tmin)
  if (tl==4) tmin <- paste("00:00z01jan",tmin,sep="")
  if (tl==7) tmin <- paste("00:00z01",tmin,sep="")
  if (tl==9) tmin <- paste("00:00z",tmin,sep="")
  if (tl==12)tmin <- paste(substr(tmin,1,2),":00",substr(tmin,3,12),sep="")
  if (date.lang()=="it") tmin <- en2it.date(tmin)
  if (date.lang()=="en") tmin <- it2en.date(tmin)
  timemin <- paste(substr(tmin,1,5),":00",sep="")
  datemin <- substr(tmin,7,15)
  datemin <- as.character(as.Date(datemin,format="%d%b%Y"),format="%m/%d/%y")
  start <- chron(datemin, timemin, format = c(dates = "m/d/y", times = "h:m:s"))
  
  sl <- nchar(dt)
  steptyp <- substr(dt,(sl-1),(sl))
  nstep   <- as.integer(substr(dt,1,(sl-2)))
  if (steptyp=="yr") {
    dum <- seq.dates(start,by="years",length.=(nt*nstep))
    dum <- dum[seq(1,length=nstep,by=nt)]
    dum <- chron(dum,out.format=date.fmt)
    tdim <- as.character(dum)
  } 
  if (steptyp=="mo") {
    dum <- seq.dates(start,by="months",length.=(nt*nstep))
    dum <- dum[seq(1,length=nt,by=nstep)]
    dum <- chron(dum,out.format=date.fmt)
    tdim <- as.character(dum)
  } 
  if (steptyp=="dy") {tdim <- as.character(chron(start+seq(0,length=nt,by=nstep       ), out.format=date.fmt))}
  if (steptyp=="hr") {tdim <- as.character(chron(start+seq(0,length=nt,by=(nstep/24)  ), out.format=date.fmt))}
  if (steptyp=="mn") {tdim <- as.character(chron(start+seq(0,length=nt,by=(nstep/1440)), out.format=date.fmt))}
  
  # reads GrADS binary file
  rm(dum)
  bits=4
  dum <- readBin(filebin,"numeric",n=nv*nx*ny*nt*nz,endian="little",size=bits)
  is.na(dum) <- dum==rmis
  dat <- array (dum, dim=c(nx,ny,nz,nv,nt), 
                dimnames=list(x=xdim,y=ydim,z=zdim,var=vdim,t=tdim))
  return(dat)
}

# Funzione R per la lettura di file in formato estra_qaria.
read.qaria <- function(file) {
  hourly <- scan(file,skip=2,what="",nlines=1)[4]=="hh"
  dat <- read.fwf(file = file,
                  na.strings=c("    -9999.0","     -999.0"),
                  widths=c(4,rep(3,2+as.integer(hourly)),11),
                  header=FALSE,skip=3)
  header=scan(file, skip=2, what="", n=4++as.integer(hourly))
  names(dat) <- header
  strTime <- paste(dat$aaaa,dat$mm,dat$gg,sep="-")
  if(hourly) {
    strTime <- paste(strTime," ",dat$hh,":00",sep="")
  } else {
    strTime <- paste(strTime," 00:00",sep="")    
  }
  Time <- as.POSIXct(strTime,tz="Africa/Algiers")
  out <- data.frame(Time=Time, Value=dat[,ncol(dat)])
  names(out)[2] <- names(dat)[ncol(dat)]
  return(out)
}

## legge output di query SQL, tra cui le anagrafiche di stazione
read.sql <- function(file) {
  check <- substr(scan(file,what="",sep="\n",blank.lines.skip = F),1,3)
  nsk <- which(check=="---")[1]
  row <- scan(file,skip=nsk-1,nlines=1,what="")
  hea <- scan(file,skip=nsk-2,nlines=1,what="")
  data <- read.fwf(file,widths=nchar(row)+1,skip=nsk)
  colnames(data) <- hea
  last <- which(rowSums(!is.na(data))<1)[1]-1
  last <- min(last, nrow(data), na.rm=T)
  data <- data[1:last,]
  return(data)
}

## legge una lista di file qaria e un file di anagrafica
## e assembla tutto in una tabella "long"
qaria2long <- function(datafiles,
                       anafile,
                       codes=NULL) {
  nf <- length(datafiles)
  ana <- read.sql(anafile)
  for (i in 1:nf) {
    if(is.null(codes)) {
      code <- substr(datafiles[i],
                      rev(gregexpr("_",datafiles[i])[[1]])[1]+1,
                      rev(gregexpr("\\.",datafiles[i])[[1]])[1]-1)
    } else {
      code <- codes[i]
    }
    dat <- read.qaria(datafiles[i])
    idx <- which(as.numeric(code)==as.numeric(as.character(ana$CODE)))
    if(length(idx)==1) {
      Dat <- data.frame(dat, 
                        Name=rep(trim(as.character(ana[idx,"NAME"])),nrow(dat)),
                        Municipality=rep(trim(as.character(ana[idx,"MUNICIPALITY"])),nrow(dat)),
                        Code=rep(code,nrow(dat)),
                        Lat=rep(ana[idx,"LAT"],nrow(dat)),
                        Lon=rep(ana[idx,"LON"],nrow(dat)),
                        Elev=rep(ana[idx,"ELEV"],nrow(dat)),
                        Type=rep(ana[idx,"CLA"],nrow(dat)))
      if(!exists("Tab")) {
        Tab <- Dat
      } else {
        Tab <- rbind(Tab,Dat)
      }
    }
  }
  Tab$Elev[Tab$Elev < -100] <- NA
  Tab <- droplevels(Tab)
  return(Tab)
}


read.field <- function(file, coords.col=1:2, data.col=3, coords.fact=1, ...) {
  dat <- read.table(file,...)
  out <- list(coords=list(x=dat[,coords.col[1]]*coords.fact,
                          y=dat[,coords.col[2]]*coords.fact),
              data=dat[,data.col])
  return(out)
}