R/plot.airshedM.R

In tianjialiu/HyAirshed-RPackage: FRP in Daily and Seasonal HYSPLIT-derived Airsheds

#' Create Plots for HYSPLIT-derived Monthly Airshed
#' 
#' Create plots for monthly airsheds with daily HYSPLIT back trajectories for a given point location
#' @usage plot.airshedM(loc, lon, lat, month_select, year_select,
#'    tYears=year_select[1]:year_select[length(year_select)],
#'    hy_alt, alt.adj=0, sPix=0.25, nPix=100, borderName="ne_10m_admin_0_countries",
#'    x.adj=FALSE, zoom=4.5, x.just=0, y.just=0, nKM=150,
#'    traj.col=c("#a6cee3","#fb9a99","#cab2d6"), traj.lwd=0.25,
#'    col.adj=c(0,0,0,0), scale.adj=c(0,0), bar.adj=c(0,0,0,0),
#'    bar2.adj=c(0,0,0,0), caption=NA, caption.adj=c(0,0),
#'    nAlt.adj=c(0,0), legend.lab=hy_alt, legend.adj=c(0,0),
#'    txt_size=0.75, ...)
#' @param loc character. Name of location, starting point, used in naming subdirectories and files
#' @param lon numeric. Longitude of starting point
#' @param lat numeric. Latitude of starting point
#' @param month_select numeric vector. Months to be run
#' @param year_select numeric vector. Years to be run
#' @param tYears numeric vector. Total range of years in study, ignoring missing years [default = year_select[1]:year_select[length(year_select)], which is the years from the first year to last year in year_select
#' @param hy_alt character string or vector. Name(s) of hysplit heights to be run
#' @param alt.adj numeric. Adjustment if more than one location [default is 0]
#' @param sPix numeric. Size of pixel, in degrees [default is 0.3]
#' @param nPix numeric. Number of pixels in horizontal and vertical direction [default is 100]
#' @param borderName character. Name of shapefile located in the "shapefiles" folder of project home directory to plot country/sub-country borders [default is "ne_10m_admin_0_countries"]
#' @param x.adj logical. Adjust x axis latitude to 0-360? default is FALSE (-180 to 180)
#' @param zoom numeric. Window of the plot, in degrees, from the lon, lat of the starting location in all 4 directions [default is 4.5]
#' @param x.just numeric. Adjust window laterally, positive is right and negative is left [default is 0]
#' @param y.just numeric. Adjust window vertically, positive is right and negative is left [default is 0]
#' @param traj.col character string or vector in html or R color name. Colors of trajectories [default is c("#a6cee3","#fb9a99","#cab2d6")]
#' @param traj.lwd numeric. Width of trajectory lines [default is 0.3]
#' @param nKM numeric. Length of scale bar, in kilometers [default is 150000]
#' @param col.adj numeric vector, length=4. Adjust color bar position/size (format is c(x.min, x.max, y.min, y.max) adjustment) [default is c(0,0,0,0)]
#' @param scale.adj numeric vector, length=2. Adjust scale bar position (format is c(x, y) adjustment), [default is c(0,0)]
#' @param bar.adj numeric vector, length=4. Adjust timestamp bar position/size (format is c(x.min, x.max, y.min, y.max) adjustment), [default is c(0,0,0,0)]
#' @param caption character or expression. Name of caption, [default is NA], which is no caption
#' @param caption.adj numeric vector, length=2. Adjust caption position (format is c(x, y) adjustment) [default is c(0,0)]
#' @param nAlt.adj numeric vector, length=2. Adjust number of trajectories label position (format is c(x, y) adjustment) [default is c(0,0)]
#' @param legend.lab character vector. Labels for back trajectories [default is hy_alt]
#' @param legend.adj numeric vector, length=2. Adjust legend position (format is c(x, y) adjustment), [default is c(0,0)]
#' @param txt_size numeric. Adjust text size of plot labels, based on size of axis tick labels [default is 0.75]
#' @param ... see global arguments: pointDir, ask_home
#' @keywords hysplit airshed monthly plot
#' @export
#' @examples 
#' #creates plots for monthly delhi airshed for Oct-Nov, 2007-2013
#' plot.airshedM(loc="delhi",lon=77.2090,lat=28.6139,month_select=10:11,year_select=2007:2013,hy_alt="0.5 km",x.just=1.5,y.just=-1.5,caption="Delhi Monthly Airshed")

