R scripts

R Functions for common operations

#=================================================================================
# Compute dominant heights
#=================================================================================
calc_hd <- function(anagUT,cav,ipso,id_ut="id_ut",sup_ha="sup_ha",d130="d130",htot="htot"){
# Computes dominant height (h_dom) for each area unit (id_ut = Unità Territoriale)
#input: 3 dataframe, registry (anagUT), dbh tally (cav) and height tally (ipso),

#       where area units are identified always with the same name (the value of <id_ut>)
    if(length(unique(anagUT[,id_ut]))!=nrow(anagUT)) print("ERROR id_ut is not key in <anagUT>")
    anagUT$n_dom <- anagUT[,sup_ha]*100
    hdom_df<-NULL
    for(ut in anagUT[,id_ut]){
        n_dom<-anagUT$n_dom[anagUT[,id_ut]==ut]
        cav0 <- sort(cav[cav[,id_ut]==ut,d130])
        if (lenght(cav0) < n_dom) d_dom <- NA else {
            s <- sum(cav0[1:n_dom]^2)+(n_dom-floor(n_dom))*cav0[n_dom]^2
            d_dom<-sqrt(s/n_dom)   #diam di area basimetrica (dominante) media
            }
        ipso0<-ipso[ipso[,id_ut]==ut,c(id_ut,d130,htot)]
        if (nrow(ipso0) < 3) h_dom <- NA else {
            names(ipso0) <- c("id_ut","d130","htot")
            fit<-lm(htot ~ I(sqrt(d130)), data=ipso0)
            h_dom<-predict(fit,data.frame(d130=d_dom))
            }
        hdom_df<-rbind(hdom_df,data.frame(id_ut=ut,b_0=fit$coefficients[1],b_1=fit$coefficients[2],d_dom=d_dom,h_dom=h_dom))
        }
    hdom_df1 <- merge(anagUT[,c(id_ut,"n_dom")],hdom_df,by.x=id_ut,by.y="id_ut")
    # Tavola delle altezze dominanti per ogni <UT>
    return(hdom_df1)
    }
#=================================================================================