PIHM.AnalysisR:

#' PIHM Analysis project.
#' Developed by Lele Shu( lele.shu at gmail.com  lzs157 at psu.edu ) 
#' Created by Wed Apr 15 20:25:45 EDT 2015
#'  Current version is for PIHM-MF or PIHM v2.4;
#' 

minfo <- function(bak,mapplot=PIHM.triplot,if.plot=TRUE){
    path=ModelInfopath;
    theFile <- file.path(path, paste(projectname,".",'minfo.txt',sep=''));
    mesh=readmesh(bak=bak);
    riv=readriv(bak=bak);
    K=suppressWarnings(Kcalib(calib.bak=bak,path=ModelInfopath,quiet=TRUE) );
    soil=readsoil(bak=FALSE)$soil;
    geol=readgeol(bak=FALSE)$geol;
    att=readatt(bak=bak);
#    lc=readlc(bak=bak)
    lc=att[,4]
    nm=mesh$size[1]
    np=mesh$size[2]
    nr=riv$River$size;
    no=length(riv$River$outlets);

    nsoil=nrow(soil);
    usoil=sort(unique(att[,2]));
    unsoil=length(usoil);

    ngeol=nrow(geol);
    ugeol=sort(unique(att[,3]));
    ungeol=length(ugeol)

    nlc=nrow(lc);
    ulc=sort(unique(att[,4]));
    unlc=length(ulc);

        
    para<-readpara(bak=bak);    
    start<-para[[which(grepl('^start',tolower(names(para))))]];
    end<-para[[which(grepl('^end',tolower(names(para))))]];
    dt <- para[[which(grepl('^rivflx1',tolower(names(para))))]];


    calib<-readcalib(bak=bak);
    
    write(paste('Name of Project = \t',projectname),file=theFile,append=FALSE);
    write(paste('From\t',start,'To\t',end ),file=theFile,append=TRUE);
    write('',file=theFile,append=TRUE);


    write(paste('Size of Mesh = \t',nm),file=theFile,append=TRUE);
    write(paste('Size of points = \t',np),file=theFile,append=TRUE);
    write(paste('Size of river = \t',nr),file=theFile,append=TRUE);
    write(paste('Size of outlets = \t',no),file=theFile,append=TRUE);
    write('',file=theFile,append=TRUE);

    write(paste('Size of soil = \t',unsoil,"/",nsoil),file=theFile,append=TRUE);
    write(paste('Size of geol = \t',ungeol,"/",ngeol),file=theFile,append=TRUE); 
    write(paste('Size of LC = \t',unlc,"/",nlc),file=theFile,append=TRUE);
    
        write('',file=theFile,append=TRUE);
    write(paste('Unique of soil = \t',paste(usoil,collapse=' ')),file=theFile,append=TRUE);
    write(paste('Unique of geol = \t',paste(ugeol,collapse=' ')),file=theFile,append=TRUE);
    write(paste('Unique of LC = \t',paste(ulc,collapse=' ')),file=theFile,append=TRUE);
    
        write('\nm/day',file=theFile,append=TRUE);
    write(paste(colnames(K$cKsoil),collapse=' '),file=theFile,append=TRUE);
    suppressWarnings(write(K$cKsoil[nrow(K$cKsoil),],file=theFile,append=TRUE))

        write('\nm/day',file=theFile,append=TRUE);
    write(paste(colnames(K$cKgeol),collapse=' '),file=theFile,append=TRUE);
    suppressWarnings(write(K$cKgeol[nrow(K$cKgeol),],file=theFile,append=TRUE))

        write('\nm/day',file=theFile,append=TRUE);
    suppressWarnings(write.table(K$Krivm,file=theFile,row.names=TRUE,col.names=TRUE,,quote=FALSE,append=TRUE))


    
 if(if.plot){
    PIHM.triplot(data=att[,2],fn=paste(projectname,'_Soil_map.png',sep=''),name='Soil Type',title=paste('Soil Map of', projectname), path=path);
    PIHM.triplot(data=att[,3],fn=paste(projectname,'_Geol_map.png',sep=''),name='Geol Type',title=paste('Geology Map of', projectname), path=path);
    PIHM.triplot(data=att[,4],fn=paste(projectname,'_LC_map.png',sep=''),name='Land Cover Type',title=paste('Land Cover Map of', projectname), path=path);

    mapplot(data=soildepth(),fn=paste(projectname,'_Soil_Depth.png',sep=''),name='Soil depth(m)',title=paste('Soil depth of', projectname), path=path);
    
    mapplot(data=soil$soil[att[,2],'INFK']*86400 * calib['KINF'] ,fn=paste(projectname,'_Kinf.png',sep=''),name='K infilration(m/day)',title=paste('Conductivity of Infiltration', projectname), path=path);
    mapplot(data=geol$geol[att[,3],'SATHK']*86400* calib['KSATH'] ,fn=paste(projectname,'_KH.png',sep=''),name='K Horizontal(m/day)',title=paste('Horizontal Conductivity of', projectname), path=path);
    mapplot(data=geol$geol[att[,3],'SATVK']*86400* calib['KSATV'] ,fn=paste(projectname,'_KV.png',sep=''),name='K Vertical(m/day)',title=paste('Vertical Conductivity of', projectname), path=path);
    
    #VegFrac
    mapplot(data=lc[att[,4],2]* calib['VEGFRAC'] * 100 ,fn=paste(projectname,'_VegFrac.png',sep=''),name='VEGFRAC(%)',title=paste('Vegetation fraction of', projectname), path=path);
    #LAI max
    mapplot(data=lc[att[,4],9] ,fn=paste(projectname,'_LAI_max.png',sep=''),name='LAI max(m2/m2)',title=paste('Max LAI of', projectname), path=path);
    #Roughness
    mapplot(data=lc[att[,4],16] ,fn=paste(projectname,'_Roughness.png',sep=''),name='Roughness(m)',title=paste('Roughness of', projectname), path=path);

 }



    

    
}
#============ #============
#============ #============


