Форумы GIS-Lab.info

Код: Выделить всё

fn="burned_dnestr.csv"
p=read.table(fn,header=TRUE,sep=";",dec=",",as.is=TRUE)
names(p)
str(p)
# [1] "month"      "DNBR_03.04" "DNBR_02.03" "DNBR_01.02"

old.par=par(mfrow=c(1,3))
hist(p$DNBR_03.04,breaks=100,col="gray")
hist(p$DNBR_02.03,breaks=100,col="gray")
hist(p$DNBR_01.02,breaks=100,col="gray")
par(old.par)

# https://www.r-bloggers.com/2011/08/fitting-mixture-distributions-with-the-r-package-mixtools/
library(mixtools)

#========================================================================================================================
# --- распределение
#========================================================================================================================
cur.x = p$DNBR_03.04
cur.mu = c(0,0.2)
cur.sigma = c(0.05,0.05)
cur.fit1 = normalmixEM(cur.x, lambda = NULL, mu = cur.mu, sigma = cur.sigma, k = 2, 
                      mean.constr = NULL, sd.constr = NULL,
                      epsilon = 1e-08, maxit = 1000, maxrestarts=20, 
                      verb = TRUE, fast=FALSE, ECM = FALSE,
                      arbmean = TRUE, arbvar = TRUE) 

hist(cur.x,breaks=100,col="gray",prob=TRUE)
xx=seq(-0.1,0.4,len=500)
for(i in 1:length(cur.fit1$lambda)) {
  ff=cur.fit1$lambda[i]*dnorm(xx,cur.fit1$mu[i],cur.fit1$sigma[i])
  points(xx,ff,type="l",col="red",lwd=3)
}

#-----------------------------------------------------------------------------------------
cur.x = p$DNBR_02.03
cur.mu = c(0,0.05,0.2)
cur.sigma = c(0.05,0.05,0.05)
cur.fit2 = normalmixEM(cur.x, lambda = NULL, mu = cur.mu, sigma = cur.sigma, k = 2, 
                      mean.constr = NULL, sd.constr = NULL,
                      epsilon = 1e-08, maxit = 1000, maxrestarts=20, 
                      verb = TRUE, fast=FALSE, ECM = FALSE,
                      arbmean = TRUE, arbvar = TRUE) 
hist(cur.x,breaks=100,col="gray",prob=TRUE)
xx=seq(-0.4,0.4,len=500)
for(i in 1:length(cur.fit2$lambda)) {
  ff=cur.fit2$lambda[i]*dnorm(xx,cur.fit2$mu[i],cur.fit2$sigma[i])
  points(xx,ff,type="l",col="red",lwd=3)
}

#-----------------------------------------------------------------------------------------
cur.x = p$DNBR_01.02
cur.mu = c(0,0.2)
cur.sigma = c(0.05,0.05)
cur.fit1 = normalmixEM(cur.x, lambda = NULL, mu = cur.mu, sigma = cur.sigma, k = 2, 
                      mean.constr = NULL, sd.constr = NULL,
                      epsilon = 1e-08, maxit = 1000, maxrestarts=20, 
                      verb = TRUE, fast=FALSE, ECM = FALSE,
                      arbmean = TRUE, arbvar = TRUE) 

hist(cur.x,breaks=100,col="gray",prob=TRUE)
xx=seq(-0.1,0.4,len=500)
for(i in 1:length(cur.fit1$lambda)) {
  ff=cur.fit1$lambda[i]*dnorm(xx,cur.fit1$mu[i],cur.fit1$sigma[i])
  points(xx,ff,type="l",col="red",lwd=3)
}

#========================================================================================================================
# --- Месяцы
#========================================================================================================================
table(p$month)
#  April     Feb Kontrol   March 
#  14855    3381    2211    8308

# Горела растительность в феврале, марте, апреле.
# Поля DNBR_03-04 - результат расчета по формуле определения выгоревших и контрольных участков по месяцам когда наблюдались пожары.
# Выгорел в феврале (DNBR_01-02) участок и надобно показать, что значения, которые рассчитаны для данного участка отличаются от других негоревших и горевших. 
# Т.е. в феврале сравниваем группу month- "Feb" с остальными.
# В марте -"March" (DNBR_02-03) и "April" ( DNBR_03-04) с остальными.
# Сравнение несгоревших участков между собой нужно для дальнейшей оценки развития растительности в зависимости от локальных условий.

cur.month = "Feb"
ind = (p$month == cur.month)
old.par=par(mfrow=c(1,3))
hist(p$DNBR_03.04[ind],breaks=100,col="gray")
hist(p$DNBR_02.03,breaks=100,col="gray")
hist(p$DNBR_01.02,breaks=100,col="gray")
par(old.par)

Код: Выделить всё

month.list=sort(unique(p$month))
var.list = c("DNBR_03.04", "DNBR_02.03", "DNBR_01.02")
old.par=par(mfrow=c(3,4))
for(i.var in 1:3) {
  cur.var = var.list[i.var]
  cur.x=p[,cur.var]
  for(i.month in 1:length(month.list)) {
    cur.month = month.list[i.month]
    cur.sel = p$month == cur.month
    hist(cur.x[cur.sel],breaks=50,col="gray",main=sprintf("%s, %s",cur.var,cur.month))
  }
}
par(old.par)