plot.airshedM <- function(loc, lon, lat, month_select, year_select,
                          tYears=year_select[1]:year_select[length(year_select)],
                          hy_alt, alt.adj=0, sPix=0.25, nPix=100, borderName="ne_10m_admin_0_countries",
                          x.adj=FALSE, zoom=4.5, x.just=0, y.just=0, nKM=150,
                          traj.col=c("#a6cee3","#fb9a99","#cab2d6"),
                          traj.lwd=0.25, col.adj=c(0,0,0,0), scale.adj=c(0,0),
                          bar.adj=c(0,0,0,0), bar2.adj=c(0,0,0,0), caption=NA,
                          caption.adj=c(0,0), nAlt.adj=c(0,0), legend.lab=hy_alt,
                          legend.adj=c(0,0), txt_size=0.75,
                          pointDir=paste0("./",loc,"/"), ask_home=TRUE) {
  
  options(error=NULL)
  
  #folder locations
  Home <- getwd()
  
  if (ask_home=="TRUE") {
    DirAns <- readline(paste("Is the project home directory:",Home,"? (enter y/n) "))
    if (DirAns!="y") {stop("Set project home directory with setwd()")}
  }
  
  #INPUT: hysplit gis & txt files
  #subfolders = iYear (e.g. 2007) -> month (e.g. oct)
  trajHome <- paste0(pointDir,"hysplit/")
  ifelse(!file.exists(trajHome),stop("missing hysplit coordinates folder"),FALSE)
  
  #INPUT: airsheds
  airDir <- paste0("./",loc,"/airshed_monthly/")
  airHome <- paste0(airDir,"shapefiles/")
  
  #OUTPUT: airshed plots
  outputHome <- paste0(airDir,"plots/")
  ifelse(!file.exists(outputHome),dir.create(outputHome,recursive=T),FALSE)
  
  border <- readShapePoly(paste0("./shapefiles/",borderName),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
  print("Okay, wait for function to finish.")
  
  mdeg <- function(reflat) {
    rlat = reflat*pi/180;
    mlat = 111412.84*cos(rlat)-93.5*cos(3*rlat)
    return(mlat)
  }
  
  gt <- GridTopology(c(lon-(sPix*nPix/2),lat-(sPix*nPix/2)),c(sPix,sPix),c(nPix,nPix))
  grd <- SpatialGrid(gt,proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
  spix <- as(grd,"SpatialPixels")
  spol <- as(spix,"SpatialPolygons")
  spol$coor <- coordinates(spol)
  
  for (iYear in year_select) {
    for (iMonth in month_select) {
      x.labs <- seq(-180,180,2); x.labs_s = seq(-180,180,0.5)
      y.labs <- seq(-180,180,2); y.labs_s = seq(-180,180,0.5)
      
      x.labsN <- c()
      for(i in 1:length(x.labs)) {
        if (x.adj==FALSE) {
          if (x.labs[i]>0) {x.labsN[i] <- paste0(x.labs[i],"°E")}
          if (x.labs[i]<0) {x.labsN[i] <- paste0(abs(x.labs[i]),"°W")}
          if (x.labs[i]==0) {x.labsN[i] <- 0}
        }
        
        if (x.adj==TRUE) {
          if (x.labs[i]>180) {x.labsN[i] <- paste0(x.labs[i]-360,"°W")}
          if (x.labs[i]<180) {x.labsN[i] <- paste0(x.labs[i],"°E")}
          if (x.labs[i]==180) {x.labsN[i] <- 0}
        }
      }
      
      y.labsN <- c()
      for(i in 1:length(y.labs)) {
        if (y.labs[i]>0) {y.labsN[i] <- paste0(y.labs[i],"°N")}
        if (y.labs[i]<0) {y.labsN[i] <- paste0(abs(y.labs[i]),"°S")}
        if (y.labs[i]==0) {y.labsN[i] <- 0}
      }
      
      box <- c(lon-zoom-x.just,lon+zoom-x.just,lat-zoom-y.just,lat+zoom-y.just)
      city_coor <- SpatialPoints(as.matrix(cbind(lon,lat)),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
      
      for (iMonth in 1:length(month_select)) {
        mon.actual <- month_select[iMonth] #month number
        mon.name <- tolower(month.abb)[mon.actual] #month number, starting from 1
        
        quartz(w=5.8,h=5.4);par(mar=c(2.2,2,3,1.2))
        plot(border,xlim=c(box[1],box[2]),ylim=c(box[3],box[4]),lwd=0.7,border="gray50");box()
        axis(s=1,tck=0.02,labels=NA,at=x.labs);axis(s=1,tck=0.01,labels=NA,at=x.labs_s);
        axis(s=1,lwd=0,line=-0.9,cex.axis=txt_size,las=1,x.labsN,at=x.labs)
        axis(s=2,tck=0.02,labels=NA,at=y.labs);axis(s=2,tck=0.01,labels=NA,at=y.labs_s)
        axis(s=2,lwd=0,line=-0.7,cex.axis=txt_size,las=3,labels=y.labsN,at=y.labs)
        axis(s=3,tck=0.02,labels=NA,at=x.labs);axis(s=3,tck=0.01,labels=NA,at=x.labs_s);
        axis(s=4,tck=0.02,labels=NA,at=y.labs);axis(s=4,tck=0.01,labels=NA,at=y.labs_s)
        
        airshedM <- readShapePoly(paste0(airHome,iYear,"/airshed",iYear*100+mon.actual),proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
        polyM <- gUnaryUnion(airshedM)
        
        #colormap for pixels
        color <- c("transparent",rev(paste0("gray",seq(1,99,1))),1)
        pixcol <- c()
        for (iPix in 1:length(airshedM)) {
          pixcol[iPix] = color[round(airshedM$den[iPix]*100/max(airshedM$den,na.rm=T))]
        }
        
        plot(polyM,add=T,border=1,lwd=2)
        plot(airshedM,add=T,col=pixcol,lwd=0.05)
        
        if ((iYear %% 4)!=0) { #non-leap year
          st.jul = c(1,32,60,91,121,152,182,213,244,274,305,335)
          end.jul = c(st.jul[-1]-1,365)
        }
        
        if ((iYear %% 4)==0) { #leap year
          st.jul = c(1,32,61,92,122,153,183,214,245,275,306,336)
          end.jul = c(st.jul[-1]-1,366)
        }
        
        for (iMonth in month_select) {
          
          TrajName2 <- list()
          for (iAltitude in 1:length(hy_alt)) {
            iAltitudeName <- hy_alt[iAltitude]
            shape.lat <- read.table(paste0(trajHome,loc,"_HysplitLat_",iYear,"_",iAltitudeName,".csv"),sep=",")
            shape.lon <- read.table(paste0(trajHome,loc,"_HysplitLon_",iYear,"_",iAltitudeName,".csv"),sep=",")
            
            nTraj <- shape.lat[,1][st.jul[iMonth]:end.jul[iMonth]]
            nTraj <- length(nTraj[!is.na(nTraj)])
            
            for(iDay in st.jul[iMonth]:end.jul[iMonth]) {
              if (length(shape.lat[iDay,][!is.na(shape.lat[iDay,])])>0) {
                shape.coor <- cbind(as.numeric(shape.lat[iDay,]),as.numeric(shape.lon[iDay,]))
                shape.coor <- apply(shape.coor,2,na.omit)
                shape <- SpatialPoints(shape.coor,proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
                traj <- as(shape,"SpatialLines")
                
                plot(shape,cex=0.1,lwd=traj.lwd,add=T,col=traj.col[iAltitude])
                plot(traj,lwd=0.3,add=T,col=traj.col[iAltitude])
              }
            }
          }
        }
        
        plot(city_coor,col=2,add=T,cex=1,pch=4,lwd=2)
        x <- 1:10; y <- 1:10; z <- outer(x,y,"+")
        obj <- list(x=x,y=y,z=z)
        ipos <- c(0.134,0.164,0.12,0.42)
        
        image.plot(obj,add=T,legend.only=T,zlim=c(0,100),col=color[-1],axes=F,
                   smallplot=ipos+col.adj,axis.args=list(at=seq(0,100,20),labels=NA,tck=-0.5))
        image.plot(obj,add=T,legend.only=T,zlim=c(0,100),col=color[-1],axes=F,
                   smallplot=ipos+col.adj,axis.args=list(at=seq(0,100,20),line=-0.4,labels=seq(0,100,20),cex.axis=0.7,lwd=0),
                   legend.args=list(text="Trajectories (% of Total)",side=2,font=2,line=0.3,cex=0.75,adj=0.5))
        
        par(new=T);plot(NA,NA,xlim=c(0,10),ylim=c(0,10),axes=F,xlab="",ylab="")
        
        spos <- c(1.5,0.1)+scale.adj
        
        segments(spos[1],spos[2],spos[1]+nKM*1000/mdeg(spos[2]),spos[2],lwd=1.2)
        segments(spos[1],spos[2]-0.12,spos[1],spos[2]+0.12,lwd=1.2)
        segments(spos[1]+nKM*1000/mdeg(spos[2]),spos[2]-0.12,spos[1]+nKM*1000/mdeg(spos[2]),spos[2]+0.12,lwd=1.2)
        text(mean(c(spos[1],spos[1]+nKM*1000/mdeg(spos[2]))),spos[2]+0.35,paste(nKM,"km"),cex=txt_size)
        
        bar <- c(-0.4,10.4,11.15,11.7)+bar.adj
        rect(bar[1],bar[3],bar[2],bar[4],xpd=T)
        rect(bar[1]+(bar[2]-bar[1])/length(tYears)*(iYear-min(tYears)),bar[3],bar[1]+(bar[2]-bar[1])/length(tYears)*(iYear-(min(tYears)-1)),bar[4],xpd=NA,col=1)
        text(bar[1]+(bar[2]-bar[1])/length(tYears)*(iYear-(min(tYears)-0.5)),mean(bar[3],bar[4])+0.5*(bar[4]-bar[3]),iYear,col="white",xpd=NA,cex=txt_size+0.05,adj=0.5)
        
        bar2 <- c(-0.4,10.4,10.525,11.075)+bar2.adj
        rect(bar2[1],bar2[3],bar2[2],bar2[4],xpd=T)
        rect(bar2[1]+(bar2[2]-bar2[1])/12*(mon.actual-1),bar2[3],bar2[1]+(bar2[2]-bar2[1])/12*(mon.actual),bar2[4],xpd=NA,col=1)
        text(bar2[1]+(bar2[2]-bar2[1])/12*(mon.actual-0.5),mean(bar2[3],bar2[4])+0.5*(bar2[4]-bar2[3]),month.abb[mon.actual],col="white",xpd=NA,cex=txt_size+0.05,adj=0.5)
        
        apos <- c(10.45,-1.2)+nAlt.adj
        text(apos[1],apos[2],paste0("Number of Trajectories Per Altitude = ",nTraj),xpd=NA,cex=0.6,font=3,adj=1)
        
        cpos <- c(0.05,9.9)+caption.adj
        if (is.na(caption)==FALSE) {
          text(cpos[1],cpos[2],caption,xpd=NA,cex=txt_size,font=2,adj=c(0,1))
        }
        
        lpos <- c(8.15,9.3)+legend.adj
        legend(lpos[1],lpos[2]-(length(hy_alt)-1)*0.2,legend=paste(legend.lab,""),lty=1,col=traj.col,cex=txt_size-0.05,title=expression(paste(bold("Altitudes"))),title.adj=0.25,yjust=0.5,bg="white")
        
        quartz.save(paste0(outputHome,"airshed",iYear*100+mon.actual,".pdf"),type="pdf")
        graphics.off()
      }
    }
    cat(iYear,": year",which(year_select==iYear),"of",length(year_select),as.character(Sys.time()),'\n')
  }
  print("Complete!")
}