landcovertable <-function(bak=TRUE,class=FALSE, output=ModelInfopath, att)
{
    #lc=readlc(bak=bak);
    if (missing(att)){
        att=readatt(bak=bak);
    }
    lc=att[,4];
   
    nlc=nrow(lc);
    ulc=sort(unique(att[,4]));
    unlc=length(ulc);

    iarea=readarea(bak=bak);
    area=sum(iarea)
    lca=matrix(0,ncol=2,nrow=unlc)
    for (i in 1:unlc){
        lca[i,1] = sum(iarea[which(att[,4]==ulc[i])]);
        lca[i,2] = 100 * lca[i,1]/area
    }
    lca=rbind(lca,colSums(lca))
    cname=c('Area(m2)', 'Percentage')
    colnames(lca)=cname
    rname=c(as.character(ulc), 'sum')
    rownames(lca)=rname
    
    theFile <- file.path(path=output, paste(projectname,".",'lctable.txt',sep=''));
    write.table(lca,file=theFile,append=FALSE);

    if (class){
        
        cl=matrix(0,ncol=2,nrow=5); 
        type=list(c(23,24), 
                  c(20,21,22),
                  c(40,41,42),
                  c(80,81,82)
                  )
        # 1-Urban 23 24;
        # 2-SubUrban 20 21 22
        # 3-Forest 40 41 42
        # 4-Agriculture 80 81 82
        # 5-Other  11 30x 90x
        i=1;
        id=which(ulc>=23 & ulc<30)
        cl[i,]=colSums(lca[id,]);
        ids=id;
        
        i=2;
        id=which(ulc>=20 & ulc<23)
        cl[i,]=colSums(lca[id,]);
        ids=c(ids,id)

        i=3;
        id=which(ulc>=40 & ulc<50)
        cl[i,]=colSums(lca[id,]);
        ids=c(ids,id)

        i=4;
        id=which(ulc>=80 & ulc<90)
        cl[i,]=colSums(lca[id,]);
        ids=c(ids,id)

        i=5;
        ss=1:unlc
        id=ss[!(ss %in% ids)]
        cl[i,]=colSums(lca[id,]);
        
        cl=rbind(cl,colSums(cl));

        rname=c('Urban','SubUrban', 'Forest', 'Agricultrue', 'Other','SUM')
        colnames(cl)=cname;
        rownames(cl)=rname;
        write.table(cl,file=theFile,append=TRUE);
        ret=cl;
    }else{
        ret = lca
    }

    return(ret);
    
}



getelev <- function(bak=TRUE, bed =FALSE){
    mesh <- readmesh(bak=bak);
    msh <- mesh$mesh;
    pts <- mesh$points;
    nabr<-mesh$mesh[,5:7]   #nabor or each cell.
    node<-mesh$mesh[,2:4]    #vetex of each cell.
    
    xv       =cbind(pts[node[,1],2],pts[node[,2],2],pts[node[,3],2]); #end with v, means vertex of triangles.
    yv       =cbind(pts[node[,1],3],pts[node[,2],3],pts[node[,3],3]);
    zbedv    =cbind(pts[node[,1],4],pts[node[,2],4],pts[node[,3],4]);
    zsurfv   =cbind(pts[node[,1],5],pts[node[,2],5],pts[node[,3],5]);
    
    xc      =rowMeans(xv);      #end of c, means centronid of the triangles.
    yc      =rowMeans(yv);
    zbedc   =rowMeans(zbedv);
    zsurfc  =rowMeans(zsurfv);
    if ( !bed){
        return(zsurfc);
    }else{
        return(zbedc); 
    }
    
}

soildepth <- function(bak=TRUE){ 
    mesh <- readmesh(bak=bak);
    msh <- mesh$mesh;
    pts <- mesh$points;

    nabr<-mesh$mesh[,5:7]   #nabor or each cell.
    node<-mesh$mesh[,2:4]    #vetex of each cell.

    
    xv       =cbind(pts[node[,1],2],pts[node[,2],2],pts[node[,3],2]); #end with v, means vertex of triangles.
    yv       =cbind(pts[node[,1],3],pts[node[,2],3],pts[node[,3],3]);
    zbedv    =cbind(pts[node[,1],4],pts[node[,2],4],pts[node[,3],4]);
    zsurfv   =cbind(pts[node[,1],5],pts[node[,2],5],pts[node[,3],5]);

    
    sd= rowMeans(zsurfv)-rowMeans(zbedv);
    return(sd);
}