R/funkcje.r
In PawelKra/wskCI: WskaŸniki konkurencji dla drzew
#!/usr/bin/env Rscript
#autor: Pawel Krapiec
# biblioteka do obliczen odleglosci pkt od lini
require(spatstat)
# ustawienia srodowiska R
options(supen = 100)
# niezbedne zmienne do obliczen
typPow <- 'K' # K - kolowki; I - powierzchnie inne
kasNiep <- 'T' # czy kasowac pola nie bedace wykorzystywane do obliczen
# Domyslna wspolrzedna srodka powierchni kolowej
X <- 0 # wspolrzedna X powierzchni
Y <- 0 # wspolrzedna Y powierzchni
#--------------------------------
#DEFINICJE NAZW POL W PLIKU WEJSCIOWYM
nrPow <- 'nr.pow.' # nazwa kolumny z numerem powierzchni
d13 <- 'd' # nazwa kolumny z piersnica
h <- 'h' # nazwa kolumny z wysokoscia drzewa
hk <- 'hk' # nazwa kolumny z wysokoscia nasady korony
azymut <- 'azymut' # nazwa kolumny z kierunkiem na ktorym rosnie drzewo
odlDrzew <- 'L' # nazwa kolumny z odlegloscia drzewa od pkt
spadek <- 'nachylenie' # spadek powierzchni probnej
wystawa <- 'wystawa' # wystawa powierzchni
gat <- 'gat.' # gatunek drzewa (So , Św, Db, Jd, Bk, Db, Gb, itd.)
nrDrzewa <- 'l.p.' # numer porzadkowy drzewa na powierzchni
R <- 'R' # promien powierzchni kolowej
rKorony <- 'szer.kor.' # kolumna z rzeczywistym promieniem korony
rKoronyModel <- 'rKoronyModel' # kolumna z obliczonym promieniem korony
Xorg <- 'X' # wspolrzedna X drzewa
Yorg <- 'Y' # wspolrzedna Y drzewa
Zorg <- 'Z' # wspolrzedna Y drzewa
Xsr <- 'Xsr' # wspolrzedna X srodka kolowki
Ysr <- 'Ysr' # wspolrzedna Y srodka kolowki
# Zmienne programowe - nie zmieniac
Xo <- 'Xo' # obliczona wspolrzedna X
Yo <- 'Yo' # obliczona wspolrzedna Y
wspBL2XY <- function(x, y, kat, dist) {
# obliczenia dla powierzchni kolowych, przelicza L i azymut na X, Y
if (kat<=90) katN <- 90 - kat
else if (kat>90 & kat<=270) katN <- 360 - kat + 90
else if (kat>270) katN <- 180 - (kat - 270)
katrad <- katN * 0.017453292519
x_przelicz <- round(x + dist * cos(katrad), 2)
y_przelicz <- round( y + dist * sin(katrad), 2)
wyn <- c(x_przelicz , y_przelicz)
return(wyn)
}
zmienne <- function() {
# funkcja zwraca nazwy i wartosci zadeklarowanych zmiennych urzywanych do
# odczytu danych z plikow i obliczen
zm_gr <- c(typPow ,kasNiep ,X ,Y, nrPow ,d13 ,h ,hk ,azymut ,
odlDrzew ,spadek ,wystawa ,gat ,nrDrzewa ,R ,rKorony ,rKoronyModel,
Xorg ,Yorg ,Zorg, Xsr ,Ysr)
zmienne_ops <- c('typPow','kasNiep','X','Y','nrPow','d13','h',
'hk','azymut','odlDrzew','spadek','wystawa','gat',
'nrDrzewa','R','rKorony','rKoronyModel','Xorg','Yorg', 'Zorg',
'Xsr','Ysr')
print(paste(zmienne_ops, zm_gr, sep = " = "))
}
d2r <- function(wart) {
# przelicza podana wartosc na radiany
if (!is.numeric(wart)) stop('Podana wartosc nie jest liczba, nie mozna przeliczyc na radiany')
w <- wart * pi/180
return(w)
}
odlAB <- function(xa, ya, xb, yb) {
# obliczenie odleglosci miedzy 2 punktami
odl <- sqrt((xa - xb)**2 + (ya - yb)**2)
return(odl)
}
pktObwodu <- function( x1 = 0, y1 = 0, r=1, npoints = 100){
# oblicza punkty z okregu o zadanym promieniu
tt <- seq(0,2*pi,length.out = npoints)
xx <- x1 + r * cos(tt)
yy <- y1 + r * sin(tt)
return(data.frame(x = xx, y = yy))
}
rys_pol <- function(testSite, kor = T, pktGran='' , korR=rKoronyModel) {
# plotowanie rozmieszczenia drzew
# atrybuty:
#kor - rysowanie zasiegu koron drzew [TRUE, FALSE]
#korR - pole z promieniami koron drzew, domyslnie pole z obliczonymi
# modelowymi wartosciami [rKoronyModel]
#pktGran - dataframe z wczytanymi pkt graniczymi dla powierzchni innych
# niz kolowki
# ustawiamy podstawowy zasięg wykresu
zakresX = range(testSite[[Xo]])
zakresY = range(testSite[[Yo]])
if (typPow == 'K') {
maxR <- testSite[[R]][1]
} else {
maxR <- max(c(zakresX, zakresY))
}
# ustalamy ramke wykresu do umieszczenia w niej danych
plot(1, xlab = 'X', ylab = 'Y',
xlim = c(zakresX[1], zakresX[2]), ylim= c(zakresY[1], zakresY[2]), type = 'n', asp=1)
if (typPow == 'K') {
okrag = list()
for (i in 1:((maxR %/%5) + 1)) {
okrag = pktObwodu(0,0,5*i)
par(new=T)
lines(x=as.vector(okrag[['x']]),
y=as.vector(okrag[['y']]),
lty = 2, asp=1,
col=312 ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]))
}
par(new=T)
# dodanie granicy powierzchni probnej (kolowki)
okrag = pktObwodu(0, 0, testSite[[R]][1])
lines(x=as.vector(okrag[['x']]),
y=as.vector(okrag[['y']]), asp=1,
col='black' ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]))
}
if (kor == T && korR %in% names(testSite)){
# rysowanie koron drzew wg podanej kolumny
# jezeli wartosc jest liczba to rysujemy korone
for (i in 1:nrow(testSite)) {
if (is.numeric(testSite[[korR]][i])) {
korona <- list()
korona <- pktObwodu(testSite[[Xo]][i], testSite[[Yo]][i], testSite[[korR]][i])
par(new=T)
# dodanie granicy powierzchni probnej (kolowki)
lines(x=as.vector(korona[['x']]),
y=as.vector(korona[['y']]),
col='green' ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]), asp=1)
}
}
}
par(new=T)
symbols(x=as.vector(testSite[[Xo]]), y=as.vector(testSite[[Yo]]),
xlim = c(zakresX[1], zakresX[2]), ylim= c(zakresY[1], zakresY[2]),
circles = as.vector(testSite[[d13]]/100),
inches = F, ann=F, asp=1 ,
bg="steelblue2",
fg='black')
najblizsze = drzCentral(testSite)
par(new=T)
symbols(x=(najblizsze[[Xo]]), y=(najblizsze[[Yo]]),
xlim = c(zakresX[1], zakresX[2]), ylim= c(zakresY[1], zakresY[2]),
circles = (najblizsze[[d13]]/100),
inches = F, ann=F, asp=1 ,
bg="red",
fg='black')
# dodajemy granice powierzchni na podstawie dodanych pktGran
if (typPow=='I') {
if (is.data.frame(pktGran) && nrow(pktGran[pktGran[[nrPow]]==testSite[[nrPow]][1],])>0){
pkt_zak <- pktGran[pktGran[[nrPow]]==testSite[[nrPow]][1], c(Xorg, Yorg)]
} else {
print('Nie podano pkt graniczych, granica oparta zostanie na zewnetrznych drzewach')
pkt_zak <- testSite[, c(Xorg, Yorg)]
}
wyb_pkt_gran <- chull(pkt_zak)
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
wyb_pkt_gran <- pkt_zak[wyb_pkt_gran, c(Xorg, Yorg)]
par(new=T)
# dodanie granicy powierzchni probnej
lines(x=as.vector(wyb_pkt_gran[[Xorg]]),
y=as.vector(wyb_pkt_gran[[Yorg]]),
col='black' ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]), asp=1)
}
text(x=as.vector(testSite[[Xo]]), y=(testSite[[Yo]]),
labels=testSite[[nrDrzewa]], asp=1,
xlim = c(zakresX[1], zakresX[2]), ylim = c(zakresY[1], zakresY[2]), pos=3, cex= 0.7)
}
rys_dh <- function(testSite) {
testSite = na.omit(testSite)
zakd13 = range(testSite[[d13]])
zakh = range(testSite[[h]])
plot(testSite[[d13]], testSite[[h]],
xlim = c(zakd13[1]-1, zakd13[2]+1), ylim= c(zakh[1]-1, zakh[2]+1),
col = 'blue', pch=19, cex=1, lty='solid', xlab="d13", ylab="h")
najblizsze = drzCentral(testSite)
par(new=T)
plot(najblizsze[[d13]], najblizsze[[h]],
xlim = c(zakd13[1]-1, zakd13[2]+1), ylim= c(zakh[1]-1, zakh[2]+1),
col = 'red', pch=19, cex=1, lty='solid', xlab='', ylab = '')
text(testSite[[d13]], testSite[[h]], labels = testSite[[nrDrzewa]],
xlim = c(zakd13[1]-1, zakd13[2]+1), ylim= c(zakh[1]-1, zakh[2]+1),
cex= 0.7, pos=3)
}
rys <- function(testSite, pktGran = '') {
par(mfcol=c(1,2))
rys_dh(testSite)
rys_pol(testSite, pktGran=pktGran)
}
G <- function(testSite, typ='A', powierzchnia, typ_pow ='m') {
# oblicza piersnicowe pole przekroju d-stanu dla roznych opcji
# typ liczenia:
# A - liczy dla wszystkich drzew, przelicza na g/ha
# M - liczy dla najgrubszego drzewa na powierzchni i przelicza na g/ha
# D - liczy dla 100 najgrubszych drzew na ha (przelicza na g/ha)
# W - liczy dla podanego drzewa i przelicz jego gi na ha
# powierzchnia - wartosc domyslnie podawana w metrach, inne wartosci [ha, km2]
# typ_pow - jednostka powierzchni w jakiej zostala przekazana do funkcji
# trigger dla powierzchni w ktorych nie sa wpisane wszystkie piersnice
dalej = 1
if (length(testSite[[d13]]) != length(na.omit(testSite[[d13]]))) {
print(testSite)
print(paste('brak wpisanych piersnic w powierzchni nr. ', testSite[[nrPow]][1], sep= ""))
print('Dla tej powierzchnie nie beda liczone statystyki')
dalej = 0
}
# przeliczenie powierzchni na odpowiednie wartosci
if (typ_pow == 'm') {
pow1<- powierzchnia / 10000 # przeliczamy powierzchnie z metrow na ha
} else if (typ_pow == 'km2'){
pow1<- powierzchnia * 100
} else if (typ_pow == 'ha'){
pow1<- powierzchnia # zostawiamy w ha
} else {
print('nierozpoznana jednostka powierzchni, dostepne sa: [m ,ha, km2]')
dalej <- 0
}
if (typ == 'D' | typ == 'A' & dalej == 1) {
#wybieramy zakres piersnic i liczymy dla nich G
if (typ == "A") {
zakd13 <- range(testSite[[d13]])
pocz <- zakd13[1] %/% 2 - 1
kon <- (zakd13[2] %/% 2) +1
ts <- testSite
} else if (typ == "D") {
ts <- testSite[order(-testSite[[d13]]),]
iledrz <- ile100( nrow(testSite), pow1)
if (nrow(ts) < iledrz) {
iledrz <- nrow(ts)
}
ts <- ts[1:iledrz, ]
zakd13 <- range(ts[[d13]])
pocz <- zakd13[1] %/% 2 - 1
kon <- (zakd13[2] %/% 2) + 1
}
g <- 0
d_n <- 0
for (i in pocz:kon) {
dd <- subset(ts,
ts[[d13]] >(i*2 ) & ts[[d13]]<= ((1+i)*2) )
# jezeli w tym przedziale nie ma zadnych piersnic idziemy dalej
if (length(dd[[d13]]>0)) {
gi <-(( (i+1)**2)*3.1415)/10000
d_n <- d_n + length(dd[[d13]])
g <- g + (gi * length(dd[[d13]]))
}
}
if (is.numeric(pow1)& pow1> 0) {
return((g/d_n)/pow1) # nie przeliczone na powierzchnie ani na nic
}
} else if (typ == 'M' & dalej == 1 ) { #"M"
ts <- testSite[order(-testSite[[d13]]),]
g <- (3.1415 * (ts[[d13]][1]/2)**2)/10000
return(g)
} else if (typ == 'W' & dalej == 1 & nrow(testSite) == 1) { # 'W'
g <- (3.1415 * (testSite[[d13]][1]/2)**2)/10000
return(g)
}
}
drzCentral <- function (testSite) {
#wyznacza drzewa centralne na powierzchni - funkcja testowa
#do modeli odleglosci wyznaczane sa w inny sposob
if (typPow == 'K') {
zakodl = range(testSite[[odlDrzew]]) #wartosc min i max odl drzew od srodka
rKol = testSite[[R]][1] # spisujemy promien kolowki do zmiennej
# sprawdzenie czy drzewo centralne spelnia minimalne warunki
if (rKol - zakodl[1] < 5) {
print('Odleglosc drzewa centralnego od granicy powierzchni jest mniejsza niz 5 m')
}
# ustalamy minimalna odleglos drzewa centralnego od granicy kolowki
if (rKol > 9.999) {
odlGran = 8
} else {
odlGran = rKol - 2
}
#wybieramy drzewa spelniajace warunek odleglosci od granicy
dcent = subset(testSite, (rKol - testSite[[odlDrzew]]) > odlGran )
} else {
dcent <- subset(testSite, odlDoGran > 8)
}
# zwracamy dataframe z drzewami spełniającymi kryteria odległosciowe
return(dcent)
}
ile100 <- function(n, pow) {
#zwraca ilosc drzew ktore nalezy wziasc z powierzchni aby wyliczyc ilosc
#drzew na ha
drz100 = 100*pow
if (drz100 < 1) {
drz100 <- 1
}
#zwracamy zaookraglona wartosc int
return(round(drz100, 0))
}
zapiszWynik <- function( wypis, nazwaPliku ) {
#zapisuje wynik w postaci pliku csv
write.table( wypis, file=nazwaPliku, col.names = T, row.names=F, sep=";")
}
wczytajPktGran <- function(f) {
# wczytuje pkt na podstawie ktorych wyznacza jest nastepnie granica pow.
data_pktgr <- read.table(f, header = T, sep=';')
return(data_pktgr)
}
wczytajPow <- function(fileIn, pktGran = data.frame()){
# kolumny niezbedne do obliczen
neededColumns <- c(nrPow, d13, h, hk, azymut, odlDrzew, gat, nrDrzewa, R,
spadek, wystawa, rKorony, Xorg, Yorg, Zorg, Xsr, Ysr)
# wczytaj dane z pliku
data<-read.table(fileIn, header = T, sep=";")
nAll <- nrow(data[1]) # liczba rekordow w pliku
nSites <- length(unique(data[nrPow])) # ile powierzchni w pliku
listNames <- unique(data[nrPow]) # lista nazw powierzchni
# usuwanie zbednych kolumn z danych wejsciowych
if (kasNiep == "T") {
data <- data[,(names(data) %in% neededColumns)]
}
Xo <- c()
Yo <- c()
rKoronyModel <- c()
# wskazniki do obliczania szerokosci korony wg modelu
a <- 0
b <- 0
for (i in 1:nAll) {
if (typPow == 'K') {
# obliczamyt wspolrzedne przy kolowkach
if (Xsr %in% names(data) & Ysr %in% names(data)){
X <- data[[Xsr]][i]
Y <- data[[Ysr]][i]
}
przel <- wspBL2XY(X, Y, data[[azymut]][i], data[[odlDrzew]][i])
Xo <- c(Xo, przel[1])
Yo <- c(Yo, przel[2])
} else {
# przepisujemy wspolrzedne orginalne jako obliczone dla zgodnosci
Xo <- c(Xo, data[[Xorg]][i])
Yo <- c(Yo, data[[Yorg]][i])
}
if (grepl(data[[d13]][i], ",")){
print("Dane liczbowe oddzielone są przecinkami, powiny być kropki")
break
}
#TODO: dodanie wspolczynnikow dla innych gatunkow
if (data[i,gat] == 'So') {
a<--0.1797
b<-0.6267
} else if(data[i,gat] == 'Jd') {
a<-0.092
b<-0.538
} else if(data[[gat]][i] == 'Św') {
a<--0.3232
b<-0.6441
} else if(data[i,gat] == 'Db') {
a<--0.3973
b<-0.7328
} else if(data[i,gat] == 'Bk') {
a<-0.1366
b<-0.6183
} else if(data[i,gat] == 'Gb') { # wartosci skopiowane z BK
a<-0.1366
b<-0.6183
} else if(data[i,gat] == 'Js') {
a<--0.1366
b<-0.7328
} else { # jezeli zadne z powyzszych to brak wspolczynnikow
a <- NA
b <- NA
}
if (!is.na(a) & data[[d13]][i] > 0 ){
rKoronyModel <- c(rKoronyModel , round((0.5 * exp(a + b* log(data[[d13]][i]))),2))
} else { # jezeli gatunek jest inny niz powyzsze - brak wspolczynnikow do obliczen
rKoronyModel <- c(rKoronyModel, 0 )
}
}
#dodanie obliczonych koordynatowa do orginalnej tabeli
data <- cbind(data, Xo, Yo)
# dodanie modelowych wielkosci korony
data <- cbind(data, rKoronyModel)
Xo = 'Xo' # nazwa kolumny z obliczonymi wsp X dla kazdej powierzchni
Yo = 'Yo' # nazwa kolumny z obliczonymi wsp Y dla kazdej powierzchni
rKoronyModel <- 'rKoronyModel' # nazwa kolumny z wyliczonym promieniem korony
b = list()
for (i in (1:length(listNames[,1]))) {
if (sapply(subset(data, data[nrPow]== listNames[i,1]), is.numeric)[[d13]] == T){
t1 = subset(data, data[nrPow] == listNames[i,1])
# ustawiamy nazwy wierszy na bardziej sensowne
rownames(t1) <- paste(t1[[nrPow]], t1[[nrDrzewa]], sep='_')
# oblicz odleglosc drzew do granicy powierzchni
if(typPow == 'K'){
# czy kolowka ma podane orginalne wspolrzedne srodka
if (Xorg %in% names(t1) & Yorg %in% names(t1)){
Xt <- t1[[Xsr]][1]
Yt <- t1[[Ysr]][1]
} else { # jesli nie ustawiamy na 0
Xt <- 0
Yt <- 0
}
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], pktObwodu(Xt, Yt, t1[[R]][1]))
} else {
if (is.data.frame(pktGran)== T){
if (nrow(pktGran[pktGran[[nrPow]]==t1[[nrPow]][1], ]) > 0){
pkt_zak <- pktGran[pktGran[[nrPow]]==t1[[nrPow]][1], c(Xorg, Yorg)]
wyb_pkt_gran <- chull(pktGran[pktGran[[nrPow]]==t1[[nrPow]][1], c(Xorg, Yorg)])
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
#pkt_zak[wyb_pkt_gran, c('X', 'Y')]
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], pkt_zak[wyb_pkt_gran, c(Xorg, Yorg)])
} else {
print('Brak danych dla pkt pomiarowych granicy powierzchni')
print('Obliczam odleglosci do granicy wyznaczonej na podstawie drzew ')
wyb_pkt_gran <- chull(t1[c(Xo, Yo)])
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
#pkt_zak[wyb_pkt_gran, c(Xo, Yo)]
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], t1[wyb_pkt_gran, c(Xo, Yo)])
}
} else {
print('Brak danych dla pkt pomiarowych granicy powierzchni')
print('Obliczam odleglosci do granicy wyznaczonej na podstawie drzew ')
wyb_pkt_gran <- chull(t1[c(Xo, Yo)])
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
#pkt_zak[wyb_pkt_gran, c(Xo, Yo)]
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], t1[wyb_pkt_gran, c(Xo, Yo)])
}
}
#obliczamy korekte wysokosci drzewa w zaleznosci od spadku na kolowce
if(typPow == 'K' & spadek %in% names(t1) & wystawa %in% names(t1)){
poprH <- oblPoprWys(t1, t1[[wystawa]][1], t1[[spadek]][1])
} else if (typPow == 'I' & Zorg %in% names(t1)){
roz <- range(t1[[Zorg]])
srodekH <- (roz[1] + roz[2]) / 2
# obliczamy przesuniecia wysokosciowe drzew
#transform(t1[[Zorg]], odlDoGran = Zorg - srodekH)
poprH <- sapply(t1[[Zorg]], function(x){return(round(x-srodekH, 2))})
} else {
#jezeli nie ma odpowiednich danych wstawiamy 0
poprH <- c(rep.int(0, nrow(t1)))
}
b[[toString(listNames[i,1])]] <- cbind(t1, odlDoGran, poprH)
} else {
print(paste('Brak wszystkich piersnic - pow. nr.: ', toString(listNames[i,1]) ))
}
}
odlDoGran <- 'odlDoGran'
return(b)
}
obliczCzWsp <- function(r, R, d) {
# oblicza czesc wspolna dwoch przecinajacych sie okregow
if (r>R) {
# zamieniamy promienie miejsacami
tt <- r
r <- R
R <- tt
} else if (r == R) {
R <- R + 0.0001
}
if (r+d < R) {
# jezeli drzewo jest calkowicie pod drzewem centralnym to wchodzi cala pow
wynik <- pi * r*r
} else {
# jezeli drzewa nachodza na siebie tylko w malej czesci to liczymy
cz1 = r**2*acos((d**2+r*r-R*R)/(2*d*r))
cz2 <- R*R*acos((d*d+R*R-r*r)/(2*d*R))
cz3 <- 0.5*sqrt((-d+r+R)*(d+r-R)*(d-r+R)*(d+r+R))
wynik <- cz1 + cz2 -cz3
if (is.nan(cz1) & is.nan(cz2) & is.nan(cz3)){
print(paste(r, R, d))
}
# aleternatywny sposob liczenia - zostawiony na wypadek przyszlych modyfikacji
#x1<-(d**2-r*r+R*R)/(2*d)
#print(paste('x1= ',x1))
#x2<-abs(d-x1)
#print(paste('x2= ', x2))
#y<- sqrt(R*R- (x1*x1))
#print(paste('y= ', y))
#a1<- (R**2)*acos(x1/R) - x1 * y
#print(paste('a1= ', a1))
#a2<- r**2 * acos(x2/r) - x2 * y
#print(paste('a2= ', a2))
#if (x1>d){
#a2<- pi * r**2 - a2
#print(paste('a2_podm = ',a2))
#}
#area <- a1+a2
#print(paste('pow= ', area))
}
return(wynik)
}
oblOdlGran <- function (X, Y, pktGran) {
# zwraca liste z odleglosciami kolejnych drzew od granicy powierzchni
# pktKran - dwuwymiarowa maciez z pkt opisujacymi granice powierzchni
drz <- cbind( X, Y)
odlGran <- apply(distppll(drz, cbind(pktGran[-nrow(pktGran), ], pktGran[-1,])),1,min)
return(round(odlGran, 2))
}
oblPoprWys <- function(testSite, wystawa = 0, nachyl = 0) {
# oblicza wysokosci odziomkow na powierzchni kolowej na podstawie nachylenia
# i wystawy podanych w stopniach
if (Xsr %in% names(testSite) & Ysr %in% names(testSite)){
xpocz <- testSite[[Xsr]][1]
ypocz <- testSite[[Ysr]][1]
} else {
xpocz <- 0
ypocz <- 0
}
x_zakres <- range(testSite[[Xo]])
y_zakres <- range(testSite[[Yo]])
if (wystawa <= 90 ) {
pkt1 <- c(x_zakres[1], y_zakres[1])
} else if (wystawa > 90 & wystawa <= 180) {
pkt1 <- c(x_zakres[1], y_zakres[2])
} else if (wystawa >180 & wystawa <= 270) {
pkt1 <- c(x_zakres[2], y_zakres[2])
} else if (wystawa > 270 ){
pkt1 <- c( x_zakres[2], y_zakres[1])
}
pkt2 <- wspBL2XY(pkt1[1], pkt1[2], wystawa, 1500)
# obliczamy dlugosc wektora
odlPP <- odlAB(pkt1[1],pkt1[2], pkt2[1], pkt2[2])
# obliczamy spadek dla srodka kolowki
u0 <- (( xpocz - pkt1[1]) * (pkt2[1] - pkt1[1]) + (ypocz - pkt1[2]) *(pkt2[2] - pkt1[2])) /
((pkt1[1] - pkt2[1])**2 + (pkt1[2] - pkt2[2])**2 )
spadek0 <- round(tan(d2r(nachyl)) * u0 * odlPP, 2)
#obliczamy spadki
ulist <- c()
for (i in 1:nrow(testSite)) {
u <- ((testSite[[Xo]][i] - pkt1[1]) * (pkt2[1] - pkt1[1]) + (testSite[[Yo]][i] - pkt1[2]) *(pkt2[2] - pkt1[2])) /
((pkt1[1] - pkt2[1])**2 + (pkt1[2] - pkt2[2])**2 )
#obliczamy poprawke wysokosci drzewa na spadek i ustawiamy 0 na srodku kolowki
spadek <- round((tan(d2r(nachyl)) * u * odlPP) - spadek0, 2)
ulist <- rbind(ulist, spadek)
}
return(ulist)
}
modeleCI <- function(testSite, debug = 1) {
#dla podanej powierzchni zwraca obliczenia dla drzew centralnych wyzanczonych
#na podstawie wybranych modeli odleglosci.
#Dostepne modele:
#odleglosciowe: 3m, 4m, 5m, 6m, 7m, 8m, 9m, 10,
#odleglosci zaleznej od piersnic: S8, S6
#odl. zaleznej od korony: 25s, 35pr
#odl. zalezne od wysokosci drzewa: 02H, 03H, 04H, 05H
#odl. zalezne od d13: 015D, 020D, 025D, 030D
#sr. odl. Oa
#4 najblizsze
#odwroconego stozka, zalezne od a i kata rozwarcia stozka
#modele stalego r od 3 do 10 metrow
testSite$lacz<-paste(testSite[[nrPow]], testSite[[nrDrzewa]], sep='_')
wyn <- testSite
# deklaracja nazw dla poszczegolnych tablic konkurencji i modeli
konk_list_names <- c('s8_konk', 's6_konk', 's25_konk', 'pr35_konk', 'sm25_konk', 'prm35_konk',
'H02_konk', 'H03_konk', 'H04_konk', 'H05_konk', 'D015_konk',
'D020_konk', 'D025_konk', 'D030_konk', 'n2_konk', 'n4_konk', 's_konk', 'sm_konk')
st_list <- c('0150', '0160', '0170', '0450', '0460', '0470', '0650',
'0660', '0670', '06650', '06660', '06670', 'nk50', 'nk60', 'nk70')
model_list_names <- list('s8','s6','s25','pr35','sm25','prm35',
'H02','H03','H04','H05','D015','D020','D025',
'D030','n2','n4', 's', 'sm')
for (j in 3:10) {
dcent <- subset(testSite, testSite[['odlDoGran']] > j)
odl <- j
if (nrow(dcent)>0){
mod <- wskCI(dcent, testSite, odl)
names(mod) <- paste(paste('m', toString(j),sep=""), names(mod), sep='_')
mod$lacz <- rownames(mod)
print(paste('Obliczam model o stalym promieniu dla:', j, 'm', 'drzew cent:', nrow(mod) ))
wyn <- merge(mod, wyn, all.y=T)
}
}
# modele oparte o piersnice
dcent <- subset(testSite, testSite[['odlDoGran']] > 7 & h>0)
if (nrow(dcent) == 0){
print(paste('Pow.',testSite[[nrPow]][1] , '- Brak drzew spelniajacych kryteria drzewa centralnego', sep=' '))
return(data.frame())
} else {
print(paste('Do analizy wybrano:', nrow(dcent), '/', nrow(testSite), 'drzew'))
#zerowanie wartosci modeli dla kazdego drzewa
mod_s8 <- data.frame()
mod_s6 <- data.frame()
mod_s25 <- data.frame()
mod_pr35 <- data.frame()
mod_sm25 <- data.frame()
mod_prm35 <- data.frame()
mod_H02 <- data.frame()
mod_H03 <- data.frame()
mod_H04 <- data.frame()
mod_H05 <- data.frame()
mod_D015 <- data.frame()
mod_D020 <- data.frame()
mod_D025 <- data.frame()
mod_D030 <- data.frame()
mod_n2 <- data.frame()
mod_n4 <- data.frame()
mod_s <- data.frame()
mod_sm <- data.frame()
mod_st <- list()
# grupujemy modele w liste
model_list <- list(mod_s8,mod_s6,mod_s25,mod_pr35,mod_sm25,mod_prm35,
mod_H02,mod_H03,mod_H04,mod_H05,mod_D015,mod_D020,mod_D025,
mod_D030,mod_n2,mod_n4, mod_s, mod_sm)
#print(dcent)
for (j in 1:nrow(dcent)) {
# obliczenia dla kazdego drzewa centralnego po kolei
print(' ')
print(paste('Obliczam drzewo centralne nr:', j, '/', nrow(dcent), '(', dcent[j,nrDrzewa], ')'))
xmin <- dcent[[Xo]][j]-40
xmax <- dcent[[Xo]][j]+40
ymin <- dcent[[Yo]][j]-40
ymax <- dcent[[Yo]][j]+40
# pobieramy wycinek powierzchnie o boku 40m do obliczen ze stozkiem
wycinek20 <- subset(testSite, Xo< xmax & Xo> xmin & Yo< ymax & Yo> ymin)
# obliczamy odleglosci do tego drzewa centralnego
odltemp <- c()
for (ii in 1:nrow(wycinek20)) {
odlobl <- odlAB(dcent[[Xo]][j],
dcent[[Yo]][j],
wycinek20[[Xo]][ii],
wycinek20[[Yo]][ii])
odltemp <- c(odltemp,round(odlobl, 2) )
}
# dopisujemy odleglosci do wycinka i sortujemy od najmniejszej odleglosci
wycinek20t <- cbind(wycinek20, odltemp)
odltemp<-'odltemp'
wycinek20 <- wycinek20t[order(wycinek20t[[odltemp]]),]
Oa <- min(wycinek20[[odltemp]][2])*2
#print(wycinek20)
# zerowanie grup konkurencyjnych dla poszczegolnych drzew
s6_konk <- data.frame()
s8_konk <- data.frame()
s25_konk <- data.frame()
sm25_konk <- data.frame()
pr35_konk <- data.frame()
prm35_konk <- data.frame()
H02_konk<- data.frame()
H03_konk<- data.frame()
H04_konk<- data.frame()
H05_konk<- data.frame()
D015_konk<- data.frame()
D020_konk<- data.frame()
D025_konk<- data.frame()
D030_konk<- data.frame()
n2_konk <- data.frame()
n4_konk <- data.frame()
s_konk <- data.frame()
sm_konk <- data.frame()
for (jj in 1:nrow(wycinek20)) {
# obliczamy odleglosc do drzewa konkurencyjnego
odlRobocza <- odlAB(dcent[[Xo]][j],
dcent[[Yo]][j],
wycinek20[[Xo]][jj],
wycinek20[[Yo]][jj])
# obliczamy grupe konkurencyjna dla modelu s6
if ((dcent[[d13]][j] + wycinek20[[d13]][jj]) / 6 > odlRobocza) {
s6_konk <- rbind(s6_konk, wycinek20[jj, ])
}
# obliczamy grupe konkurencyjna dla modelu s8
if ((dcent[[d13]][j] + wycinek20[[d13]][jj]) / 8 > odlRobocza) {
s8_konk <- rbind(s8_konk, wycinek20[jj, ])
}
# czy drzewo posiada pomiary szerokosci korony
if (rKorony %in% names(wycinek20) & is.numeric(wycinek20[[rKorony]][jj])) {
if ( wycinek20[[rKorony]][jj] > 0){
if (odlRobocza < 2.5 * 2 * dcent[[rKorony]][j]) {
s25_konk <- rbind(s25_konk, wycinek20[jj, ])
}
#obliczamy grupe konkurencyjna dla modelu pr35
if (odlRobocza < 3.5 * dcent[[rKorony]][j]) {
pr35_konk <- rbind(pr35_konk, wycinek20[jj, ])
}
# obliczamy grupe konkurencyjna dla modelu s
if (odlRobocza < dcent[[rKorony]][j] + wycinek20[[rKorony]][jj] &
!is.na(dcent[[hk]][j]) & !is.na(wycinek20[[h]][jj]) & !is.na(wycinek20[[hk]][jj]) &
is.numeric(wycinek20[[h]][jj]) & is.numeric(wycinek20[[hk]][jj]) &
is.numeric(dcent[[h]][j])){
if ((dcent[[hk]][j]) <=
(wycinek20[[h]][jj] - dcent[['poprH']][j] + wycinek20[['poprH']][jj])){
s_konk <- rbind(s_konk, wycinek20[jj,])
}
}
}
}
# czy drzewo posiada modelowe szerokosci korony
if (rKoronyModel %in% names(wycinek20) & wycinek20[[rKoronyModel]][jj] > 0 &
is.numeric(wycinek20[[rKoronyModel]][jj])) {
if (odlRobocza < 2.5 * 2 * dcent[[rKoronyModel]][j]) {
sm25_konk <- rbind(sm25_konk, wycinek20[jj, ])
}
#obliczamy grupe konkurencyjna dla modelu m25s
if (odlRobocza < 3.5 * dcent[[rKoronyModel]][j]) {
prm35_konk <- rbind(prm35_konk, wycinek20[jj, ])
}
# obliczamy grupe konkurencyjna dla modelu sm
if (odlRobocza < dcent[[rKoronyModel]][j] + wycinek20[[rKoronyModel]][jj] &
!is.na(dcent[[hk]][j]) & !is.na(wycinek20[[h]][jj]) & !is.na(wycinek20[[hk]][jj]) &
is.numeric(wycinek20[[h]][jj]) & is.numeric(wycinek20[[hk]][jj]) &
is.numeric(dcent[[h]][j])){
if ((dcent[[hk]][j]) <=
(wycinek20[[h]][jj] - dcent[['poprH']][j] + wycinek20[['poprH']][jj])){
sm_konk <- rbind(sm_konk, wycinek20[jj,])
}
}
}
# obliczamy grupe konkurencyjna w zaleznosci od wysokosci drzewa centralnego
if (!is.na(dcent[[h]][j]) & is.numeric(dcent[[h]][j]) &
!is.na(wycinek20[[h]][j]) & is.numeric(wycinek20[[h]][j]) & wycinek20[[h]][j]>0) {
if (odlRobocza < dcent[[h]][j] * 0.2) {
H02_konk <- rbind(H02_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.3) {
H03_konk <- rbind(H03_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.4) {
H04_konk <- rbind(H04_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.5) {
H05_konk <- rbind(H05_konk, wycinek20[jj,])
}
}
# grupy konkurencyjne oparte na piersnicach
if (odlRobocza < dcent[[h]][j] * 0.15) {
D015_konk <- rbind(D015_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.20) {
D020_konk <- rbind(D020_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.25) {
D025_konk <- rbind(D025_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.30) {
D030_konk <- rbind(D030_konk, wycinek20[jj,])
}
# model sredniej odleglosci Oa
if (Oa > 0) {
if (odlRobocza < Oa) {
n2_konk <- rbind(n2_konk, wycinek20[jj,])
}
}
}
#wybieramy grupe konkurencyjna dla modelu 4n
ile_wyb <- 1
if (nrow(wycinek20) >= 5 & hk %in% names(wycinek20) & h %in% names(wycinek20)) {
for (l in 2:nrow(wycinek20)) {
if (is.numeric(wycinek20[[h]][l]) & is.numeric(wycinek20[[hk]][l])&
is.numeric(dcent[[h]][j]) & !is.na(dcent[[hk]][j])&
!is.na(wycinek20[[h]][l]) & !is.na(wycinek20[[hk]][l])) {
if ((dcent[[hk]][j]) <=
(wycinek20[[h]][l] - dcent[['poprH']][j] + wycinek20[['poprH']][l]) & ile_wyb <5){
n4_konk <- rbind(n4_konk, wycinek20[l,])
}
# jezeli drzewo nie weszlo to i tak je liczymy
ile_wyb <- ile_wyb + 1
}
}
}
# konkurencja dla modeli zwiazanych ze stozkiem
wys_stozka <- c(0.1, 0.4, 0.6, 0.66)
# obliczamy tangens kata alfa
tan50 <- tan(d2r(90-(50/2)))
tan60 <- tan(d2r(90-(60/2)))
tan70 <- tan(d2r(90-(70/2)))
lista_k <- c(tan50, tan60, tan70)
st_konk <- list()
nk_trig <- -999 # trigger dla nasady korony
if (h %in% names(dcent) & hk %in% names(dcent)){
for (l in 1:nrow(wycinek20)) {
if (is.numeric(wycinek20[[h]][l]) & is.numeric(wycinek20[[hk]][l])&
is.numeric(dcent[[h]][j]) & !is.na(dcent[[hk]][j])&
!is.na(wycinek20[[h]][l]) & !is.na(wycinek20[[hk]][l])) {
for (ll in 1:length(lista_k)){
for (kk in 1:length(wys_stozka)) {
# skladamy nazwe grupy konkurencyjnej i dodajemy pod nia dane konkurencji
naz_temp <- paste(0,
unlist(strsplit(toString(wys_stozka[kk]),
split='.', fixed=T))[2],
toString((ll+4)*10), sep="")
#sprawdzamy czy warunek jest dodatni lub zerowy
if (((wycinek20[[h]][l] + wycinek20[['poprH']][l] - dcent[['poprH']][j]) -
(wys_stozka[kk]* (dcent[[h]][j])))/
lista_k[ll] >0){
war_wys <- ( ((wycinek20[[h]][l] + wycinek20[['poprH']][l]- dcent[['poprH']][j]) -
(wys_stozka[kk]* (dcent[[h]][j])))/
lista_k[ll])
} else {
war_wys <- 0.0000000000000001
}
if (wycinek20[[odltemp]][l] < war_wys & wycinek20[[odltemp]][l] >0 ) {
st_konk[[naz_temp]] <- rbind(st_konk[[naz_temp]], wycinek20[l,])
}
# obliczenia dla nasady korony
nazkor_temp <- paste('nk', toString((ll+4)*10), sep="")
if (((wycinek20[[h]][l] - dcent[['poprH']][j]+ wycinek20[['poprH']][l] )-
(dcent[[hk]][j]))>0 ){
war_wys_nk <- (((wycinek20[[h]][l] - dcent[['poprH']][j]+ wycinek20[['poprH']][l] )-
(dcent[[hk]][j]))/
lista_k[ll])
} else {
war_wys_nk <- 0.00000000000001
}
#print(paste(nazkor_temp, wycinek20[[odltemp]][l], war_wys_nk))
if (wycinek20[[odltemp]][l] < war_wys_nk & wycinek20[[odltemp]][l] >0 &
nk_trig != ll) {
st_konk[[nazkor_temp]] <- rbind(st_konk[[nazkor_temp]],wycinek20[l,])
# trigger pomijajacy jedna petlekk
nk_trig <- ll
}
}
}
}
}
}
# obliczenia dla standardowych modeli
konk_list <- list(s8_konk, s6_konk, s25_konk, pr35_konk, sm25_konk, prm35_konk,
H02_konk, H03_konk, H04_konk, H05_konk, D015_konk,
D020_konk, D025_konk, D030_konk, n2_konk, n4_konk, s_konk, sm_konk)
names(konk_list) <- konk_list_names
for (ii in 1:length(konk_list)) {
if (is.data.frame(konk_list[[konk_list_names[ii]]]) && !nrow(konk_list[[konk_list_names[ii]]])==0) {
model_list[[ii]] <- rbind(model_list[[ii]], wskCI(dcent[j,], konk_list[[konk_list_names[ii]]], odlMod=range(konk_list[[konk_list_names[ii]]][['odltemp']])[2]))
if (debug == 1) {
print(paste(' Model ', konk_list_names[ii], ', drzew konk: ', nrow(konk_list[[konk_list_names[ii]]]), ', max odl: ', range(konk_list[[konk_list_names[ii]]][['odltemp']])[2], sep=''))
}
}
}
# obliczenia dla modeli ze stozkiem
for (item in st_list) {
if ((is.data.frame(st_konk[[item]]))&& !nrow(st_konk[[item]])==0) {
mod_st[[item]] <- rbind(rbind(mod_st[[item]], wskCI(dcent[j,], st_konk[[item]], odlMod=range(st_konk[[item]][['odltemp']])[2])))
if (debug == 1) {
print(paste(' Model ', item , ', drzew konk: ', nrow(st_konk[[item]]), ', max odl: ', range(st_konk[[item]][['odltemp']])[2], sep=''))
}
}
}
}
for (model in 1:length(model_list)) {
if (nrow(model_list[[model]]) >0) {
names(model_list[[model]]) <- paste(model_list_names[[model]], names(model_list[[model]]), sep='_')
model_list[[model]]$lacz <- rownames(model_list[[model]])
wyn <- merge(model_list[[model]], wyn, all.y=T, all.x=T)
}
}
for (model in names(mod_st)) {
if (nrow(mod_st[[model]]) > 0) {
names(mod_st[[model]]) <- paste(model, names(mod_st[[ model ]]), sep='_')
mod_st[[model]]$lacz <- rownames(mod_st[[model]])
wyn <- merge(mod_st[[model]], wyn, all.y=T, all.x=T)
}
}
# zamieniamy wszystkie NA na 0
wyn[is.na(wyn)] <- 0
# zwracamy ostateczna tabele
return(wyn)
}
}
wskCI <- function(dcent, testSite, odlMod = 6) {
# funkcja oblicza wskazniki konkurencjii w dla drzew centralnych zapisanych jako
#dcent (data.frame) spelniajacych kryteria odleglosci od granicy powierzchni.
#Wskazniki oblczane sa na podstawie drzew konkurencyjnych testSite (data.frame)
#ktore wystepuja w podanej odleglosci od drzewa centralnego - odlMod (numeric)
# wyznaczamy drzewa centralne - dla celow testowych
if (!is.data.frame(dcent)){
dcent <- drzCentral(testSite)
}
# lista drzew bez wpisanej wysokosci, na koncu wypiszemy komunikat
brakh <-c()
if (nrow(dcent) > 0) {
ll_names <- list('NrDrz', 'NrPow', 'CI','CI2' ,'CI3','CI4' ,'CI5' ,'CI6' ,'CI7','CI8' ,
'CI9','CI10','CI11','CI12','CI13','CI14','CI15','CI16','CI17',
'CI18','CI19','CI20','CI21','CI22','CI23','CI24','CI25','CI26',
'AP', 'CI28', 'CI29', 'BAL','BALMOD', 'CI32', 'CI33',
'SDI', 'CI35', 'CI36', 'CI37',
'ME3mod', 'BR1mod', 'BR2mod', 'BR3mod', 'BR4mod', 'BR5mod', 'PK16mod', 'CCmod',
'ME3rzecz', 'BR1rzecz', 'BR2rzecz', 'BR3rzecz', 'BR4rzecz', 'BR5rzecz',
'PK16rzecz', 'CCrzecz', 'Konkurenci' )
ll<- data.frame()
for (i in 1:nrow(dcent)) {
CI <- 0
CI2 <- 0
CI3 <- 0
CI4 <- 0
CI5 <- 0
CI6 <- 0
CI7 <- 0
CI8 <- 0
CI9 <- 0
CI10 <- 0
CI11 <- 0
CI12 <- 0
CI13 <- 0
CI14 <- 0
CI15 <- 0
CI16 <- 0
CI17 <- 0
CI18 <- 0
CI19 <- 0
CI20 <- 0
CI21 <- 0
CI22 <- 0
CI23 <- 0
CI24 <- 0
CI25 <- 0
CI26 <- 0
BAL <- 0
sCI35 <- 0
nCI35 <- 0
CI35 <- 0
CI36 <- 0
gCI37 <- 0
CI37 <- 0
CI32 <- 0
CI33 <- 0
Oirzecz <- 0
Oimodel <- 0
di <- 0
ME3rzecz <- 0
BRskl <- 0
BR1rzecz <- 0
BR2rzecz <- 0
BR3rzecz <- 0
BR4rzecz <- 0
BR5rzecz <- 0
CCrzecz <- 0
ME3mod <- 0
BRskl <- 0
AP <- 0
CI28 <- 0
CI29 <- 0
SDI <- 0
BALMOD <- 0
BR1mod <- 0
BR2mod <- 0
BR3mod <- 0
BR4mod <- 0
BR5mod <- 0
CCmod <- 0
PK16mod <- 0
PK16rzecz <- 0
SDIsklad <- 0
SDItrig <- 0
ilekonkurentow <- 0
if (odlMod == 999) {
pr_pow <- 40
} else {
pr_pow <- odlMod
}
# obliczamy powierzchnie wybranego wycinka pow. probnej dla drzewa centralego
pow <- 3.1415*pr_pow**2
# czy w danych wystepuja rzeczywiste pomiary korony drzewa
koronaLiczenie <- 1
if (rKorony %in% colnames(testSite)) {
koronaLiczenie <- 2
}
for (j in 1:nrow(testSite)) {
# odleglosi miedzy drzewami w celu wyznaczeniea drzew konkurencyjnych
# w zasiegu analizowanego drzewa
odl<-odlAB(dcent[[Xo]][i],
dcent[[Yo]][i],
testSite[[Xo]][j],
testSite[[Yo]][j])
# promienie konkurencji
#if (odl < 3) {r <- 3} else {r <- as.integer(odl)}
r <- pr_pow
if (odl > 0 & odl < pr_pow ) {
ilekonkurentow <- ilekonkurentow + 1
# podzbior do obliczenia SDI
if (SDItrig == 0) {
SDIsklad <- testSite[j,]
SDItrig <- 1
} else {
SDIsklad <- rbind(SDIsklad, testSite[j,])
}
#piersnice wieksze od d centralnego i modyfikacje
if (testSite[[d13]][j] > dcent[[d13]][i]) {
BAL <- BAL + G(testSite[j,], powierzchnia=pow, typ="W")
CI11 <- CI11 + atan(testSite[[d13]][j]/odl)**-1
CI25 <- CI25 + ( testSite[[d13]][j]-dcent[[d13]][i] ) / odl
CI36 <- CI36 + (testSite[[d13]][j] / dcent[[d13]][i])**0.5
} else if (dcent[[d13]][i] == testSite[[d13]][j]){
CI36 <- CI36 + (testSite[[d13]][j] / dcent[[d13]][i])**1
} else if (dcent[[d13]][i] > testSite[[d13]][j]){
CI36 <- CI36 + (testSite[[d13]][j] / dcent[[d13]][i])**2
}
#wysokosci - wieksze i wszystkie inne o ile istnieja
if (!is.na(testSite[[h]][j]) & !is.na(dcent[[h]][i])) {
if (testSite[[h]][j] > dcent[[h]][i]) {
CI33 <- CI33 + G(testSite[j,], powierzchnia=pow, typ='W')
CI12 <- CI12 + atan(testSite[[d13]][j]/odl)**-1
CI15 <- CI15 + atan(testSite[[h]][j]/odl)**-1
CI19 <- CI19 + atan(abs((testSite[[h]][j]-dcent[[h]][i])/odl))**-1
}
if (testSite[[h]][j]>dcent[[h]][i]/2) {
CI14 <- CI14 + atan(testSite[[h]][j]/odl)**-1
}
if (testSite[[h]][j]>dcent[[h]][i]*0.8) {
CI18 <- CI18 + atan(abs((testSite[[h]][j] - 0.8*dcent[[h]][i])/odl))**-1
}
CI2 <- CI2 + testSite[[h]][j]/(dcent[[h]][i]*odl)
CI13 <- CI13 + atan(testSite[[h]][j]/odl)**(-1)
CI17 <- CI17 + (testSite[[h]][j]/dcent[[h]][i])*atan(d2r(testSite[[h]][j]/odl))**(-1)
} else { # jak nie istnieja to pomijamy i wypisujemy
if (!is.na(testSite[[h]][j])) {
if (!(testSite[[nrDrzewa]][j] %in% brakh)) {
brakh <- c(brakh, testSite[[nrDrzewa]][j] )
}
} else{
if (!(dcent[[nrDrzewa]][i] %in% brakh)) {
brakh <- c(brakh, dcent[[nrDrzewa]][i] )
}
}
}
di <- di + testSite[[d13]][j]
# obliczenia zalezne tylko i wylacznie od odleglosci
CI <- CI + testSite[[d13]][j] / (dcent[[d13]][i] * odl)
CI3 <- CI3 + testSite[[d13]][j]/(dcent[[d13]][i]*(odl+1))
CI4 <- CI4 + abs((testSite[[d13]][j]/dcent[[d13]][i])*exp(-(16*odl)/(testSite[[d13]][j]+dcent[[d13]][i])))
CI5 <- CI5 + G(testSite[j,],powierzchnia=pow, typ='W')/(odl/r)
CI6 <- CI6 + (G(testSite[j,],powierzchnia=pow, typ='W')/
G(dcent[i,],powierzchnia=pow, typ='W'))/(odl/r)
CI7 <- CI7 + (testSite[[d13]][j]/dcent[[d13]][i])/(odl/r)
CI8 <- CI8 + (testSite[[d13]][j]/dcent[[d13]][i])/sqrt(odl/r)
CI9 <- CI9 + (testSite[[d13]][j]/dcent[[d13]][i])/(odl/r)**2
CI10 <- CI10 + atan(testSite[[d13]][j]/odl)**(-1)
CI16 <- CI16 + (testSite[[d13]][j]/dcent[[d13]][i])*atan(testSite[[d13]][j]/odl)**(-1)
CI20 <- CI20 + (testSite[[d13]][j] / odl)
CI21 <- CI21 + ( testSite[[d13]][j]/dcent[[d13]][i] ) / odl**2
CI22 <- CI22 + ( testSite[[d13]][j]/dcent[[d13]][i] )**2 / odl
CI23 <- CI23 + odl
CI24 <- CI24 + testSite[[d13]][j] / odl
#niezalezne od odleglosci - teoretycznie
CI32 <- CI32 + G(testSite[j,],powierzchnia=pow, typ='W')
sCI35 <- sCI35 + G(testSite[j,],powierzchnia=pow, typ='W')
nCI35 <- nCI35 + 1
gCI37 <- gCI37 + G(testSite[j,],powierzchnia=pow, typ='W')
# wskazniki zwiazane z korona drzewa
if (odl< ( testSite[[rKoronyModel]][j]+dcent[[rKoronyModel]][i] ) &
dcent[[rKoronyModel]][i] >0 & testSite[[rKoronyModel]][j] > 0 ) {
Oi <- obliczCzWsp(dcent[[rKoronyModel]][i], testSite[[rKoronyModel]][j], odl)
#print(paste( testSite[[rKoronyModel]][j]+dcent[[rKoronyModel]][i], odl ))
Oimodel <- Oimodel + (Oi/(3.1415*dcent[[rKoronyModel]][i]**2))
ME3mod <- ME3mod + (Oi/(3.1415*dcent[[rKoronyModel]][i]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BRsklm <- (Oi/(3.1415*dcent[[rKoronyModel]][i]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BR1mod <- BR1mod + BRsklm**1
BR2mod <- BR2mod + BRsklm**1.5
BR3mod <- BR3mod + BRsklm**2
BR4mod <- BR4mod + BRsklm**2.5
BR5mod <- BR5mod + BRsklm**3
CCmod <- CCmod +(3.1415*testSite[[rKoronyModel]][j]**2)
#print(CCmod)
} else if (koronaLiczenie == 2 ){
if (odl< ( testSite[[rKorony]][j]+dcent[[rKorony]][i]) &
!(is.na(testSite[[rKorony]][j])) &
!(is.na(dcent[[rKorony]][i]) ) ) {
Oir <- obliczCzWsp(dcent[[rKorony]][i], testSite[[rKorony]][j], odl)
Oirzecz <- Oirzecz + (Oir/(3.1415*dcent[[rKorony]][i]**2))
ME3rzecz <- ME3rzecz + (Oir/(3.1415*testSite[[rKorony]][j]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BRsklr <- (Oir/(3.1415*dcent[[rKorony]][i]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BR1rzecz <- BR1rzecz + BRsklr**1
BR2rzecz <- BR2rzecz + BRsklr**1.5
BR3rzecz <- BR3rzecz + BRsklr**2
BR4rzecz <- BR4rzecz + BRsklr**2.5
BR5rzecz <- BR5rzecz + BRsklr**3
CCrzecz <- CCrzecz +(3.1415*testSite[[rKorony]][j]**2)
}
}
}
}
if (ilekonkurentow > 0) {
CI26 <- 0
# specjalna petla dla wskaznika CI26
for (j in 1:nrow(testSite)) {
odl=odlAB(dcent[[Xo]][i],
dcent[[Yo]][i],
testSite[[Xo]][j],
testSite[[Yo]][j])
if (odl > 0 & odl < pr_pow ) {
CI26 <- CI26 + 3.1415*((odl*dcent[[d13]][i])/
(dcent[[d13]][i]+testSite[[d13]][j])**2*
((testSite[[d13]][j]/odl)/(di/CI23)))
}
}
#CI36 obliczane jest w petli - PAMIETAC O DODANIU DO WYPISU
CI35 <- G(dcent[i,],powierzchnia=pow, typ='W') / (sCI35/nCI35 )
CI37 <- G(dcent[i,],powierzchnia=pow, typ='W') / (gCI37)
# sorutujemy drzewa po wysokosci na stanowisku
Hwyb <- testSite[order(-SDIsklad[[h]]),]
if (nrow(SDIsklad) < ile100( nrow(SDIsklad), pow/10000)) {
iledrz <- nrow(SDIsklad)
} else {
iledrz <- ile100( nrow(SDIsklad), pow/10000)
}
# liczymy wysokosc z drzew dominujacych
H100 <- mean(Hwyb[1:iledrz,h])
RS <- sqrt(pow/nrow(SDIsklad))/H100
BALMOD <- (1-(1-BAL/G(SDIsklad, powierzchnia=pow)))/RS
# obliczenie kwadratu przecietnej piersnicy d-stanu
SDIsklad <- as.data.frame(SDIsklad)
if (nrow(SDIsklad) > 0) {
zakd13 <- range(SDIsklad[[d13]])
pocz <- zakd13[1] %/% 2
kon <- (zakd13[2] %/% 2) +1
d_n <- 0
dlicz<- 0
# kwadrat przecietnej piersnicy
for (k in pocz:kon) {
dd <- subset(SDIsklad,
SDIsklad[[d13]] > (k*2 - 0.01) & SDIsklad[[d13]]< ((1+k)*2))
# jezeli w tym przedziale nie ma zadnych piersnic idziemy dalej
if (length(dd[[d13]]>0)) {
dlicz <- dlicz + length(dd[[d13]])*(1 + k * 2)
d_n <- d_n + length(dd[[d13]])
}
}
dg <- (dlicz/d_n) # kwadrat przecietnej piersnicy d-stanu
SDI <- nrow(SDIsklad)* (25 / dg)**( -1.605 )
# korzystamy z wybranych drzew konkurencyjnych do SDIsklad
CI28 <- G(dcent[i,],powierzchnia=pow, typ='W')/G(SDIsklad, powierzchnia=pow, typ='M')
if (nrow(SDIsklad)>0) {
CI29 <- G(dcent[i,],powierzchnia=pow, typ='W')/G(SDIsklad, powierzchnia=pow, typ='D')
} else {
CI29 <- 0
}
AP<-(G(dcent[i,], typ='W', powierzchnia=pow)/
(G(SDIsklad, powierzchnia=pow)*nrow(SDIsklad)))
} else {
SDI <- 0
CI28 <- 0
CI29 <- 0
AP <- 0
}
# wskazniki zalezne od odleglosci (krony drzew)
PK16mod <- Oimodel
PK16rzecz <- Oirzecz
}
l <- c(dcent[[nrDrzewa]][i], dcent[[nrPow]][i], CI,CI2 ,CI3,CI4 ,CI5 ,CI6,
CI7,CI8 ,CI9,CI10,CI11,CI12,CI13,CI14,CI15,CI16,CI17,CI18,
CI19,CI20,CI21,CI22,CI23,CI24,CI25,CI26,AP, CI28, CI29,
BAL,BALMOD, CI32, CI33, SDI, CI35, CI36, CI37,
ME3mod, BR1mod, BR2mod, BR3mod, BR4mod, BR5mod, PK16mod, CCmod,
ME3rzecz, BR1rzecz, BR2rzecz, BR3rzecz, BR4rzecz, BR5rzecz,
PK16rzecz, CCrzecz, ilekonkurentow )
# testowanie
#sp <- 0
#if (sp == 0) sp <- ll
#for (hh in 1:54) {
#print(paste(sp[hh], l[hh]))
#}
#print(length(l))
ll <- rbind(ll, l)
}
#if (length(brakh) > 0 ) {
#for (ii in 1:length(brakh)) {
#print(paste(c('Pow nr.', testSite[[nrPow]][1], 'drz nr.', brakh[ii], 'brak wysokosci - pomijam'), collapse = ' '))
#}
#}
# dodajemy nazwy kolumn i wierszay do dataframe'u
if (nrow(ll) > 0 ) {
names(ll) <- ll_names
nazTemp<-c(ll[,1], 'a')
rownames(ll) <- paste(dcent[[nrPow]][1], nazTemp[-length(nazTemp)], sep='_')
#usuwamy numer drzewa i nr powierzchni
ll[1:2] <-list( NULL)
}
return(ll)
}
}
dolaczAB <- function(A,B) {
# funkja laczy 2 data.frame'y w jeden sprawdzajac poprawnosc danych w
# kazdej kolumnie
if (is.data.frame(A) & nrow(A) >0 & is.data.frame(B) & nrow(B) >0) {
a.names <- names(A)
b.names <- names(B)
all.names <- union(a.names,b.names)
#print(paste("Number of columns:",length(all.names)))
a.type <- NULL
for (i in 1:ncol(A)) {
a.type[i] <- typeof(A[,i])
}
b.type <- NULL
for (i in 1:ncol(B)) {
b.type[i] <- typeof(B[,i])
}
a_b.names <- names(A)[!names(A)%in%names(B)]
b_a.names <- names(B)[!names(B)%in%names(A)]
if (length(a_b.names)>0 | length(b_a.names)>0){
print("Columns in data frame A but not in data frame B:")
print(length(a_b.names))
print("Columns in data frame B but not in data frame A:")
print(length(b_a.names))
} else if(a.names==b.names & a.type==b.type){
C <- rbind(A,B)
return(C)
}
C <- list()
for(i in 1:length(all.names)) {
l.a <- all.names[i]%in%a.names
pos.a <- match(all.names[i],a.names)
typ.a <- a.type[pos.a]
l.b <- all.names[i]%in%b.names
pos.b <- match(all.names[i],b.names)
typ.b <- b.type[pos.b]
if(l.a & l.b) {
if(typ.a==typ.b) {
vec <- c(A[,pos.a],B[,pos.b])
} else {
warning(c("Type mismatch in variable named: ",all.names[i],"\n"))
vec <- try(c(A[,pos.a],B[,pos.b]))
}
} else if (l.a) {
vec <- c(A[,pos.a],rep(NA,nrow(B)))
} else {
vec <- c(rep(NA,nrow(A)),B[,pos.b])
}
C[[i]] <- vec
}
names(C) <- all.names
C <- as.data.frame(C)
return(C)
} else if (is.data.frame(A) & nrow(A) >=0 & is.data.frame(B) & nrow(B) >=0) {
C <- rbind(A, B)
return(C)
} else {
print('Bledne typy danych do laczenia, oba zestawy danych musza byc w data.frame')
}
}
PawelKra/wskCI documentation built on May 8, 2019, 1:28 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#!/usr/bin/env Rscript
#autor: Pawel Krapiec
# biblioteka do obliczen odleglosci pkt od lini
require(spatstat)
# ustawienia srodowiska R
options(supen = 100)
# niezbedne zmienne do obliczen
typPow <- 'K' # K - kolowki; I - powierzchnie inne
kasNiep <- 'T' # czy kasowac pola nie bedace wykorzystywane do obliczen
# Domyslna wspolrzedna srodka powierchni kolowej
X <- 0 # wspolrzedna X powierzchni
Y <- 0 # wspolrzedna Y powierzchni
#--------------------------------
#DEFINICJE NAZW POL W PLIKU WEJSCIOWYM
nrPow <- 'nr.pow.' # nazwa kolumny z numerem powierzchni
d13 <- 'd' # nazwa kolumny z piersnica
h <- 'h' # nazwa kolumny z wysokoscia drzewa
hk <- 'hk' # nazwa kolumny z wysokoscia nasady korony
azymut <- 'azymut' # nazwa kolumny z kierunkiem na ktorym rosnie drzewo
odlDrzew <- 'L' # nazwa kolumny z odlegloscia drzewa od pkt
spadek <- 'nachylenie' # spadek powierzchni probnej
wystawa <- 'wystawa' # wystawa powierzchni
gat <- 'gat.' # gatunek drzewa (So , Św, Db, Jd, Bk, Db, Gb, itd.)
nrDrzewa <- 'l.p.' # numer porzadkowy drzewa na powierzchni
R <- 'R' # promien powierzchni kolowej
rKorony <- 'szer.kor.' # kolumna z rzeczywistym promieniem korony
rKoronyModel <- 'rKoronyModel' # kolumna z obliczonym promieniem korony
Xorg <- 'X' # wspolrzedna X drzewa
Yorg <- 'Y' # wspolrzedna Y drzewa
Zorg <- 'Z' # wspolrzedna Y drzewa
Xsr <- 'Xsr' # wspolrzedna X srodka kolowki
Ysr <- 'Ysr' # wspolrzedna Y srodka kolowki
# Zmienne programowe - nie zmieniac
Xo <- 'Xo' # obliczona wspolrzedna X
Yo <- 'Yo' # obliczona wspolrzedna Y
wspBL2XY <- function(x, y, kat, dist) {
# obliczenia dla powierzchni kolowych, przelicza L i azymut na X, Y
if (kat<=90) katN <- 90 - kat
else if (kat>90 & kat<=270) katN <- 360 - kat + 90
else if (kat>270) katN <- 180 - (kat - 270)
katrad <- katN * 0.017453292519
x_przelicz <- round(x + dist * cos(katrad), 2)
y_przelicz <- round( y + dist * sin(katrad), 2)
wyn <- c(x_przelicz , y_przelicz)
return(wyn)
}
zmienne <- function() {
# funkcja zwraca nazwy i wartosci zadeklarowanych zmiennych urzywanych do
# odczytu danych z plikow i obliczen
zm_gr <- c(typPow ,kasNiep ,X ,Y, nrPow ,d13 ,h ,hk ,azymut ,
odlDrzew ,spadek ,wystawa ,gat ,nrDrzewa ,R ,rKorony ,rKoronyModel,
Xorg ,Yorg ,Zorg, Xsr ,Ysr)
zmienne_ops <- c('typPow','kasNiep','X','Y','nrPow','d13','h',
'hk','azymut','odlDrzew','spadek','wystawa','gat',
'nrDrzewa','R','rKorony','rKoronyModel','Xorg','Yorg', 'Zorg',
'Xsr','Ysr')
print(paste(zmienne_ops, zm_gr, sep = " = "))
}
d2r <- function(wart) {
# przelicza podana wartosc na radiany
if (!is.numeric(wart)) stop('Podana wartosc nie jest liczba, nie mozna przeliczyc na radiany')
w <- wart * pi/180
return(w)
}
odlAB <- function(xa, ya, xb, yb) {
# obliczenie odleglosci miedzy 2 punktami
odl <- sqrt((xa - xb)**2 + (ya - yb)**2)
return(odl)
}
pktObwodu <- function( x1 = 0, y1 = 0, r=1, npoints = 100){
# oblicza punkty z okregu o zadanym promieniu
tt <- seq(0,2*pi,length.out = npoints)
xx <- x1 + r * cos(tt)
yy <- y1 + r * sin(tt)
return(data.frame(x = xx, y = yy))
}
rys_pol <- function(testSite, kor = T, pktGran='' , korR=rKoronyModel) {
# plotowanie rozmieszczenia drzew
# atrybuty:
#kor - rysowanie zasiegu koron drzew [TRUE, FALSE]
#korR - pole z promieniami koron drzew, domyslnie pole z obliczonymi
# modelowymi wartosciami [rKoronyModel]
#pktGran - dataframe z wczytanymi pkt graniczymi dla powierzchni innych
# niz kolowki
# ustawiamy podstawowy zasięg wykresu
zakresX = range(testSite[[Xo]])
zakresY = range(testSite[[Yo]])
if (typPow == 'K') {
maxR <- testSite[[R]][1]
} else {
maxR <- max(c(zakresX, zakresY))
}
# ustalamy ramke wykresu do umieszczenia w niej danych
plot(1, xlab = 'X', ylab = 'Y',
xlim = c(zakresX[1], zakresX[2]), ylim= c(zakresY[1], zakresY[2]), type = 'n', asp=1)
if (typPow == 'K') {
okrag = list()
for (i in 1:((maxR %/%5) + 1)) {
okrag = pktObwodu(0,0,5*i)
par(new=T)
lines(x=as.vector(okrag[['x']]),
y=as.vector(okrag[['y']]),
lty = 2, asp=1,
col=312 ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]))
}
par(new=T)
# dodanie granicy powierzchni probnej (kolowki)
okrag = pktObwodu(0, 0, testSite[[R]][1])
lines(x=as.vector(okrag[['x']]),
y=as.vector(okrag[['y']]), asp=1,
col='black' ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]))
}
if (kor == T && korR %in% names(testSite)){
# rysowanie koron drzew wg podanej kolumny
# jezeli wartosc jest liczba to rysujemy korone
for (i in 1:nrow(testSite)) {
if (is.numeric(testSite[[korR]][i])) {
korona <- list()
korona <- pktObwodu(testSite[[Xo]][i], testSite[[Yo]][i], testSite[[korR]][i])
par(new=T)
# dodanie granicy powierzchni probnej (kolowki)
lines(x=as.vector(korona[['x']]),
y=as.vector(korona[['y']]),
col='green' ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]), asp=1)
}
}
}
par(new=T)
symbols(x=as.vector(testSite[[Xo]]), y=as.vector(testSite[[Yo]]),
xlim = c(zakresX[1], zakresX[2]), ylim= c(zakresY[1], zakresY[2]),
circles = as.vector(testSite[[d13]]/100),
inches = F, ann=F, asp=1 ,
bg="steelblue2",
fg='black')
najblizsze = drzCentral(testSite)
par(new=T)
symbols(x=(najblizsze[[Xo]]), y=(najblizsze[[Yo]]),
xlim = c(zakresX[1], zakresX[2]), ylim= c(zakresY[1], zakresY[2]),
circles = (najblizsze[[d13]]/100),
inches = F, ann=F, asp=1 ,
bg="red",
fg='black')
# dodajemy granice powierzchni na podstawie dodanych pktGran
if (typPow=='I') {
if (is.data.frame(pktGran) && nrow(pktGran[pktGran[[nrPow]]==testSite[[nrPow]][1],])>0){
pkt_zak <- pktGran[pktGran[[nrPow]]==testSite[[nrPow]][1], c(Xorg, Yorg)]
} else {
print('Nie podano pkt graniczych, granica oparta zostanie na zewnetrznych drzewach')
pkt_zak <- testSite[, c(Xorg, Yorg)]
}
wyb_pkt_gran <- chull(pkt_zak)
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
wyb_pkt_gran <- pkt_zak[wyb_pkt_gran, c(Xorg, Yorg)]
par(new=T)
# dodanie granicy powierzchni probnej
lines(x=as.vector(wyb_pkt_gran[[Xorg]]),
y=as.vector(wyb_pkt_gran[[Yorg]]),
col='black' ,xlim = c(zakresX[1], zakresX[2]),
ylim= c(zakresY[1], zakresY[2]), asp=1)
}
text(x=as.vector(testSite[[Xo]]), y=(testSite[[Yo]]),
labels=testSite[[nrDrzewa]], asp=1,
xlim = c(zakresX[1], zakresX[2]), ylim = c(zakresY[1], zakresY[2]), pos=3, cex= 0.7)
}
rys_dh <- function(testSite) {
testSite = na.omit(testSite)
zakd13 = range(testSite[[d13]])
zakh = range(testSite[[h]])
plot(testSite[[d13]], testSite[[h]],
xlim = c(zakd13[1]-1, zakd13[2]+1), ylim= c(zakh[1]-1, zakh[2]+1),
col = 'blue', pch=19, cex=1, lty='solid', xlab="d13", ylab="h")
najblizsze = drzCentral(testSite)
par(new=T)
plot(najblizsze[[d13]], najblizsze[[h]],
xlim = c(zakd13[1]-1, zakd13[2]+1), ylim= c(zakh[1]-1, zakh[2]+1),
col = 'red', pch=19, cex=1, lty='solid', xlab='', ylab = '')
text(testSite[[d13]], testSite[[h]], labels = testSite[[nrDrzewa]],
xlim = c(zakd13[1]-1, zakd13[2]+1), ylim= c(zakh[1]-1, zakh[2]+1),
cex= 0.7, pos=3)
}
rys <- function(testSite, pktGran = '') {
par(mfcol=c(1,2))
rys_dh(testSite)
rys_pol(testSite, pktGran=pktGran)
}
G <- function(testSite, typ='A', powierzchnia, typ_pow ='m') {
# oblicza piersnicowe pole przekroju d-stanu dla roznych opcji
# typ liczenia:
# A - liczy dla wszystkich drzew, przelicza na g/ha
# M - liczy dla najgrubszego drzewa na powierzchni i przelicza na g/ha
# D - liczy dla 100 najgrubszych drzew na ha (przelicza na g/ha)
# W - liczy dla podanego drzewa i przelicz jego gi na ha
# powierzchnia - wartosc domyslnie podawana w metrach, inne wartosci [ha, km2]
# typ_pow - jednostka powierzchni w jakiej zostala przekazana do funkcji
# trigger dla powierzchni w ktorych nie sa wpisane wszystkie piersnice
dalej = 1
if (length(testSite[[d13]]) != length(na.omit(testSite[[d13]]))) {
print(testSite)
print(paste('brak wpisanych piersnic w powierzchni nr. ', testSite[[nrPow]][1], sep= ""))
print('Dla tej powierzchnie nie beda liczone statystyki')
dalej = 0
}
# przeliczenie powierzchni na odpowiednie wartosci
if (typ_pow == 'm') {
pow1<- powierzchnia / 10000 # przeliczamy powierzchnie z metrow na ha
} else if (typ_pow == 'km2'){
pow1<- powierzchnia * 100
} else if (typ_pow == 'ha'){
pow1<- powierzchnia # zostawiamy w ha
} else {
print('nierozpoznana jednostka powierzchni, dostepne sa: [m ,ha, km2]')
dalej <- 0
}
if (typ == 'D' | typ == 'A' & dalej == 1) {
#wybieramy zakres piersnic i liczymy dla nich G
if (typ == "A") {
zakd13 <- range(testSite[[d13]])
pocz <- zakd13[1] %/% 2 - 1
kon <- (zakd13[2] %/% 2) +1
ts <- testSite
} else if (typ == "D") {
ts <- testSite[order(-testSite[[d13]]),]
iledrz <- ile100( nrow(testSite), pow1)
if (nrow(ts) < iledrz) {
iledrz <- nrow(ts)
}
ts <- ts[1:iledrz, ]
zakd13 <- range(ts[[d13]])
pocz <- zakd13[1] %/% 2 - 1
kon <- (zakd13[2] %/% 2) + 1
}
g <- 0
d_n <- 0
for (i in pocz:kon) {
dd <- subset(ts,
ts[[d13]] >(i*2 ) & ts[[d13]]<= ((1+i)*2) )
# jezeli w tym przedziale nie ma zadnych piersnic idziemy dalej
if (length(dd[[d13]]>0)) {
gi <-(( (i+1)**2)*3.1415)/10000
d_n <- d_n + length(dd[[d13]])
g <- g + (gi * length(dd[[d13]]))
}
}
if (is.numeric(pow1)& pow1> 0) {
return((g/d_n)/pow1) # nie przeliczone na powierzchnie ani na nic
}
} else if (typ == 'M' & dalej == 1 ) { #"M"
ts <- testSite[order(-testSite[[d13]]),]
g <- (3.1415 * (ts[[d13]][1]/2)**2)/10000
return(g)
} else if (typ == 'W' & dalej == 1 & nrow(testSite) == 1) { # 'W'
g <- (3.1415 * (testSite[[d13]][1]/2)**2)/10000
return(g)
}
}
drzCentral <- function (testSite) {
#wyznacza drzewa centralne na powierzchni - funkcja testowa
#do modeli odleglosci wyznaczane sa w inny sposob
if (typPow == 'K') {
zakodl = range(testSite[[odlDrzew]]) #wartosc min i max odl drzew od srodka
rKol = testSite[[R]][1] # spisujemy promien kolowki do zmiennej
# sprawdzenie czy drzewo centralne spelnia minimalne warunki
if (rKol - zakodl[1] < 5) {
print('Odleglosc drzewa centralnego od granicy powierzchni jest mniejsza niz 5 m')
}
# ustalamy minimalna odleglos drzewa centralnego od granicy kolowki
if (rKol > 9.999) {
odlGran = 8
} else {
odlGran = rKol - 2
}
#wybieramy drzewa spelniajace warunek odleglosci od granicy
dcent = subset(testSite, (rKol - testSite[[odlDrzew]]) > odlGran )
} else {
dcent <- subset(testSite, odlDoGran > 8)
}
# zwracamy dataframe z drzewami spełniającymi kryteria odległosciowe
return(dcent)
}
ile100 <- function(n, pow) {
#zwraca ilosc drzew ktore nalezy wziasc z powierzchni aby wyliczyc ilosc
#drzew na ha
drz100 = 100*pow
if (drz100 < 1) {
drz100 <- 1
}
#zwracamy zaookraglona wartosc int
return(round(drz100, 0))
}
zapiszWynik <- function( wypis, nazwaPliku ) {
#zapisuje wynik w postaci pliku csv
write.table( wypis, file=nazwaPliku, col.names = T, row.names=F, sep=";")
}
wczytajPktGran <- function(f) {
# wczytuje pkt na podstawie ktorych wyznacza jest nastepnie granica pow.
data_pktgr <- read.table(f, header = T, sep=';')
return(data_pktgr)
}
wczytajPow <- function(fileIn, pktGran = data.frame()){
# kolumny niezbedne do obliczen
neededColumns <- c(nrPow, d13, h, hk, azymut, odlDrzew, gat, nrDrzewa, R,
spadek, wystawa, rKorony, Xorg, Yorg, Zorg, Xsr, Ysr)
# wczytaj dane z pliku
data<-read.table(fileIn, header = T, sep=";")
nAll <- nrow(data[1]) # liczba rekordow w pliku
nSites <- length(unique(data[nrPow])) # ile powierzchni w pliku
listNames <- unique(data[nrPow]) # lista nazw powierzchni
# usuwanie zbednych kolumn z danych wejsciowych
if (kasNiep == "T") {
data <- data[,(names(data) %in% neededColumns)]
}
Xo <- c()
Yo <- c()
rKoronyModel <- c()
# wskazniki do obliczania szerokosci korony wg modelu
a <- 0
b <- 0
for (i in 1:nAll) {
if (typPow == 'K') {
# obliczamyt wspolrzedne przy kolowkach
if (Xsr %in% names(data) & Ysr %in% names(data)){
X <- data[[Xsr]][i]
Y <- data[[Ysr]][i]
}
przel <- wspBL2XY(X, Y, data[[azymut]][i], data[[odlDrzew]][i])
Xo <- c(Xo, przel[1])
Yo <- c(Yo, przel[2])
} else {
# przepisujemy wspolrzedne orginalne jako obliczone dla zgodnosci
Xo <- c(Xo, data[[Xorg]][i])
Yo <- c(Yo, data[[Yorg]][i])
}
if (grepl(data[[d13]][i], ",")){
print("Dane liczbowe oddzielone są przecinkami, powiny być kropki")
break
}
#TODO: dodanie wspolczynnikow dla innych gatunkow
if (data[i,gat] == 'So') {
a<--0.1797
b<-0.6267
} else if(data[i,gat] == 'Jd') {
a<-0.092
b<-0.538
} else if(data[[gat]][i] == 'Św') {
a<--0.3232
b<-0.6441
} else if(data[i,gat] == 'Db') {
a<--0.3973
b<-0.7328
} else if(data[i,gat] == 'Bk') {
a<-0.1366
b<-0.6183
} else if(data[i,gat] == 'Gb') { # wartosci skopiowane z BK
a<-0.1366
b<-0.6183
} else if(data[i,gat] == 'Js') {
a<--0.1366
b<-0.7328
} else { # jezeli zadne z powyzszych to brak wspolczynnikow
a <- NA
b <- NA
}
if (!is.na(a) & data[[d13]][i] > 0 ){
rKoronyModel <- c(rKoronyModel , round((0.5 * exp(a + b* log(data[[d13]][i]))),2))
} else { # jezeli gatunek jest inny niz powyzsze - brak wspolczynnikow do obliczen
rKoronyModel <- c(rKoronyModel, 0 )
}
}
#dodanie obliczonych koordynatowa do orginalnej tabeli
data <- cbind(data, Xo, Yo)
# dodanie modelowych wielkosci korony
data <- cbind(data, rKoronyModel)
Xo = 'Xo' # nazwa kolumny z obliczonymi wsp X dla kazdej powierzchni
Yo = 'Yo' # nazwa kolumny z obliczonymi wsp Y dla kazdej powierzchni
rKoronyModel <- 'rKoronyModel' # nazwa kolumny z wyliczonym promieniem korony
b = list()
for (i in (1:length(listNames[,1]))) {
if (sapply(subset(data, data[nrPow]== listNames[i,1]), is.numeric)[[d13]] == T){
t1 = subset(data, data[nrPow] == listNames[i,1])
# ustawiamy nazwy wierszy na bardziej sensowne
rownames(t1) <- paste(t1[[nrPow]], t1[[nrDrzewa]], sep='_')
# oblicz odleglosc drzew do granicy powierzchni
if(typPow == 'K'){
# czy kolowka ma podane orginalne wspolrzedne srodka
if (Xorg %in% names(t1) & Yorg %in% names(t1)){
Xt <- t1[[Xsr]][1]
Yt <- t1[[Ysr]][1]
} else { # jesli nie ustawiamy na 0
Xt <- 0
Yt <- 0
}
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], pktObwodu(Xt, Yt, t1[[R]][1]))
} else {
if (is.data.frame(pktGran)== T){
if (nrow(pktGran[pktGran[[nrPow]]==t1[[nrPow]][1], ]) > 0){
pkt_zak <- pktGran[pktGran[[nrPow]]==t1[[nrPow]][1], c(Xorg, Yorg)]
wyb_pkt_gran <- chull(pktGran[pktGran[[nrPow]]==t1[[nrPow]][1], c(Xorg, Yorg)])
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
#pkt_zak[wyb_pkt_gran, c('X', 'Y')]
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], pkt_zak[wyb_pkt_gran, c(Xorg, Yorg)])
} else {
print('Brak danych dla pkt pomiarowych granicy powierzchni')
print('Obliczam odleglosci do granicy wyznaczonej na podstawie drzew ')
wyb_pkt_gran <- chull(t1[c(Xo, Yo)])
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
#pkt_zak[wyb_pkt_gran, c(Xo, Yo)]
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], t1[wyb_pkt_gran, c(Xo, Yo)])
}
} else {
print('Brak danych dla pkt pomiarowych granicy powierzchni')
print('Obliczam odleglosci do granicy wyznaczonej na podstawie drzew ')
wyb_pkt_gran <- chull(t1[c(Xo, Yo)])
wyb_pkt_gran <- c(wyb_pkt_gran, wyb_pkt_gran[1])
#pkt_zak[wyb_pkt_gran, c(Xo, Yo)]
odlDoGran <- oblOdlGran(t1[[Xo]], t1[[Yo]], t1[wyb_pkt_gran, c(Xo, Yo)])
}
}
#obliczamy korekte wysokosci drzewa w zaleznosci od spadku na kolowce
if(typPow == 'K' & spadek %in% names(t1) & wystawa %in% names(t1)){
poprH <- oblPoprWys(t1, t1[[wystawa]][1], t1[[spadek]][1])
} else if (typPow == 'I' & Zorg %in% names(t1)){
roz <- range(t1[[Zorg]])
srodekH <- (roz[1] + roz[2]) / 2
# obliczamy przesuniecia wysokosciowe drzew
#transform(t1[[Zorg]], odlDoGran = Zorg - srodekH)
poprH <- sapply(t1[[Zorg]], function(x){return(round(x-srodekH, 2))})
} else {
#jezeli nie ma odpowiednich danych wstawiamy 0
poprH <- c(rep.int(0, nrow(t1)))
}
b[[toString(listNames[i,1])]] <- cbind(t1, odlDoGran, poprH)
} else {
print(paste('Brak wszystkich piersnic - pow. nr.: ', toString(listNames[i,1]) ))
}
}
odlDoGran <- 'odlDoGran'
return(b)
}
obliczCzWsp <- function(r, R, d) {
# oblicza czesc wspolna dwoch przecinajacych sie okregow
if (r>R) {
# zamieniamy promienie miejsacami
tt <- r
r <- R
R <- tt
} else if (r == R) {
R <- R + 0.0001
}
if (r+d < R) {
# jezeli drzewo jest calkowicie pod drzewem centralnym to wchodzi cala pow
wynik <- pi * r*r
} else {
# jezeli drzewa nachodza na siebie tylko w malej czesci to liczymy
cz1 = r**2*acos((d**2+r*r-R*R)/(2*d*r))
cz2 <- R*R*acos((d*d+R*R-r*r)/(2*d*R))
cz3 <- 0.5*sqrt((-d+r+R)*(d+r-R)*(d-r+R)*(d+r+R))
wynik <- cz1 + cz2 -cz3
if (is.nan(cz1) & is.nan(cz2) & is.nan(cz3)){
print(paste(r, R, d))
}
# aleternatywny sposob liczenia - zostawiony na wypadek przyszlych modyfikacji
#x1<-(d**2-r*r+R*R)/(2*d)
#print(paste('x1= ',x1))
#x2<-abs(d-x1)
#print(paste('x2= ', x2))
#y<- sqrt(R*R- (x1*x1))
#print(paste('y= ', y))
#a1<- (R**2)*acos(x1/R) - x1 * y
#print(paste('a1= ', a1))
#a2<- r**2 * acos(x2/r) - x2 * y
#print(paste('a2= ', a2))
#if (x1>d){
#a2<- pi * r**2 - a2
#print(paste('a2_podm = ',a2))
#}
#area <- a1+a2
#print(paste('pow= ', area))
}
return(wynik)
}
oblOdlGran <- function (X, Y, pktGran) {
# zwraca liste z odleglosciami kolejnych drzew od granicy powierzchni
# pktKran - dwuwymiarowa maciez z pkt opisujacymi granice powierzchni
drz <- cbind( X, Y)
odlGran <- apply(distppll(drz, cbind(pktGran[-nrow(pktGran), ], pktGran[-1,])),1,min)
return(round(odlGran, 2))
}
oblPoprWys <- function(testSite, wystawa = 0, nachyl = 0) {
# oblicza wysokosci odziomkow na powierzchni kolowej na podstawie nachylenia
# i wystawy podanych w stopniach
if (Xsr %in% names(testSite) & Ysr %in% names(testSite)){
xpocz <- testSite[[Xsr]][1]
ypocz <- testSite[[Ysr]][1]
} else {
xpocz <- 0
ypocz <- 0
}
x_zakres <- range(testSite[[Xo]])
y_zakres <- range(testSite[[Yo]])
if (wystawa <= 90 ) {
pkt1 <- c(x_zakres[1], y_zakres[1])
} else if (wystawa > 90 & wystawa <= 180) {
pkt1 <- c(x_zakres[1], y_zakres[2])
} else if (wystawa >180 & wystawa <= 270) {
pkt1 <- c(x_zakres[2], y_zakres[2])
} else if (wystawa > 270 ){
pkt1 <- c( x_zakres[2], y_zakres[1])
}
pkt2 <- wspBL2XY(pkt1[1], pkt1[2], wystawa, 1500)
# obliczamy dlugosc wektora
odlPP <- odlAB(pkt1[1],pkt1[2], pkt2[1], pkt2[2])
# obliczamy spadek dla srodka kolowki
u0 <- (( xpocz - pkt1[1]) * (pkt2[1] - pkt1[1]) + (ypocz - pkt1[2]) *(pkt2[2] - pkt1[2])) /
((pkt1[1] - pkt2[1])**2 + (pkt1[2] - pkt2[2])**2 )
spadek0 <- round(tan(d2r(nachyl)) * u0 * odlPP, 2)
#obliczamy spadki
ulist <- c()
for (i in 1:nrow(testSite)) {
u <- ((testSite[[Xo]][i] - pkt1[1]) * (pkt2[1] - pkt1[1]) + (testSite[[Yo]][i] - pkt1[2]) *(pkt2[2] - pkt1[2])) /
((pkt1[1] - pkt2[1])**2 + (pkt1[2] - pkt2[2])**2 )
#obliczamy poprawke wysokosci drzewa na spadek i ustawiamy 0 na srodku kolowki
spadek <- round((tan(d2r(nachyl)) * u * odlPP) - spadek0, 2)
ulist <- rbind(ulist, spadek)
}
return(ulist)
}
modeleCI <- function(testSite, debug = 1) {
#dla podanej powierzchni zwraca obliczenia dla drzew centralnych wyzanczonych
#na podstawie wybranych modeli odleglosci.
#Dostepne modele:
#odleglosciowe: 3m, 4m, 5m, 6m, 7m, 8m, 9m, 10,
#odleglosci zaleznej od piersnic: S8, S6
#odl. zaleznej od korony: 25s, 35pr
#odl. zalezne od wysokosci drzewa: 02H, 03H, 04H, 05H
#odl. zalezne od d13: 015D, 020D, 025D, 030D
#sr. odl. Oa
#4 najblizsze
#odwroconego stozka, zalezne od a i kata rozwarcia stozka
#modele stalego r od 3 do 10 metrow
testSite$lacz<-paste(testSite[[nrPow]], testSite[[nrDrzewa]], sep='_')
wyn <- testSite
# deklaracja nazw dla poszczegolnych tablic konkurencji i modeli
konk_list_names <- c('s8_konk', 's6_konk', 's25_konk', 'pr35_konk', 'sm25_konk', 'prm35_konk',
'H02_konk', 'H03_konk', 'H04_konk', 'H05_konk', 'D015_konk',
'D020_konk', 'D025_konk', 'D030_konk', 'n2_konk', 'n4_konk', 's_konk', 'sm_konk')
st_list <- c('0150', '0160', '0170', '0450', '0460', '0470', '0650',
'0660', '0670', '06650', '06660', '06670', 'nk50', 'nk60', 'nk70')
model_list_names <- list('s8','s6','s25','pr35','sm25','prm35',
'H02','H03','H04','H05','D015','D020','D025',
'D030','n2','n4', 's', 'sm')
for (j in 3:10) {
dcent <- subset(testSite, testSite[['odlDoGran']] > j)
odl <- j
if (nrow(dcent)>0){
mod <- wskCI(dcent, testSite, odl)
names(mod) <- paste(paste('m', toString(j),sep=""), names(mod), sep='_')
mod$lacz <- rownames(mod)
print(paste('Obliczam model o stalym promieniu dla:', j, 'm', 'drzew cent:', nrow(mod) ))
wyn <- merge(mod, wyn, all.y=T)
}
}
# modele oparte o piersnice
dcent <- subset(testSite, testSite[['odlDoGran']] > 7 & h>0)
if (nrow(dcent) == 0){
print(paste('Pow.',testSite[[nrPow]][1] , '- Brak drzew spelniajacych kryteria drzewa centralnego', sep=' '))
return(data.frame())
} else {
print(paste('Do analizy wybrano:', nrow(dcent), '/', nrow(testSite), 'drzew'))
#zerowanie wartosci modeli dla kazdego drzewa
mod_s8 <- data.frame()
mod_s6 <- data.frame()
mod_s25 <- data.frame()
mod_pr35 <- data.frame()
mod_sm25 <- data.frame()
mod_prm35 <- data.frame()
mod_H02 <- data.frame()
mod_H03 <- data.frame()
mod_H04 <- data.frame()
mod_H05 <- data.frame()
mod_D015 <- data.frame()
mod_D020 <- data.frame()
mod_D025 <- data.frame()
mod_D030 <- data.frame()
mod_n2 <- data.frame()
mod_n4 <- data.frame()
mod_s <- data.frame()
mod_sm <- data.frame()
mod_st <- list()
# grupujemy modele w liste
model_list <- list(mod_s8,mod_s6,mod_s25,mod_pr35,mod_sm25,mod_prm35,
mod_H02,mod_H03,mod_H04,mod_H05,mod_D015,mod_D020,mod_D025,
mod_D030,mod_n2,mod_n4, mod_s, mod_sm)
#print(dcent)
for (j in 1:nrow(dcent)) {
# obliczenia dla kazdego drzewa centralnego po kolei
print(' ')
print(paste('Obliczam drzewo centralne nr:', j, '/', nrow(dcent), '(', dcent[j,nrDrzewa], ')'))
xmin <- dcent[[Xo]][j]-40
xmax <- dcent[[Xo]][j]+40
ymin <- dcent[[Yo]][j]-40
ymax <- dcent[[Yo]][j]+40
# pobieramy wycinek powierzchnie o boku 40m do obliczen ze stozkiem
wycinek20 <- subset(testSite, Xo< xmax & Xo> xmin & Yo< ymax & Yo> ymin)
# obliczamy odleglosci do tego drzewa centralnego
odltemp <- c()
for (ii in 1:nrow(wycinek20)) {
odlobl <- odlAB(dcent[[Xo]][j],
dcent[[Yo]][j],
wycinek20[[Xo]][ii],
wycinek20[[Yo]][ii])
odltemp <- c(odltemp,round(odlobl, 2) )
}
# dopisujemy odleglosci do wycinka i sortujemy od najmniejszej odleglosci
wycinek20t <- cbind(wycinek20, odltemp)
odltemp<-'odltemp'
wycinek20 <- wycinek20t[order(wycinek20t[[odltemp]]),]
Oa <- min(wycinek20[[odltemp]][2])*2
#print(wycinek20)
# zerowanie grup konkurencyjnych dla poszczegolnych drzew
s6_konk <- data.frame()
s8_konk <- data.frame()
s25_konk <- data.frame()
sm25_konk <- data.frame()
pr35_konk <- data.frame()
prm35_konk <- data.frame()
H02_konk<- data.frame()
H03_konk<- data.frame()
H04_konk<- data.frame()
H05_konk<- data.frame()
D015_konk<- data.frame()
D020_konk<- data.frame()
D025_konk<- data.frame()
D030_konk<- data.frame()
n2_konk <- data.frame()
n4_konk <- data.frame()
s_konk <- data.frame()
sm_konk <- data.frame()
for (jj in 1:nrow(wycinek20)) {
# obliczamy odleglosc do drzewa konkurencyjnego
odlRobocza <- odlAB(dcent[[Xo]][j],
dcent[[Yo]][j],
wycinek20[[Xo]][jj],
wycinek20[[Yo]][jj])
# obliczamy grupe konkurencyjna dla modelu s6
if ((dcent[[d13]][j] + wycinek20[[d13]][jj]) / 6 > odlRobocza) {
s6_konk <- rbind(s6_konk, wycinek20[jj, ])
}
# obliczamy grupe konkurencyjna dla modelu s8
if ((dcent[[d13]][j] + wycinek20[[d13]][jj]) / 8 > odlRobocza) {
s8_konk <- rbind(s8_konk, wycinek20[jj, ])
}
# czy drzewo posiada pomiary szerokosci korony
if (rKorony %in% names(wycinek20) & is.numeric(wycinek20[[rKorony]][jj])) {
if ( wycinek20[[rKorony]][jj] > 0){
if (odlRobocza < 2.5 * 2 * dcent[[rKorony]][j]) {
s25_konk <- rbind(s25_konk, wycinek20[jj, ])
}
#obliczamy grupe konkurencyjna dla modelu pr35
if (odlRobocza < 3.5 * dcent[[rKorony]][j]) {
pr35_konk <- rbind(pr35_konk, wycinek20[jj, ])
}
# obliczamy grupe konkurencyjna dla modelu s
if (odlRobocza < dcent[[rKorony]][j] + wycinek20[[rKorony]][jj] &
!is.na(dcent[[hk]][j]) & !is.na(wycinek20[[h]][jj]) & !is.na(wycinek20[[hk]][jj]) &
is.numeric(wycinek20[[h]][jj]) & is.numeric(wycinek20[[hk]][jj]) &
is.numeric(dcent[[h]][j])){
if ((dcent[[hk]][j]) <=
(wycinek20[[h]][jj] - dcent[['poprH']][j] + wycinek20[['poprH']][jj])){
s_konk <- rbind(s_konk, wycinek20[jj,])
}
}
}
}
# czy drzewo posiada modelowe szerokosci korony
if (rKoronyModel %in% names(wycinek20) & wycinek20[[rKoronyModel]][jj] > 0 &
is.numeric(wycinek20[[rKoronyModel]][jj])) {
if (odlRobocza < 2.5 * 2 * dcent[[rKoronyModel]][j]) {
sm25_konk <- rbind(sm25_konk, wycinek20[jj, ])
}
#obliczamy grupe konkurencyjna dla modelu m25s
if (odlRobocza < 3.5 * dcent[[rKoronyModel]][j]) {
prm35_konk <- rbind(prm35_konk, wycinek20[jj, ])
}
# obliczamy grupe konkurencyjna dla modelu sm
if (odlRobocza < dcent[[rKoronyModel]][j] + wycinek20[[rKoronyModel]][jj] &
!is.na(dcent[[hk]][j]) & !is.na(wycinek20[[h]][jj]) & !is.na(wycinek20[[hk]][jj]) &
is.numeric(wycinek20[[h]][jj]) & is.numeric(wycinek20[[hk]][jj]) &
is.numeric(dcent[[h]][j])){
if ((dcent[[hk]][j]) <=
(wycinek20[[h]][jj] - dcent[['poprH']][j] + wycinek20[['poprH']][jj])){
sm_konk <- rbind(sm_konk, wycinek20[jj,])
}
}
}
# obliczamy grupe konkurencyjna w zaleznosci od wysokosci drzewa centralnego
if (!is.na(dcent[[h]][j]) & is.numeric(dcent[[h]][j]) &
!is.na(wycinek20[[h]][j]) & is.numeric(wycinek20[[h]][j]) & wycinek20[[h]][j]>0) {
if (odlRobocza < dcent[[h]][j] * 0.2) {
H02_konk <- rbind(H02_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.3) {
H03_konk <- rbind(H03_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.4) {
H04_konk <- rbind(H04_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.5) {
H05_konk <- rbind(H05_konk, wycinek20[jj,])
}
}
# grupy konkurencyjne oparte na piersnicach
if (odlRobocza < dcent[[h]][j] * 0.15) {
D015_konk <- rbind(D015_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.20) {
D020_konk <- rbind(D020_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.25) {
D025_konk <- rbind(D025_konk, wycinek20[jj,])
}
if (odlRobocza < dcent[[h]][j] * 0.30) {
D030_konk <- rbind(D030_konk, wycinek20[jj,])
}
# model sredniej odleglosci Oa
if (Oa > 0) {
if (odlRobocza < Oa) {
n2_konk <- rbind(n2_konk, wycinek20[jj,])
}
}
}
#wybieramy grupe konkurencyjna dla modelu 4n
ile_wyb <- 1
if (nrow(wycinek20) >= 5 & hk %in% names(wycinek20) & h %in% names(wycinek20)) {
for (l in 2:nrow(wycinek20)) {
if (is.numeric(wycinek20[[h]][l]) & is.numeric(wycinek20[[hk]][l])&
is.numeric(dcent[[h]][j]) & !is.na(dcent[[hk]][j])&
!is.na(wycinek20[[h]][l]) & !is.na(wycinek20[[hk]][l])) {
if ((dcent[[hk]][j]) <=
(wycinek20[[h]][l] - dcent[['poprH']][j] + wycinek20[['poprH']][l]) & ile_wyb <5){
n4_konk <- rbind(n4_konk, wycinek20[l,])
}
# jezeli drzewo nie weszlo to i tak je liczymy
ile_wyb <- ile_wyb + 1
}
}
}
# konkurencja dla modeli zwiazanych ze stozkiem
wys_stozka <- c(0.1, 0.4, 0.6, 0.66)
# obliczamy tangens kata alfa
tan50 <- tan(d2r(90-(50/2)))
tan60 <- tan(d2r(90-(60/2)))
tan70 <- tan(d2r(90-(70/2)))
lista_k <- c(tan50, tan60, tan70)
st_konk <- list()
nk_trig <- -999 # trigger dla nasady korony
if (h %in% names(dcent) & hk %in% names(dcent)){
for (l in 1:nrow(wycinek20)) {
if (is.numeric(wycinek20[[h]][l]) & is.numeric(wycinek20[[hk]][l])&
is.numeric(dcent[[h]][j]) & !is.na(dcent[[hk]][j])&
!is.na(wycinek20[[h]][l]) & !is.na(wycinek20[[hk]][l])) {
for (ll in 1:length(lista_k)){
for (kk in 1:length(wys_stozka)) {
# skladamy nazwe grupy konkurencyjnej i dodajemy pod nia dane konkurencji
naz_temp <- paste(0,
unlist(strsplit(toString(wys_stozka[kk]),
split='.', fixed=T))[2],
toString((ll+4)*10), sep="")
#sprawdzamy czy warunek jest dodatni lub zerowy
if (((wycinek20[[h]][l] + wycinek20[['poprH']][l] - dcent[['poprH']][j]) -
(wys_stozka[kk]* (dcent[[h]][j])))/
lista_k[ll] >0){
war_wys <- ( ((wycinek20[[h]][l] + wycinek20[['poprH']][l]- dcent[['poprH']][j]) -
(wys_stozka[kk]* (dcent[[h]][j])))/
lista_k[ll])
} else {
war_wys <- 0.0000000000000001
}
if (wycinek20[[odltemp]][l] < war_wys & wycinek20[[odltemp]][l] >0 ) {
st_konk[[naz_temp]] <- rbind(st_konk[[naz_temp]], wycinek20[l,])
}
# obliczenia dla nasady korony
nazkor_temp <- paste('nk', toString((ll+4)*10), sep="")
if (((wycinek20[[h]][l] - dcent[['poprH']][j]+ wycinek20[['poprH']][l] )-
(dcent[[hk]][j]))>0 ){
war_wys_nk <- (((wycinek20[[h]][l] - dcent[['poprH']][j]+ wycinek20[['poprH']][l] )-
(dcent[[hk]][j]))/
lista_k[ll])
} else {
war_wys_nk <- 0.00000000000001
}
#print(paste(nazkor_temp, wycinek20[[odltemp]][l], war_wys_nk))
if (wycinek20[[odltemp]][l] < war_wys_nk & wycinek20[[odltemp]][l] >0 &
nk_trig != ll) {
st_konk[[nazkor_temp]] <- rbind(st_konk[[nazkor_temp]],wycinek20[l,])
# trigger pomijajacy jedna petlekk
nk_trig <- ll
}
}
}
}
}
}
# obliczenia dla standardowych modeli
konk_list <- list(s8_konk, s6_konk, s25_konk, pr35_konk, sm25_konk, prm35_konk,
H02_konk, H03_konk, H04_konk, H05_konk, D015_konk,
D020_konk, D025_konk, D030_konk, n2_konk, n4_konk, s_konk, sm_konk)
names(konk_list) <- konk_list_names
for (ii in 1:length(konk_list)) {
if (is.data.frame(konk_list[[konk_list_names[ii]]]) && !nrow(konk_list[[konk_list_names[ii]]])==0) {
model_list[[ii]] <- rbind(model_list[[ii]], wskCI(dcent[j,], konk_list[[konk_list_names[ii]]], odlMod=range(konk_list[[konk_list_names[ii]]][['odltemp']])[2]))
if (debug == 1) {
print(paste(' Model ', konk_list_names[ii], ', drzew konk: ', nrow(konk_list[[konk_list_names[ii]]]), ', max odl: ', range(konk_list[[konk_list_names[ii]]][['odltemp']])[2], sep=''))
}
}
}
# obliczenia dla modeli ze stozkiem
for (item in st_list) {
if ((is.data.frame(st_konk[[item]]))&& !nrow(st_konk[[item]])==0) {
mod_st[[item]] <- rbind(rbind(mod_st[[item]], wskCI(dcent[j,], st_konk[[item]], odlMod=range(st_konk[[item]][['odltemp']])[2])))
if (debug == 1) {
print(paste(' Model ', item , ', drzew konk: ', nrow(st_konk[[item]]), ', max odl: ', range(st_konk[[item]][['odltemp']])[2], sep=''))
}
}
}
}
for (model in 1:length(model_list)) {
if (nrow(model_list[[model]]) >0) {
names(model_list[[model]]) <- paste(model_list_names[[model]], names(model_list[[model]]), sep='_')
model_list[[model]]$lacz <- rownames(model_list[[model]])
wyn <- merge(model_list[[model]], wyn, all.y=T, all.x=T)
}
}
for (model in names(mod_st)) {
if (nrow(mod_st[[model]]) > 0) {
names(mod_st[[model]]) <- paste(model, names(mod_st[[ model ]]), sep='_')
mod_st[[model]]$lacz <- rownames(mod_st[[model]])
wyn <- merge(mod_st[[model]], wyn, all.y=T, all.x=T)
}
}
# zamieniamy wszystkie NA na 0
wyn[is.na(wyn)] <- 0
# zwracamy ostateczna tabele
return(wyn)
}
}
wskCI <- function(dcent, testSite, odlMod = 6) {
# funkcja oblicza wskazniki konkurencjii w dla drzew centralnych zapisanych jako
#dcent (data.frame) spelniajacych kryteria odleglosci od granicy powierzchni.
#Wskazniki oblczane sa na podstawie drzew konkurencyjnych testSite (data.frame)
#ktore wystepuja w podanej odleglosci od drzewa centralnego - odlMod (numeric)
# wyznaczamy drzewa centralne - dla celow testowych
if (!is.data.frame(dcent)){
dcent <- drzCentral(testSite)
}
# lista drzew bez wpisanej wysokosci, na koncu wypiszemy komunikat
brakh <-c()
if (nrow(dcent) > 0) {
ll_names <- list('NrDrz', 'NrPow', 'CI','CI2' ,'CI3','CI4' ,'CI5' ,'CI6' ,'CI7','CI8' ,
'CI9','CI10','CI11','CI12','CI13','CI14','CI15','CI16','CI17',
'CI18','CI19','CI20','CI21','CI22','CI23','CI24','CI25','CI26',
'AP', 'CI28', 'CI29', 'BAL','BALMOD', 'CI32', 'CI33',
'SDI', 'CI35', 'CI36', 'CI37',
'ME3mod', 'BR1mod', 'BR2mod', 'BR3mod', 'BR4mod', 'BR5mod', 'PK16mod', 'CCmod',
'ME3rzecz', 'BR1rzecz', 'BR2rzecz', 'BR3rzecz', 'BR4rzecz', 'BR5rzecz',
'PK16rzecz', 'CCrzecz', 'Konkurenci' )
ll<- data.frame()
for (i in 1:nrow(dcent)) {
CI <- 0
CI2 <- 0
CI3 <- 0
CI4 <- 0
CI5 <- 0
CI6 <- 0
CI7 <- 0
CI8 <- 0
CI9 <- 0
CI10 <- 0
CI11 <- 0
CI12 <- 0
CI13 <- 0
CI14 <- 0
CI15 <- 0
CI16 <- 0
CI17 <- 0
CI18 <- 0
CI19 <- 0
CI20 <- 0
CI21 <- 0
CI22 <- 0
CI23 <- 0
CI24 <- 0
CI25 <- 0
CI26 <- 0
BAL <- 0
sCI35 <- 0
nCI35 <- 0
CI35 <- 0
CI36 <- 0
gCI37 <- 0
CI37 <- 0
CI32 <- 0
CI33 <- 0
Oirzecz <- 0
Oimodel <- 0
di <- 0
ME3rzecz <- 0
BRskl <- 0
BR1rzecz <- 0
BR2rzecz <- 0
BR3rzecz <- 0
BR4rzecz <- 0
BR5rzecz <- 0
CCrzecz <- 0
ME3mod <- 0
BRskl <- 0
AP <- 0
CI28 <- 0
CI29 <- 0
SDI <- 0
BALMOD <- 0
BR1mod <- 0
BR2mod <- 0
BR3mod <- 0
BR4mod <- 0
BR5mod <- 0
CCmod <- 0
PK16mod <- 0
PK16rzecz <- 0
SDIsklad <- 0
SDItrig <- 0
ilekonkurentow <- 0
if (odlMod == 999) {
pr_pow <- 40
} else {
pr_pow <- odlMod
}
# obliczamy powierzchnie wybranego wycinka pow. probnej dla drzewa centralego
pow <- 3.1415*pr_pow**2
# czy w danych wystepuja rzeczywiste pomiary korony drzewa
koronaLiczenie <- 1
if (rKorony %in% colnames(testSite)) {
koronaLiczenie <- 2
}
for (j in 1:nrow(testSite)) {
# odleglosi miedzy drzewami w celu wyznaczeniea drzew konkurencyjnych
# w zasiegu analizowanego drzewa
odl<-odlAB(dcent[[Xo]][i],
dcent[[Yo]][i],
testSite[[Xo]][j],
testSite[[Yo]][j])
# promienie konkurencji
#if (odl < 3) {r <- 3} else {r <- as.integer(odl)}
r <- pr_pow
if (odl > 0 & odl < pr_pow ) {
ilekonkurentow <- ilekonkurentow + 1
# podzbior do obliczenia SDI
if (SDItrig == 0) {
SDIsklad <- testSite[j,]
SDItrig <- 1
} else {
SDIsklad <- rbind(SDIsklad, testSite[j,])
}
#piersnice wieksze od d centralnego i modyfikacje
if (testSite[[d13]][j] > dcent[[d13]][i]) {
BAL <- BAL + G(testSite[j,], powierzchnia=pow, typ="W")
CI11 <- CI11 + atan(testSite[[d13]][j]/odl)**-1
CI25 <- CI25 + ( testSite[[d13]][j]-dcent[[d13]][i] ) / odl
CI36 <- CI36 + (testSite[[d13]][j] / dcent[[d13]][i])**0.5
} else if (dcent[[d13]][i] == testSite[[d13]][j]){
CI36 <- CI36 + (testSite[[d13]][j] / dcent[[d13]][i])**1
} else if (dcent[[d13]][i] > testSite[[d13]][j]){
CI36 <- CI36 + (testSite[[d13]][j] / dcent[[d13]][i])**2
}
#wysokosci - wieksze i wszystkie inne o ile istnieja
if (!is.na(testSite[[h]][j]) & !is.na(dcent[[h]][i])) {
if (testSite[[h]][j] > dcent[[h]][i]) {
CI33 <- CI33 + G(testSite[j,], powierzchnia=pow, typ='W')
CI12 <- CI12 + atan(testSite[[d13]][j]/odl)**-1
CI15 <- CI15 + atan(testSite[[h]][j]/odl)**-1
CI19 <- CI19 + atan(abs((testSite[[h]][j]-dcent[[h]][i])/odl))**-1
}
if (testSite[[h]][j]>dcent[[h]][i]/2) {
CI14 <- CI14 + atan(testSite[[h]][j]/odl)**-1
}
if (testSite[[h]][j]>dcent[[h]][i]*0.8) {
CI18 <- CI18 + atan(abs((testSite[[h]][j] - 0.8*dcent[[h]][i])/odl))**-1
}
CI2 <- CI2 + testSite[[h]][j]/(dcent[[h]][i]*odl)
CI13 <- CI13 + atan(testSite[[h]][j]/odl)**(-1)
CI17 <- CI17 + (testSite[[h]][j]/dcent[[h]][i])*atan(d2r(testSite[[h]][j]/odl))**(-1)
} else { # jak nie istnieja to pomijamy i wypisujemy
if (!is.na(testSite[[h]][j])) {
if (!(testSite[[nrDrzewa]][j] %in% brakh)) {
brakh <- c(brakh, testSite[[nrDrzewa]][j] )
}
} else{
if (!(dcent[[nrDrzewa]][i] %in% brakh)) {
brakh <- c(brakh, dcent[[nrDrzewa]][i] )
}
}
}
di <- di + testSite[[d13]][j]
# obliczenia zalezne tylko i wylacznie od odleglosci
CI <- CI + testSite[[d13]][j] / (dcent[[d13]][i] * odl)
CI3 <- CI3 + testSite[[d13]][j]/(dcent[[d13]][i]*(odl+1))
CI4 <- CI4 + abs((testSite[[d13]][j]/dcent[[d13]][i])*exp(-(16*odl)/(testSite[[d13]][j]+dcent[[d13]][i])))
CI5 <- CI5 + G(testSite[j,],powierzchnia=pow, typ='W')/(odl/r)
CI6 <- CI6 + (G(testSite[j,],powierzchnia=pow, typ='W')/
G(dcent[i,],powierzchnia=pow, typ='W'))/(odl/r)
CI7 <- CI7 + (testSite[[d13]][j]/dcent[[d13]][i])/(odl/r)
CI8 <- CI8 + (testSite[[d13]][j]/dcent[[d13]][i])/sqrt(odl/r)
CI9 <- CI9 + (testSite[[d13]][j]/dcent[[d13]][i])/(odl/r)**2
CI10 <- CI10 + atan(testSite[[d13]][j]/odl)**(-1)
CI16 <- CI16 + (testSite[[d13]][j]/dcent[[d13]][i])*atan(testSite[[d13]][j]/odl)**(-1)
CI20 <- CI20 + (testSite[[d13]][j] / odl)
CI21 <- CI21 + ( testSite[[d13]][j]/dcent[[d13]][i] ) / odl**2
CI22 <- CI22 + ( testSite[[d13]][j]/dcent[[d13]][i] )**2 / odl
CI23 <- CI23 + odl
CI24 <- CI24 + testSite[[d13]][j] / odl
#niezalezne od odleglosci - teoretycznie
CI32 <- CI32 + G(testSite[j,],powierzchnia=pow, typ='W')
sCI35 <- sCI35 + G(testSite[j,],powierzchnia=pow, typ='W')
nCI35 <- nCI35 + 1
gCI37 <- gCI37 + G(testSite[j,],powierzchnia=pow, typ='W')
# wskazniki zwiazane z korona drzewa
if (odl< ( testSite[[rKoronyModel]][j]+dcent[[rKoronyModel]][i] ) &
dcent[[rKoronyModel]][i] >0 & testSite[[rKoronyModel]][j] > 0 ) {
Oi <- obliczCzWsp(dcent[[rKoronyModel]][i], testSite[[rKoronyModel]][j], odl)
#print(paste( testSite[[rKoronyModel]][j]+dcent[[rKoronyModel]][i], odl ))
Oimodel <- Oimodel + (Oi/(3.1415*dcent[[rKoronyModel]][i]**2))
ME3mod <- ME3mod + (Oi/(3.1415*dcent[[rKoronyModel]][i]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BRsklm <- (Oi/(3.1415*dcent[[rKoronyModel]][i]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BR1mod <- BR1mod + BRsklm**1
BR2mod <- BR2mod + BRsklm**1.5
BR3mod <- BR3mod + BRsklm**2
BR4mod <- BR4mod + BRsklm**2.5
BR5mod <- BR5mod + BRsklm**3
CCmod <- CCmod +(3.1415*testSite[[rKoronyModel]][j]**2)
#print(CCmod)
} else if (koronaLiczenie == 2 ){
if (odl< ( testSite[[rKorony]][j]+dcent[[rKorony]][i]) &
!(is.na(testSite[[rKorony]][j])) &
!(is.na(dcent[[rKorony]][i]) ) ) {
Oir <- obliczCzWsp(dcent[[rKorony]][i], testSite[[rKorony]][j], odl)
Oirzecz <- Oirzecz + (Oir/(3.1415*dcent[[rKorony]][i]**2))
ME3rzecz <- ME3rzecz + (Oir/(3.1415*testSite[[rKorony]][j]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BRsklr <- (Oir/(3.1415*dcent[[rKorony]][i]**2))*( testSite[[d13]][j]/dcent[[d13]][i] )
BR1rzecz <- BR1rzecz + BRsklr**1
BR2rzecz <- BR2rzecz + BRsklr**1.5
BR3rzecz <- BR3rzecz + BRsklr**2
BR4rzecz <- BR4rzecz + BRsklr**2.5
BR5rzecz <- BR5rzecz + BRsklr**3
CCrzecz <- CCrzecz +(3.1415*testSite[[rKorony]][j]**2)
}
}
}
}
if (ilekonkurentow > 0) {
CI26 <- 0
# specjalna petla dla wskaznika CI26
for (j in 1:nrow(testSite)) {
odl=odlAB(dcent[[Xo]][i],
dcent[[Yo]][i],
testSite[[Xo]][j],
testSite[[Yo]][j])
if (odl > 0 & odl < pr_pow ) {
CI26 <- CI26 + 3.1415*((odl*dcent[[d13]][i])/
(dcent[[d13]][i]+testSite[[d13]][j])**2*
((testSite[[d13]][j]/odl)/(di/CI23)))
}
}
#CI36 obliczane jest w petli - PAMIETAC O DODANIU DO WYPISU
CI35 <- G(dcent[i,],powierzchnia=pow, typ='W') / (sCI35/nCI35 )
CI37 <- G(dcent[i,],powierzchnia=pow, typ='W') / (gCI37)
# sorutujemy drzewa po wysokosci na stanowisku
Hwyb <- testSite[order(-SDIsklad[[h]]),]
if (nrow(SDIsklad) < ile100( nrow(SDIsklad), pow/10000)) {
iledrz <- nrow(SDIsklad)
} else {
iledrz <- ile100( nrow(SDIsklad), pow/10000)
}
# liczymy wysokosc z drzew dominujacych
H100 <- mean(Hwyb[1:iledrz,h])
RS <- sqrt(pow/nrow(SDIsklad))/H100
BALMOD <- (1-(1-BAL/G(SDIsklad, powierzchnia=pow)))/RS
# obliczenie kwadratu przecietnej piersnicy d-stanu
SDIsklad <- as.data.frame(SDIsklad)
if (nrow(SDIsklad) > 0) {
zakd13 <- range(SDIsklad[[d13]])
pocz <- zakd13[1] %/% 2
kon <- (zakd13[2] %/% 2) +1
d_n <- 0
dlicz<- 0
# kwadrat przecietnej piersnicy
for (k in pocz:kon) {
dd <- subset(SDIsklad,
SDIsklad[[d13]] > (k*2 - 0.01) & SDIsklad[[d13]]< ((1+k)*2))
# jezeli w tym przedziale nie ma zadnych piersnic idziemy dalej
if (length(dd[[d13]]>0)) {
dlicz <- dlicz + length(dd[[d13]])*(1 + k * 2)
d_n <- d_n + length(dd[[d13]])
}
}
dg <- (dlicz/d_n) # kwadrat przecietnej piersnicy d-stanu
SDI <- nrow(SDIsklad)* (25 / dg)**( -1.605 )
# korzystamy z wybranych drzew konkurencyjnych do SDIsklad
CI28 <- G(dcent[i,],powierzchnia=pow, typ='W')/G(SDIsklad, powierzchnia=pow, typ='M')
if (nrow(SDIsklad)>0) {
CI29 <- G(dcent[i,],powierzchnia=pow, typ='W')/G(SDIsklad, powierzchnia=pow, typ='D')
} else {
CI29 <- 0
}
AP<-(G(dcent[i,], typ='W', powierzchnia=pow)/
(G(SDIsklad, powierzchnia=pow)*nrow(SDIsklad)))
} else {
SDI <- 0
CI28 <- 0
CI29 <- 0
AP <- 0
}
# wskazniki zalezne od odleglosci (krony drzew)
PK16mod <- Oimodel
PK16rzecz <- Oirzecz
}
l <- c(dcent[[nrDrzewa]][i], dcent[[nrPow]][i], CI,CI2 ,CI3,CI4 ,CI5 ,CI6,
CI7,CI8 ,CI9,CI10,CI11,CI12,CI13,CI14,CI15,CI16,CI17,CI18,
CI19,CI20,CI21,CI22,CI23,CI24,CI25,CI26,AP, CI28, CI29,
BAL,BALMOD, CI32, CI33, SDI, CI35, CI36, CI37,
ME3mod, BR1mod, BR2mod, BR3mod, BR4mod, BR5mod, PK16mod, CCmod,
ME3rzecz, BR1rzecz, BR2rzecz, BR3rzecz, BR4rzecz, BR5rzecz,
PK16rzecz, CCrzecz, ilekonkurentow )
# testowanie
#sp <- 0
#if (sp == 0) sp <- ll
#for (hh in 1:54) {
#print(paste(sp[hh], l[hh]))
#}
#print(length(l))
ll <- rbind(ll, l)
}
#if (length(brakh) > 0 ) {
#for (ii in 1:length(brakh)) {
#print(paste(c('Pow nr.', testSite[[nrPow]][1], 'drz nr.', brakh[ii], 'brak wysokosci - pomijam'), collapse = ' '))
#}
#}
# dodajemy nazwy kolumn i wierszay do dataframe'u
if (nrow(ll) > 0 ) {
names(ll) <- ll_names
nazTemp<-c(ll[,1], 'a')
rownames(ll) <- paste(dcent[[nrPow]][1], nazTemp[-length(nazTemp)], sep='_')
#usuwamy numer drzewa i nr powierzchni
ll[1:2] <-list( NULL)
}
return(ll)
}
}
dolaczAB <- function(A,B) {
# funkja laczy 2 data.frame'y w jeden sprawdzajac poprawnosc danych w
# kazdej kolumnie
if (is.data.frame(A) & nrow(A) >0 & is.data.frame(B) & nrow(B) >0) {
a.names <- names(A)
b.names <- names(B)
all.names <- union(a.names,b.names)
#print(paste("Number of columns:",length(all.names)))
a.type <- NULL
for (i in 1:ncol(A)) {
a.type[i] <- typeof(A[,i])
}
b.type <- NULL
for (i in 1:ncol(B)) {
b.type[i] <- typeof(B[,i])
}
a_b.names <- names(A)[!names(A)%in%names(B)]
b_a.names <- names(B)[!names(B)%in%names(A)]
if (length(a_b.names)>0 | length(b_a.names)>0){
print("Columns in data frame A but not in data frame B:")
print(length(a_b.names))
print("Columns in data frame B but not in data frame A:")
print(length(b_a.names))
} else if(a.names==b.names & a.type==b.type){
C <- rbind(A,B)
return(C)
}
C <- list()
for(i in 1:length(all.names)) {
l.a <- all.names[i]%in%a.names
pos.a <- match(all.names[i],a.names)
typ.a <- a.type[pos.a]
l.b <- all.names[i]%in%b.names
pos.b <- match(all.names[i],b.names)
typ.b <- b.type[pos.b]
if(l.a & l.b) {
if(typ.a==typ.b) {
vec <- c(A[,pos.a],B[,pos.b])
} else {
warning(c("Type mismatch in variable named: ",all.names[i],"\n"))
vec <- try(c(A[,pos.a],B[,pos.b]))
}
} else if (l.a) {
vec <- c(A[,pos.a],rep(NA,nrow(B)))
} else {
vec <- c(rep(NA,nrow(A)),B[,pos.b])
}
C[[i]] <- vec
}
names(C) <- all.names
C <- as.data.frame(C)
return(C)
} else if (is.data.frame(A) & nrow(A) >=0 & is.data.frame(B) & nrow(B) >=0) {
C <- rbind(A, B)
return(C)
} else {
print('Bledne typy danych do laczenia, oba zestawy danych musza byc w data.frame')
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.