<div dir="ltr">Boa tarde colegas programadores!<div><br></div><div>Estou tendo dificuldades em acertar o formato que devo organizar os outputs de um looping que fiz. </div><div>Tenho um banco de dados com datas e algumas variáveis climáticas. Dessas datas, preciso simular alguns cálculos para dias específicos. São 36 simulações por ano, durante um período de 34 anos, totalizando 1224 simulações.</div><div>Para alcançar esse objetivo estou trabalhando com 3 comandos for, um dentro do outro.</div><div>O primeiro seleciona as datas/posições que devem ser feitos os cálculos.</div><div>O segundo executa os devidos cálculos para as datas selecionadas e o ultimo integra os valores calculados, gerando a informação de interesse.</div><div><br></div><div>Acredito que a logica de cálculos dentro de cada for() está ok, mas a saída não está puxando os valores que o looping calcula.</div><div><br></div><div>Segue o script no corpo do email e, em anexo, um exemplo de um ano do banco de dados que está sendo utilizado.</div><div><br></div><div>Desde já, agradeço a colaboração dos colegas.</div><div><br></div><div>Abs</div><div><br></div><div>#--------------------Remover objetos salvos na memória do R--------------------#</div><div>rm(list = ls())</div><div><br></div><div>#--------------------Seleção de diretorios e arquivos de interesse--------------------#</div><div>a = read.table("x.txt", sep = "\t", header = T)<br></div><div><br></div><div>#--------------------Calculo das Variaves Astronomicas--------------------#</div><div>head(a)</div><div>a$DATE = as.Date(levels(a$DATE))[a$DATE]</div><div>ano = format(a$DATE, trim = T, '%Y'); a$ano=ano</div><div>mes = format(a$DATE, trim = T, '%m') ; a$mes=mes</div><div>dia = format(a$DATE, trim = T, '%d') ; a$dia=dia </div><div>dece = ifelse(dia<11, 1,ifelse(dia<21,2,3)); a$decendio=dece</div><div><br></div><div>lat = -7.53</div><div>corr = pi/180</div><div>decl = 23.45*sin(corr*((a$DiaJuliano-80)*360/365)); decl</div><div>hn = 1/corr*acos(-tan(corr*lat)*tan(corr*decl))</div><div>a$N = 2*hn/15</div><div>s = (1+0.033*cos(a$DiaJuliano*360/365))</div><div>t = corr*hn*sin(corr*lat)*sin(corr*decl)</div><div>u = cos(corr*lat)*cos(corr*decl)*sin(corr*hn)</div><div>a$Qo = 37.6*s*(t+u)</div><div>head(a)</div><div>x = 0.25</div><div>b = 0.50</div><div>n.est = (a$SRAD/a$Qo-x)*a$N/b; n.est</div><div>a$n<- ifelse((n.est)<0,1, n.est); head(a)</div><div>a$n.N = a$n/a$N; head(a)</div><div><br></div><div>#--------------------Transformar MJ em caloria--------------------#</div><div>a$Srad.cal = round((a$SRAD*1000000)/41800,1); head(a) #round é para delimitar o numero de casas apos a virgula</div><div>a$sem = paste0(dia, '-', mes)</div><div>a$TMED = (((a$TMAX+a$TMIN))/2)</div><div>head(a)</div><div><br></div><div>#-------------------- Determinacao das datas de simulacao--------------------#</div><div>ndc = 140</div><div>datas = 1:366 # referentes aos dia juliano</div><div>dim(a)</div><div>anos = table(a$ano); anos</div><div>potencial = matrix(NA, nrow = 366, ncol = 34 )</div><div>rownames(potencial) = 1:366</div><div>colnames(potencial) = c(1980:2013)</div><div><br></div><div>head(potencial); dim(potencial)</div><div>head(a)</div><div><br></div><div>for (i in (a$DiaJuliano)){</div><div> </div><div> epocas = a[which( a$sem == '01-01'|a$sem == '11-01' |a$sem =='21-01'|</div><div> a$sem == '01-02'|a$sem == '11-02' |a$sem =='21-02'|</div><div> a$sem == '01-03'|a$sem == '11-03' |a$sem =='21-03'|</div><div> a$sem == '01-04'|a$sem == '11-04' |a$sem =='21-04'|</div><div> a$sem == '01-05'|a$sem == '11-05' |a$sem =='21-05'|</div><div> a$sem == '01-06'|a$sem == '11-06' |a$sem =='21-06'|</div><div> a$sem == '01-07'|a$sem == '11-07' |a$sem =='21-07'|</div><div> a$sem == '01-08'|a$sem == '11-08' |a$sem =='21-08'|</div><div> a$sem == '01-09'|a$sem == '11-09' |a$sem =='21-09'|</div><div> a$sem == '01-10'|a$sem == '11-10' |a$sem =='21-10'|</div><div> a$sem == '01-11'|a$sem == '11-11' |a$sem =='21-11'|</div><div> a$sem == '01-12'|a$sem == '11-12' |a$sem =='21-12'& a$ano>=1980),]</div><div> </div><div> semear = length(epocas[,1])</div><div> </div><div> for (j in length(semear)){</div><div> </div><div> ini = as.numeric(rownames(epocas[j]))</div><div> fim = (ini+ndc)-1</div><div> </div><div> pac = a[(ini:fim),]</div><div> </div><div> #--------------------Inserindo as Produtividades Parciais--------------------#</div><div> #----------Inserindo influencia da nebulosidade----------#</div><div> pac$ctc[j] = -4.16+(0.4325*pac$TMED[j])-(0.00725*((a$TMED[j])^2))</div><div> pac$ctn[j] = -1.064+(0.173*a$TMED[j])-(0.0029*((a$TMED[j])^2))</div><div> pac$PPBc[j] = (107.2+(0.36*pac$Qo[j]*1000000)/41800)*pac$ctc[j]*pac$n.N[j]</div><div> pac$PPBn[j] = (31.7+(0.219*pac$Qo[j]*1000000)/41800)*pac$ctn[j]*(1-(pac$n.N[j]))</div><div> pac$PPB[j] = round((pac$PPBc[j])+(pac$PPBn[j]),1); head(pac)</div><div> </div><div> #--------------------Calculo da Produtividade - base seca --------------------#</div><div> pac$c.resp[j] = ifelse(a$TMED[j] < 20, 0.6, 0.5)</div><div> pac$c.colheita[j] = 0.5</div><div> umidade = 13</div><div> c.um = (1/(1-0.01*umidade))</div><div> pac$iaf[j] = 5</div><div> pac$c.iaf[j] = (0.0093+0.185*pac$iaf[j])-(0.0175*pac$iaf[j]^2)</div><div> pac$PP.Liquida_seca = NA</div><div> pac$PP.Liquida_seca[1] = round(pac$PPB[1]*pac$c.resp[1]*pac$c.colheita[1]*pac$c.iaf[1],1); head(pac) # em kg/ha</div><div> </div><div> for (w in 2:ndc){</div><div> </div><div> pac$PP.Liquida_seca[w] = round(pac$PPB[w]*pac$c.resp[w]*pac$c.colheita[w]*pac$c.iaf[w]+</div><div> pac$PP.Liquida_seca[w-1],1)</div><div> }</div><div> #--------------------Soma da PP--------------------#</div><div> pac$PP_final = round(pac$PP.Liquida_seca[w]*c.um,1)</div><div> dim(pac)</div><div> PP = pac[ndc,32]</div><div> potencial[i,j] = PP</div><div> }</div><div> print(i)</div><div>}</div><div><br></div><div>View(pac)<br></div><div>View(potencial)</div><div> </div><div><br></div><div><br clear="all"><div><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div>Yury Duarte<br></div>Engenheiro Agrônomo - ESALQ/USP<br></div></div></div>
</div></div>