From 611bc170175d0d0617732a40ecea6a17c40c0966 Mon Sep 17 00:00:00 2001 From: Timon Date: Thu, 29 Aug 2024 16:35:59 +0200 Subject: [PATCH] adjusted all functions and cleaned scripts --- ...ots planting date and harevsting data.xlsx | Bin 96127 -> 99055 bytes laravel_app/tests/__fixtures__/harvest.xlsx | Bin 22525 -> 25453 bytes r_app/1_harvest_data_EcoFarm_v2.R | 195 ------------ r_app/ci_extraction.R | 166 ++++------- r_app/ci_extraction_utils.R | 277 ++---------------- r_app/interpolate_growth_model.R | 109 +++---- r_app/mosaic_creation.R | 100 ++++--- r_app/mosaic_creation_utils.R | 204 ------------- r_app/parameters_project.R | 137 ++------- 9 files changed, 227 insertions(+), 961 deletions(-) delete mode 100644 r_app/1_harvest_data_EcoFarm_v2.R diff --git a/Current - Pivots planting date and harevsting data.xlsx b/Current - Pivots planting date and harevsting data.xlsx index 6dc445e765684a7a52dc958222a4469227939d67..f5445d7a9915bc63cb6f3f1b8e7fe6002e0236f0 100644 GIT binary patch delta 3334 zcmb7HWmFVgx1J$ofMJvzrAE3_K@^nkR=R{ChAwFs1{hjW7#Knr@>L{;80rNC=@11e zDJ7)4q`bV}z2EwN+;xB4=d69s-e<49_d4g-^ZeN%-9?jy*8!=lgh7D+1`P1$2IV*5 zKL)wEH23v>>EmQ6A}T5>CJcY+`%j)AuMT7Z67r zriTW12sd8;jWGfEF`j-vIt&d0NKVQP;4NUs+)U%ZJaGZiASaAm#5776S_J~A*C!aT zmI&OfpY1tgqK7u~*Jobo0ZUt|+|=v*Koq>`-g@B=KR0`dPIAS(>?}QEF=m>>-BT3aE%Q_xLW&5XNS69l024Bxrh|tTHE$t*iZb6dmX2w?Iax z-aRvrri!~mo^Z7rtQLRgRSie~v~tZ9gPWb{uf4scCbpAd-QH06ldqDL=>ZZoRTXja zcDX5{VN+V(3|a5{bv+hr8hX@+oQ4`clLsaH`zBKbob@`0i=C8T8V9cel+!9p^>|-CPy_9+jl8M>3A`ldEa>u&-HrIb6L_UvCTg_g(0sxVCxqd2dLScxFsRRL-emKc-w*2 z9H@~OmYM0}&iLUizsSM=NS3st?X zcATq!x+jlhM3<1Q?_g0M??!OELFSv=+KZP4UtPW zN&&cacL9c2^ZRWbXn?*nHPjN^X#EDX#QuoAVadzu(N(+wS*aH>)C?Q$quGg$Z;vrvWB~`vq`tm6+Yck6+z5yLj*~FEzLJhI<3eip9#Y5hpq6)+}KT{Nm8>n-< zOMS4Nn6rXVx|>>|)O{~nBViZ-H>+Eqk>QDSh)b8fM>Y;kVo5mrp|57(1Iv zRFQrt{IMD{zv!HfC{%{q3vWEpsoB$Ol~E z!NZPBnv#fQi!mBk=bsrXiZTopkRREMNVCYv=>@|cgrJ9gOP2I%NDa>6ne5=Q?hw7{ zkSjaIC4C2#ByUd(=L#n(xu7pXhUD=$Sz3l0ZpHhh?eaP0XFua9w8p@Nz&ru`I6pn>*bQ{0&VpI3QMu^dERQ-Ub$Ge<PQw&1x!>-}vj{{CA7z0%1b8 zJSl^|E*(hru=?(M%G5w=|8e120r+`orR9@zUTwxXbjl?IG9kUpB*FbSc3LnDUcLEE zWbNDPC9-`H*=9ZJK9ht7PbQPfe-(M2 z*b_oj_UcwZjOjR?{o`NHuZ0Y4uCH0@(9ht1i@f7m4BX;*RHvoZk5}wGd_QZXoe}rZ z+&CwWm%5<)_EX7Bxi7;_Jxs;n3aVm@R}h$Nm7X(RifAb?x^*3Al8;eurgO7}u}W62 zohWMK4weU6TSGfT8wr|7tj6be+43blvI@A;9x}dOaIP{zW^*uCZ>_m(Roo>buC(Po zOmE=Iy4-DpHV#zsdoYeTC4M9drgE2TNHQ^pbKtbKJ$d2-6BG@NJPFq{&L?D%mhKAj zC0;?^{(4EoM(huTASn!8r&-I3=IA89bLrx8YO>2yL7X*e6<0{5D^z4xG}+-H!lh`m z4vgdo86Fl!eVpm#{GG|B!Q(a0aXCV1WbZ%whH}HbcdVA6>D$HPi=SPG3w*2ePZbBY&V7x2N!dn)ZbQNby1{xsFVk3N_`Dyy@ycKvO&BGZJl z;b5LwE#F6Dn^^nd|w{@ z>9AFdsRo;(l?v6ek(oU+kCugyp6)o4iL9i~<;+f}RzoGVQ+*UOZph}b^QYhQWq+lT zs6?~=FipFw(YG}))p>8tGgcZ{6&6mx`6h#F>zX0u4||K7anD(A+rciYD_{0NEE1O|LGzX&=MR$WMW*s)37O)P79+CE z7IFh}pGK^kk70^Q&l3`MX!F|vd!WD4hQzFX-B(VnK5BL_u_?;MKgq4~)TA1!Inzl< z0U^~taTqA8MCGRjvZa#8DELP zzPL!iew_`4mTnFhrPVJw;q}oThL+YpM^VFSrdYjWxZ$Z`#wz_|cP7z=`9Nc?)DO^u zK!X0dyP{M^=cNWQ_TX{%>+*mRw(*8=KGR6B^uE#T7}RAgxA{-*j}oQ_ z3V$?!!qcc&{;qAUin%SpI;tc3#6%HJduIQ(Z@wka?7AmfeWncxWtM;XH@Dk?1Zywg zeG={eTUp$&{r|x&E8zdq8G4#PA|b$kNc6u%0^!RmAQp(3^aFxb|HG;Q0NQ_d`qydi q1q9f=^6~X>GjVejLHOFci3ta|x$0?>lKqosx>*W09wr2M6a5#B{2crM delta 3076 zcmV+f4Eyu%g$Dob1%R{x3|5-GdJ zuIj4ldU>jfy!qw0kF-acW?>R{h&sku&`Y&(aV7A6d%xKKE&j)&s zcXW?5c=+c=*R?$=|E&G!e>2+CjeeWpd=F&eSG8~NZa2DlG%fhoT&>H%;@}w7qM`^z7xkHtip!w08TeTbkiCJTRfmdm3wTNb@I+6jFew-)MZ*rn<+; z0b~`A3sZ3UsJU7ompY#3EKMiyUC$xz!kK$3x*+rULf61{akj>6f1~dhJFI83Ojstn z+J`tzvLs+vIN8sVAP6y?fm+Mio-#xbL$jJO6r(+w!sw8K^eGEzCYWC{mWF=LXr{mY z`KLF>jI5U1Hz#2c%iGagt5{j0o_@eQ}HnfX+$z^x}`;k>CCY`w_QqC zD`&y2N9PN3Wz!{{f9qPD&864GnJxb!jS2FWf;~MwT|L=XN&3KTXU+cc>7#T#Ww2v5 z{I`elrz0+@uCJ62lu#a+lGf+WH^gp9t5s!BTSTxpdYxpPnVYI zc!p1!e+T$)k``HTII9Q;t7Tv<7MK!;5bQ&CU0Jph1U9zqc`#U;Ejr#HhIiLd zG?1pH$Ib7P=zgU!O%c!S+i=e_&iGq@{|Ju#SCs8;5;Bd$En;-btH`AD?rVC<6(>W!9e) zMrFN_LDV2?#~BaH(-vvflhP{QIIDEnjId0H_Q$c#$NxBg!`P~*kdwqPztzm zfP05;v_*#?yF%9u{CF&R20ak5v@f5IELW}FQ=Lvi`Q(SA2R@Ad%kw&;Na z6SQuG|8?L~pn=wssqZVB1}MsEZNwd1*S{78WHC5qF*s;3IBGHKu;t(493|8zC_G{9 zqOGKn?3+`A325}Sk%C-=>sSVTtqnjL2Jl5`+Lebjp@0TgNbb{mM%2{ zf9jbAe`2^~Ik$iW8bUBL(4ZLal9ZDKK3YJdXdBq37Ld4>fmW-ffsx}cJSxbEiF` zD2Ey`Uo2`7TdN_Q;>9JSgJ&yitu9-uZj2sIJxTCXMeU(}p;cd)>guYs5qo8$k$9ua zXoQKn#6OOO>K#dCosqbs$)v|J!z*Z$)sJJRp5Ci8*Br*4W6(!0uO6$N3b|hTe;1l+ z6quFCYf%3;Xn6Gh{_lUj|LYFF%J8uUgwL;(aPX(W09-(zW@+zZej*-vFj4V{x85+{ zQRN))HjVjSqt@~9xt+j>WqW=KPf19_e}4yUK8ReOHq`-D9{WiSEifrr_DFdCl-sIF z)-^tg*ysZ4LCYFHw2nT}onK=Ne^LN&RSbB4E2QU$vmtvc(Vb_=n`p)Iyi78s)aJnA zLCY56*+t}8-du^Zx;H$8dl#ZsepXFlD>t?-SAm)vYor|Jau)Sm-+#;rhyGEqYIDw8 zS!r?ENE>=Ia4 zd`yyW`Jqe_=QF||ya7>h$3wZG?oT`e0=z#uUn|fg_d#%?HT>kAFic$om>VOOiCrHo zC%?fAp|dYdsZm=-T|yj4f9*jVBk4S@4W<1k!QaYh^HcB6L!L&$k_69+T$|0TD{EH$7LrwMWL_r_-;}>hU<6$X(?3~P zu{cZossx^#OV^KT$tKxmOO{Z?=E7z9@@wE&CKKx<6Kcr&t=*tDcw& z5xmAwdthnUh60VbPWrCR$&{Jl;!1og?Dct9E$M( z@|vQPiKV9F$=}`HTxVH=Lr8h#JwONcJo&nTSHaUYKUoTD2|a{rt$y1y&XKRM&zHN8C#G%JF6J=?eY9O0to~k^r(6|-AT)p$DN?NuN$mi2W<-0B4R1Why z4QR^Wf7EPM#dGE(9>#$;2aAn=;UWLoch&r=Vl>U0kAV=ZHa3)by2}n;xGw zC6k`3p3I8h^>kvm?CHdC+tZ2Rx~CJveUB&JiuW$220Eu6zGt{pqbB%YDZooK5LIP- zoQbpBlZ8%zzCw`kgQSep?zh^mv@25$hY2~$BT?wA?|+?n`yT)R0RR6308mQ<1QY-U z00;m803iSmY&X2Qm+@)=9+wYk0SA|uaRDZO42{8x!Y~v?_XGYzc>9<(aXMC-irNb9 zXTWVpU$xM@CruP#n17$FF83a|r`f;5EMB<8iOE8_E?I z0AqD?bZ>1~ZEP=TbY*QZE_iKhP)h*EP0H6YJ0Ic4) zxjQ)+I@k*$-EAC%1iT#VS+IX4+?8idcp}V&n#uoWvh~QopnMKf=yFH`HQts~tYR!} z;Vv(ga$ij)FZ){H^vjkqAT*Zm%XqwVWtx;|o~U?;;=>j;D={M)zXJ}NUY#%Zt!QH3 zWTjXKvEeo2{d9*Z+@g=D+$A${RDFIDrk06KTTI#I_V#9)Biwa`wug|jamUQvD$RB@ z${>yt5+7V~C_pT$JXSkXyhtRlnU_FXXX-D-gGdUcMH39*^?kDznrXPd#uZOP-Lf}t0s~X$?sv&~ z1vNK;iX71{(_~^=%(!7BgI6mL5N}GYA6?_(u>l9EJh6_utfoG_qP>(#_$o(J^EU9+4bC*9v}!z;Mb=oV!`qP6LcT4Rl-SuDDwBMpf(s-aXaW5e8hxs zNPHTlR~qC=Tn?SH`!Jov5y5iR_=rN}&WeLaHY`j$y7DCsmQBMF`?ywmj^9*}6Zz!J zQqdbiDr1(=IS(XKy=0NU`O1dS9C>#Ge zG;r+bPC_6BHBOP2G9jFT9Q2%(Wpw)Hlb~V=4&n+Z6Q^$yv;oW(SHTPZ`g}%vO3HI` z>D5)O2ch9FJ?2TM)jDQ+<5jTtNFy(|eN`jwy{TmFnO3TtqIIw``n9vOSkHVkpws1V z{dh!_EWt~pvb-c*&I%I?3z*V$rA;mF(Q%r$sOePef7M?zMCudc=^jJwb<$-kEObP zCKv5F@En7Tr@=&2OZ>#XK0XaTu4VBDfmSbduCLZ|pGp%?39VmmeeiqJ0$n>#+6$b` z2w&u<+eiOOwQkoTH%2+dT-2aLSF`L>@bX*o8n$Eh!1Sff zJtw8PNAHm6-bK*yJ2o3-(+5w>4-0{hmXv^oy4_IyK>QZthD24)v184uS1|o|GZkH(^iCPafCoy)BJYtA zbx!l4zFRJ#*=Oe}@Qu${d(Q-CPSK!teoObh+-|Plz7cv&=f*x8x|%W^hB>gP6V~_4 zQ(FSfga9YPSTY=J?@fmE-AUxFY7MegKcmp{uVlbod?H-24b& z`EE#_tQ$2N9IY;R!A!<$ZQhZSLa6cWt~SI=SAqg!Mp2mU+_1}d3 z5(7hWgx7~|=L+kov3`nsxEYnPj8tTZD^cvYAF3WX0I)W$p1&i>8EhM#Abp=?3=;hy zHZVx9`;~9eB&3d7 z`dKU~)3@OU+pEbyoj{TZh|63%+hzGrDIzGt5up_6*WKA|h` zAJjr+N1ACvI*7!Bz_=W-mgGD6K{2MIckJzcCM_#S(w2aKq|u>_gD0ow5uHdrC!5Ap zi50&}tnCx&z9pT08l!%DmfIIJE##tHosDeEY=JR3-9eH>Qe!0XNeBg7x+9-K%cB5D0rHj3ydkCVjZ`;lu#M7G(ky=OWxtUR z%TCq^*RZ=)r>0eg{GrYfkqXMZ31!2^%~8pe%A|#qoL)Rw54d^QZ*-aVym6I&5sV?# zsTxA@wn*#Zk|l3ao<@>zKUy+@()vtgSh9GVv@rTpc?jw?3B@F|=QU?J<(PJwDd!~` z{OkkGR`mCVL?C}|w7Tp{=_JS4_ znNVmL>3A46%@<%@u|6)iI=*s&`n-T@HXn8zO{TRvcCM90v_%pVDnCT;4`(9j_Nh3& z#wpOnc%h0(B;$B4gU_Nm{Rj%*-12&7G)8SB_3P;sKf>bb>OnR7iS^%tpEwu1H#i?v zYpV9(6x#NSXZ5v`!dp!YGQzniaystHiYCi^9jNc5&wDGcBD8P`hDw*i?Qr4*i{7D4 zYgogq6txCw2TLfEXvON0f);jfskf=gzs;*w7sRV zbRSQy<+`1uEtG$N*riMTN=+3!QD@PVn7gt)yvt)DiF#Bb-Y#J^{E>h?P9xT?+ogL{l|^l*dKF}Y zIPcBU-wvu8QRb5jJSp%K35-+X);wlYp8hGiTMYJcjs74w9;<~jK+D(j^jva<^oelO z`0|9Lwv|<_WmhQBPe>s8GL5UaTMaD2gm16VPHTYhKPC;I?|kZmGYjoFtOlz?P#aD) zDA8uFVBCs-L^0uJ$}>%_(5N;#npz4)y10dhp;Ls5o3=dgb=Kt9JaJ@}mE1jnF;~BP z6a^{Q!3A?v(_oF8aUl*UwlbJVj4EE8(-0HaV$ZJao z*Rqo-!?w0>G{c1J&ZIdWc>c~vomt_zGWn9eu?-`qrYS#MXN=xZT2%XDW~)uJri|hh zK7Y9`^-HRVi<)HyC&YhS20*5Nek^}lC(6VtU}e7Ny&-w^mWU$xlP%+9CsWmS@G*@ z7ey0@vM{)AqWy=X<~zJMc*nuW_-$^;{X|$^5}BXAr^d2sjZe7*D$C9NZ?&y#5#?h~ zrK6#g4XsZs=p5&WvXt~^=dtM{v;L5R^d%fsUaw*rt}#FI=^x5ps{;k51H*XQH#JM?t^p+A2mCYXBS*xl9_xh>0Ob-oo?h$ycbK>#f&bMr@&81N ze=qwFKvY5Whk*$^2uhWIrtUw~_kYE$2lzJz6#2t|a1COTKWTI~>w@^Eso9_CUw5Rk A2mk;8 delta 3526 zcmV;%4LS1d#sU4W0kHWo20sA+080XsK~fZxurn=_J~AkO)DGR=a?iQCdhO;sjjJ~Z z(R$Y*n&oH(PI+zJp+m=HaE&&7J2!LiGlXOzu9k$0&D^0OrMV#_ssSb$v&pL9qc;KzviR(o}vz@?pZ@&K2sZ%1XS($Jz%pt;s+&vQ{x#*JZ$ zv9yO*U5&s?_{{x|ZrbvpEL`$TbOXJ_a6oAsL?PM#}0(I|ET ze>cP(_O4v6DF?q&(IWpnn!+c%hPao{uRv%^bBAMlZSNP>sI?PfE{L&d+9mwrNj z9^K9M0Zu$16K_?%z1#0*>Pgk$kGfvJTkiGS#a{iexG@cFvrtv_M%6brhPu&L@0zr~ z7SbBcul6*?krOhpnfKI}e1|5_63L_hQNK`s-lTfKK?t&n$C)a)eD-2vVxpo2U7M?o zsgjjxtcYe4ZLwB$Z8MWVcYeIZY&Sc9G8S$XB^H{A=bnS3AP#Ic$APB=w(Vf5z_wP3 zM;RiBA*G^>oRNp7FbX4(K4K1y1@jxmBF9P?jc0Ga|L)BxCfighiJ|C`GJY;`h%jt; zywK%ZLRnD|z}{JuQ1Q`rs7qq5I#X2*Lhx#>Ypa#Mwr#pW3;0oO3tLnNBlF{bMQSyE ztf#-E$^>}}!JeO==g;~)h#tA_H02-nA5+uI6?Te;zZ+h=-&8)3dpFZJ>(z3#(PgqU zH?pBwy1YRql8vRkux(S<5hl5eus<+$i>@@)*vK1AQ)L6Kmh#5XO_^*=YAp5jhMK52 zF2@Ulh)I5VP=nS}^5bsP4_dx|prKpF90?jy7a_DL1AL_YXWhj_+00a z``OONNusz5BgmDsqIjFPNK01j% zpqPzhpl=_JD5B)=4m?e)sm zmPsr$d+t#tJFYLj)IiS-mV)%IDx|3A@>GtPIbU#v`Mu^AugX8gtGyc@y zKY?NY6~)KffK21CM~ohL7MZmE14R)&el8;RBqkpNEVR%?tnhfBIgt3h;a}oe=2hIu z#OJu;yu*n^s=Jh8hDFMMOrQ@z6bw)iIIfE<7wVFtpdd~hy67b>C3{&M9YQD+uLYR{m1^XFBliMbL=pMqyDu4K zVakv{U0O|U5sjgy!aVXw7K%YaFq!R70i$WMNu8)h*7P&(mi3%i-A*GJ6e&r#J5vdeRo-6m^DUhswY}DHNpx%(S2bu`bN?ILqCf+EDa*evb!sCwuJj;x za9qDHvwBD71&_=t9+{Usa<%69mslqk>K+o0Q#t534ml^&v=LIFyvM2q)xL!c*|Ie% zV^g)DqOwgyc=3d`5Our<$~cX>F*O4dU0$wCDnt9D%KV9cY?9SN0}|UZ)8-ZV_!vZ- z#J14Vu%VT}I@N&0)MT_?uVjo2Yqik~6C02VVHs%mYSAwK8+~~WBF;jnVM3WMO5_l4 z53qw~MF~ru^6FZfy-iDoYfFgp@i46rON(4xcJKEm79q?ik%X(X&9fSBbN;Up<8p0- zmOfH2&ZG2ygxmyj#A3Oucx)wyly%n?DWg>km)cs73x*w>M&ykS>i=x0tDbYS+%I)B zvfgMj8lj>#@wcO)YCBVEXJp+`XL3~rt+V=e^wiO~73x}WIBP`l+R1B#J$2H5JLOuX zJuOwFT~<+zj#PnMk@ZBQFt84*vHEu!QQ8R%>Wl|}VRF*{JLCWL&bW;!mT`DMVEK^} zUKUF*SIFF?+WU4f@9~EYAU1o%K?#7(3QYC{pZK__(hGAqR8-I9%;78`%-_EQ0Hgt7 zRZGppw}J$~0dfv1t|C|Ml5$_d2Q!I-G`pD?hd?v$KVE{~IQ+hw-Czt-KuMo7;GMsa zp4Z}kRRFS?0(yXuj3I{?hRJjVLf3@#podIPuBbrUCp-lIYiqi+vl8Gg&3HiYub8n! z3XHtjVsY5ukFIP$dxG6mJFq&R-Z4OB)M8%%6!G{97j+L*lLl{fwifSn8LLqn#A-^1~ z4w53@@M|nDc-3X7gXw%5)u~#V>~9Z$zcs!nl5!>bWOMHJ%Xz`mH;*QGiy`;GEcQ$p zqxH|;KsmE~SLsxwEhtf&NbQ|} z3=$IaGw9v?mO0FAPGIXuV+!Gk=}-!l*5KfS!>HTOK2F52e41J-&F44~hrWFpKDu!`q0F+zPMs{eX!%c>W|G zc;*kKo}a_fiBphP0#L!>DvBj`q>`U_7LBC9mSXe-&w@zGg^cIY0p37i=yC$y{QwgH zb;AI9A(&c9kZPqQgsWZ9TL8CzVhLn}qe~RwW5I0x{dZlWrO0}77wJ<}+wJA)qRX$h zz%4#LK!11j>5K6|9|n_d|Pj}i*a^7JF0ctwR9q8M)~`VrE}_kN~1lMHtng; zQ-l4~_Q5~5ugLVsng+p6RAI|+2mq*XZmT_r0+(ncaUDGvPX2CxdlSb2b|B=D_fR^p zm&wNsyfO~`rAr0F;bE+**5*OLw}rnzpoHswFpf+YC~mA~+AZG8Oz9lmV3LQIm{IR8 ztN!_tgGqChLQ?3X8)RmGKm6kU( zRHd4|E?}iRa)bGq+BD)9U>jL+Na)Kr5rD}?hrXN^F2GEt)vuuCm3b8f9kX8my#DlH z-Eq=65Ij6=_8#?QR{X4|6T@LoCx+9WP7KFAofyu0JaKv>WvZh#bL!!9hD$kWLa|T? zaLo=xRvGVSV(lif&~ngMsEGJNsf^Rkx16q=sfC1`#bj0K`X5Pe{{WML4HOEb|HK-+ z0RRAf1d|;$8Ix~e1(UTiD1Wq+QES^U6ouah`wvF%vg}y4rP8D{^#X+z%GN<&75O?5 zvm_(Q%f{G$pX_#wt|{5r_GGTjJ?Gw|Pd>l+to!5>IIpZ(qCCw|0!G?ancWh7xZ7Zc z5+8(VMQ;r(k+JAv_4?I$>(5#dgb&uibpV~j7b-rln7khOhymo}{^srbIkNvU+>4?GudA}+1#U~m zG0nFH=7I@K%WlzitSE$RUu@2{P~{ornCE$hsi;cKX~A&I8PKw*IWY0}>$H!c+z2Cf zFe4MS_=n%oK(47jU^h_n4~D<5Hd5IOQ=obnla4)9M@PGvHT zuAj+hC4qik*^@lY$)7X^aN3{jI3ALUSH&686T6*HpZszqr~1h+WAz=AfejS0L@N%n zU0s<10Y8)XU;mTrUJ5NC6kMPx3IG7kGynh-000000000000000097{t0AqD?bZ>1~ zZEP=TbY*Qblb~M|3Z(zU8oU7j0DT0rNMX% #empty columns - rename(pivot_quadrant = 1, - area = 2, - variety = 3, - planting_date = 4, - Age = 5, - ratoon = 6, - Year_replanted = 7, - Harvesting_date_2020 = 8, - Harvesting_age_2020 = 9, - MT_weight_2020 = 10, - Tcha_2020 = 11, - Tchm_2020 = 12, - Harvesting_date_2021 = 13, - Harvesting_age_2021 = 14, - MT_weight_2021 = 15, - Tcha_2021 = 16, - Tchm_2021 = 17, - Harvesting_date_2022 = 18, - Harvesting_age_2022 = 19, - MT_weight_2022 = 20, - Tcha_2022 = 21, - Tchm_2022 = 22, - Harvesting_date_2023 = 23, - Harvesting_age_2023 = 24, - MT_weight_2023 = 25, - Tcha_2023 = 26, - Tchm_2023 = 27, - ) %>% - - # slice(-1) %>% #select(-age) %>% - # filter(pivot_quadrant != "Total") %>% #drop_na(pivot_quadrant) %>% - mutate(planting_date = ymd(planting_date ), - Harvesting_date_2020 = ymd(Harvesting_date_2020), - Harvesting_date_2021 = ymd(Harvesting_date_2021), - Harvesting_date_2022= ymd(Harvesting_date_2022), - Age = round(Age,0)) %>% filter(pivot_quadrant != "Total/Average")%>% - filter(pivot_quadrant != "Total") - -#copy each row and add ABCD -quadrants <- dates %>% slice(rep(1:n(), each=4)) %>% - group_by(pivot_quadrant) %>% - mutate(pivot_quadrant = paste0(pivot_quadrant, c("A", "B", "C", "D"))) %>% - filter(pivot_quadrant != "P1.8D" & pivot_quadrant != "P1.8 Q.DA"& pivot_quadrant != "P1.8 Q.DB"& pivot_quadrant != "P1.8 Q.DC") %>% - mutate(pivot_quadrant = case_when(pivot_quadrant == "P1.3 ABA" ~ "P1.3A", - pivot_quadrant == "P1.3 ABB" ~ "1.3B", - pivot_quadrant == "P1.3 ABC" ~ "1.3A", - pivot_quadrant == "P1.3 ABD" ~ "1.3A", - pivot_quadrant == "P1.3 CDA" ~ "1.3C", - pivot_quadrant == "P1.3 CDB" ~ "1.3D", - pivot_quadrant == "P1.3 CDC" ~ "1.3C", - pivot_quadrant == "P1.3 CDD" ~ "1.3C", - - pivot_quadrant == "P1.8 Q.DD" ~ "1.8D", - - pivot_quadrant == "P1.9.ABA" ~ "1.9A", - pivot_quadrant == "P1.9.ABB" ~ "1.9B", - pivot_quadrant == "P1.9.ABC" ~ "1.9A", - pivot_quadrant == "P1.9.ABD" ~ "1.9A", - pivot_quadrant == "P1.9.CDA" ~ "1.9C", - pivot_quadrant == "P1.9.CDB" ~ "1.9D", - pivot_quadrant == "P1.9.CDC" ~ "1.9C", - pivot_quadrant == "P1.9.CDD" ~ "1.9C", - - pivot_quadrant == "P2.3AA" ~ "2.3A", - pivot_quadrant == "P2.3AB" ~ "2.3A", - pivot_quadrant == "P2.3AC" ~ "2.3A", - pivot_quadrant == "P2.3AD" ~ "2.3A", - pivot_quadrant == "P2.3DA" ~ "2.3D", - pivot_quadrant == "P2.3DB" ~ "2.3D", - pivot_quadrant == "P2.3DC" ~ "2.3D", - pivot_quadrant == "P2.3DD" ~ "2.3D", - pivot_quadrant == "P2.3BCA" ~ "2.3B", - pivot_quadrant == "P2.3BCB" ~ "2.3B", - pivot_quadrant == "P2.3BCC" ~ "2.3C", - pivot_quadrant == "P2.3BCD" ~ "2.3C", - - pivot_quadrant == "P2.5 ABA" ~ "2.5A", - pivot_quadrant == "P2.5 ABB" ~ "2.5B", - pivot_quadrant == "P2.5 ABC" ~ "2.5A", - pivot_quadrant == "P2.5 ABD" ~ "2.5A", - pivot_quadrant == "P2.5 CDA" ~ "2.5C", - pivot_quadrant == "P2.5 CDB" ~ "2.5D", - pivot_quadrant == "P2.5 CDC" ~ "2.5C", - pivot_quadrant == "P2.5 CDD" ~ "2.5C", - - pivot_quadrant == "P3.1ABA" ~ "3.1A", - pivot_quadrant == "P3.1ABB" ~ "3.1B", - pivot_quadrant == "P3.1ABC" ~ "3.1A", - pivot_quadrant == "P3.1ABD" ~ "3.1A", - pivot_quadrant == "P3.1CDA" ~ "3.1C", - pivot_quadrant == "P3.1CDB" ~ "3.1D", - pivot_quadrant == "P3.1CDC" ~ "3.1C", - pivot_quadrant == "P3.1CDD" ~ "3.1C", - - pivot_quadrant == "P3.2 ABA" ~ "3.2A", - pivot_quadrant == "P3.2 ABB" ~ "3.2B", - pivot_quadrant == "P3.2 ABC" ~ "3.2A", - pivot_quadrant == "P3.2 ABD" ~ "3.2A", - pivot_quadrant == "P3.2 CDA" ~ "3.2C", - pivot_quadrant == "P3.2 CDB" ~ "3.2D", - pivot_quadrant == "P3.2 CDC" ~ "3.2C", - pivot_quadrant == "P3.2 CDD" ~ "3.2C", - - pivot_quadrant == "DL 1.3A" ~ "DL1.3", - pivot_quadrant == "DL 1.3B" ~ "DL1.3", - pivot_quadrant == "DL 1.3C" ~ "DL1.3", - pivot_quadrant == "DL 1.3D" ~ "DL1.3", - - pivot_quadrant == "DL 1.1A" ~ "DL1.1", - pivot_quadrant == "DL 1.1B" ~ "DL1.1", - pivot_quadrant == "DL 1.1C" ~ "DL1.1", - pivot_quadrant == "DL 1.1D" ~ "DL1.1", - - TRUE ~ pivot_quadrant) ) %>% unique() %>% - mutate_at("pivot_quadrant", str_replace, "P", "") %>% - mutate(pivot = pivot_quadrant) %>% - mutate_at("pivot",str_replace, "[ABCD]", "") %>% - mutate(pivot = case_when(pivot == "L1.1"~"DL1.1", - pivot == "L1.3" ~"DL1.3", - TRUE ~ pivot)) - -quadrants2 <- quadrants %>% - mutate( - season_start_2021 = case_when(!is.na(Harvesting_date_2021) ~ Harvesting_date_2021 - (Harvesting_age_2021 * 30)) , - season_end_2021 = case_when(!is.na(Harvesting_date_2021) ~ Harvesting_date_2021), - - season_start_2022 = case_when(is.na(Harvesting_date_2021) & !is.na(Harvesting_date_2022) ~ Harvesting_date_2022 - (Harvesting_age_2022 * 30), - !is.na(Harvesting_date_2021) & !is.na(Harvesting_date_2022) ~ Harvesting_date_2021 - ), - season_end_2022 = case_when(!is.na(Harvesting_date_2022) & !is.na(season_start_2022) ~ Harvesting_date_2022), - - season_start_2023 = case_when(ratoon == 0 ~ planting_date, - TRUE ~ Harvesting_date_2022), - season_start_2023 = case_when(is.na(Harvesting_date_2022) ~ Harvesting_date_2021, - TRUE ~ season_start_2023), - - season_end_2023 = case_when(!is.na(Harvesting_date_2023) ~ Harvesting_date_2023, - TRUE ~ ymd(Sys.Date())), - - season_start_2024 = case_when(!is.na(Harvesting_date_2023) ~ Harvesting_date_2023, - TRUE ~ NA), - season_end_2024 = case_when(!is.na(season_start_2024) ~ ymd(Sys.Date())), - ) - -saveRDS(quadrants2, here(harvest_dir, "harvest_data_new")) \ No newline at end of file diff --git a/r_app/ci_extraction.R b/r_app/ci_extraction.R index 6e210a1..5d37554 100644 --- a/r_app/ci_extraction.R +++ b/r_app/ci_extraction.R @@ -7,7 +7,6 @@ library(lubridate) library(exactextractr) library(readxl) - # Vang alle command line argumenten op args <- commandArgs(trailingOnly = TRUE) @@ -54,38 +53,31 @@ weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic") daily_vrt <- here(data_dir, "vrt") harvest_dir <- here(data_dir, "HarvestData") +source(here("r_app/parameters_project.R")) +source(here("r_app/ci_extraction_utils.R")) +# source(here("r_app/mosaic_creation_utils.R")) + source("parameters_project.R") source("ci_extraction_utils.R") -source("mosaic_creation_utils.R") +# source("mosaic_creation_utils.R") -dir.create(here(laravel_storage_dir)) -dir.create(here(data_dir)) -dir.create(here(extracted_CI_dir)) -dir.create(here(daily_CI_vals_dir)) -dir.create(here(cumulative_CI_vals_dir)) -dir.create(here(weekly_CI_mosaic)) -dir.create(here(daily_vrt)) -dir.create(merged_final) -dir.create(harvest_dir) +# dir.create(here(laravel_storage_dir)) +# dir.create(here(data_dir)) +# dir.create(here(extracted_CI_dir)) +# dir.create(here(daily_CI_vals_dir)) +# dir.create(here(cumulative_CI_vals_dir)) +# dir.create(here(weekly_CI_mosaic)) +# dir.create(here(daily_vrt)) +# dir.create(merged_final) +# dir.create(harvest_dir) -# end_date <- lubridate::dmy("20-6-2024") +end_date <- lubridate::dmy("28-08-2024") week <- week(end_date) -#weeks_ago = 0 -# Creating weekly mosaic -#dates <- date_list(weeks_ago) dates <- date_list(end_date, offset) print(dates) -#load pivot geojson -# pivot_sf_q <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect() -# raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif") - -# filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>% -# compact() %>% -# flatten_chr() -# head(filtered_files) raster_files <- list.files(planet_tif_folder,full.names = T, pattern = ".tif") @@ -93,39 +85,13 @@ filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = compact() %>% flatten_chr() head(filtered_files) -# filtered_files <- raster_files #for first CI extraction -# create_mask_and_crop <- function(file, field_boundaries) { -# message("starting ", file) -# CI <- rast(file) -# # names(CI) <- c("green","nir") -# message("raster loaded") -# -# CI <- CI[[2]]/CI[[1]]-1 -# # CI <- CI$nir/CI$green-1 -# message("CI calculated") -# CI <- terra::crop(CI, field_boundaries, mask = TRUE) #%>% CI_func() -# -# new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) -# writeRaster(CI, new_file, overwrite = TRUE) -# -# vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt")) -# terra::vrt(new_file, vrt_file, overwrite = TRUE) -# -# return(CI) -# } vrt_list <- list() -# for (file in raster_files) { -# v_crop <- create_mask_and_crop(file, field_boundaries) -# message(file, " processed") -# gc() -# } - for (file in filtered_files) { - v_crop <- create_mask_and_crop(file, field_boundaries) + v_crop <- create_mask_and_crop(file, field_boundaries, merged_final) emtpy_or_full <- global(v_crop, "notNA") vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt")) @@ -140,63 +106,59 @@ for (file in filtered_files) { gc() } -# Extracting CI - - -# pivot_sf_q <- st_read(here("..", "Data", "pivot_20210625.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% vect() -# pivot_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% group_by(pivot) %>% summarise() %>% vect() -# message("pivot loaded") raster_files_NEW <- list.files(merged_final,full.names = T, pattern = ".tif") -# pivots_dates0 <- readRDS(here(harvest_dir, "harvest_data_new")) %>% filter( -# pivot %in% c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", -# "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", -# "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3") -# ) - -# harvesting_data <- pivots_dates0 %>% -# select(c("pivot_quadrant", "season_start_2021", "season_end_2021", "season_start_2022", "season_end_2022", "season_start_2023", "season_end_2023", "season_start_2024", "season_end_2024")) %>% -# pivot_longer(cols = starts_with("season"), names_to = "Year", values_to = "value") %>% -# separate(Year, into = c("name", "Year"), sep = "(?<=season_start|season_end)\\_", remove = FALSE) %>% -# mutate(name = str_to_title(name)) %>% -# pivot_wider(names_from = name, values_from = value) %>% -# rename(field = pivot_quadrant) - -#If run for the firsttime, it will extract all data since the start and put into a table.rds. otherwise it will only add on to the existing table. - - # Define the path to the file -file_path <- here(cumulative_CI_vals_dir,"combined_CI_data.rds") +file_path <- here(cumulative_CI_vals_dir, "combined_CI_data.rds") # Check if the file exists if (!file.exists(file_path)) { - # Create the file with columns "column1" and "column2" - data <- data.frame(sub_field=NA, field=NA) - saveRDS(data, file_path) + # File does not exist, create it with all available data + print("combined_CI_data.rds does not exist. Preparing combined_CI_data.rds file for all available images.") + + # Extract data from all raster files + walk(raster_files_NEW, extract_rasters_daily, field_geojson = field_boundaries, quadrants = FALSE, daily_CI_vals_dir) + + # Combine all extracted data + extracted_values <- list.files(here(daily_CI_vals_dir), full.names = TRUE) + + pivot_stats <- extracted_values %>% + map(readRDS) %>% list_rbind() %>% + group_by(sub_field) + + # Save the combined data to the file + saveRDS(pivot_stats, file_path) + + print("All CI values extracted from all historic images and saved to combined_CI_data.rds.") + +} else { + # File exists, add new data + print("combined_CI_data.rds exists, adding the latest image data to the table.") + + # Filter and process the latest data + filtered_files <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>% + compact() %>% + flatten_chr() + + walk(filtered_files, extract_rasters_daily, field_geojson = field_boundaries, quadrants = TRUE, daily_CI_vals_dir) + + # Extract new values + extracted_values <- list.files(daily_CI_vals_dir, full.names = TRUE) + extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>% + compact() %>% + flatten_chr() + + pivot_stats <- extracted_values %>% + map(readRDS) %>% list_rbind() %>% + group_by(sub_field) + + # Load existing data and append new data + combined_CI_data <- readRDS(file_path) + pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data) + + # Save the updated combined data + saveRDS(pivot_stats2, file_path) + + print("All CI values extracted from the latest images and added to combined_CI_data.rds.") } - -print("combined_CI_data.rds exists, adding the latest image data to the table.") - -filtered_files <- map(dates$days_filter, ~ raster_files_NEW[grepl(pattern = .x, x = raster_files_NEW)]) %>% - compact() %>% - flatten_chr() -walk(filtered_files, extract_rasters_daily, field_geojson= field_boundaries, quadrants = TRUE, daily_CI_vals_dir) - -extracted_values <- list.files(daily_CI_vals_dir, full.names = TRUE) - -extracted_values <- map(dates$days_filter, ~ extracted_values[grepl(pattern = .x, x = extracted_values)]) %>% - compact() %>% - flatten_chr() - -pivot_stats <- extracted_values %>% - map(readRDS) %>% list_rbind() %>% - group_by(sub_field) - - -combined_CI_data <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #%>% drop_na(pivot_quadrant) -head(combined_CI_data) -pivot_stats2 <- bind_rows(pivot_stats, combined_CI_data) -# pivot_stats2 <- combined_CI_data -print("All CI values extracted from latest image.") -saveRDS(pivot_stats2, here(cumulative_CI_vals_dir,"combined_CI_data.rds")) #used to save the rest of the data into one file diff --git a/r_app/ci_extraction_utils.R b/r_app/ci_extraction_utils.R index 2bb8e52..f757480 100644 --- a/r_app/ci_extraction_utils.R +++ b/r_app/ci_extraction_utils.R @@ -1,6 +1,5 @@ # Utils for ci extraction - -date_list <- function(end_date, offset, week){ +date_list <- function(end_date, offset){ offset <- as.numeric(offset) - 1 end_date <- as.Date(end_date) start_date <- end_date - lubridate::days(offset) @@ -12,44 +11,26 @@ date_list <- function(end_date, offset, week){ return(list("week" = week, "year" = year, "days_filter" = days_filter)) } -# date_list <- function(weeks_in_the_paste){ -# week <- week(Sys.Date()- weeks(weeks_in_the_paste) ) -# year <- year(Sys.Date()- weeks(weeks_in_the_paste) ) -# days_filter <- Sys.Date() - weeks(weeks_in_the_paste) - days(0:6) -# -# return(c("week" = week, -# "year" = year, -# "days_filter" = list(days_filter))) -# -# } - -CI_func <- function(x, drop_layers = FALSE){ - CI <- x[[4]]/x[[2]]-1 - add(x) <- CI - names(x) <- c("red", "green", "blue","nir", "cloud" ,"CI") - if(drop_layers == FALSE){ - return(x) - }else{ - return(x$CI) - } -} - -mask_raster <- function(raster, fields){ - # x <- rast(filtered_files[1]) - x <- rast(raster) - emtpy_or_full <- global(x, sum) +create_mask_and_crop <- function(file, field_boundaries, merged_final_dir) { + message("starting ", file) + loaded_raster <- rast(file) + names(loaded_raster) <- c("Red", "Green", "Blue", "NIR") + message("raster loaded") - if(emtpy_or_full[1,] >= 2000000){ - names(x) <- c("red", "green", "blue","nir", "cloud") - cloud <- x$cloud - - cloud[cloud == 0 ] <- NA - x_masked <- mask(x, cloud, inverse = T) %>% crop(.,fields, mask = TRUE ) - x_masked <- x_masked %>% CI_func() - - message(raster, " masked") - return(x_masked) - } + CI <- loaded_raster$NIR / loaded_raster$Green - 1 + + loaded_raster$CI <- CI + message("CI calculated") + loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE) #%>% CI_func() + + loaded_raster[loaded_raster == 0] <- NA + new_file <- here(merged_final_dir, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) + writeRaster(loaded_raster, new_file, overwrite = TRUE) + + vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt")) + terra::vrt(new_file, vrt_file, overwrite = TRUE) + + return(loaded_raster) } date_extract <- function(file_path) { @@ -57,228 +38,16 @@ date_extract <- function(file_path) { } extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) { - # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE) - # date <- date_extract(filtered_files[1]) - # field_geojson <- sf::st_as_sf(pivot_sf_q) field_geojson <- sf::st_as_sf(field_geojson) x <- rast(file[1]) date <- date_extract(file) + pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>% st_drop_geometry() %>% rename("{date}" := mean_CI) + save_suffix <- if (quadrants){"quadrant"} else {"whole_field"} save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds")) + saveRDS(pivot_stats, save_path) } -right = function(text, num_char) { - substr(text, nchar(text) - (num_char-1), nchar(text)) -} - -extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) { - # field_names = "1.2A" - # harvesting_data = harvesting_data - # field_CI_data = pivot_stats_long - # season= 2021 - - filtered_harvesting_data <- harvesting_data %>% - filter(year == season, field %in% field_names) - - filtered_field_CI_data <- field_CI_data %>% - filter(field %in% field_names) - - # CI <- map_df(field_names, ~ { - ApproxFun <- approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value) - Dates <- seq.Date(ymd(min(filtered_field_CI_data$Date)), ymd(max(filtered_field_CI_data$Date)), by = 1) - LinearFit <- ApproxFun(Dates) - - CI <- data.frame(Date = Dates, FitData = LinearFit) %>% - left_join(., filtered_field_CI_data, by = "Date") %>% - filter(Date > filtered_harvesting_data$season_start & Date < filtered_harvesting_data$season_end) %>% - mutate(DOY = seq(1, n(), 1), - model = paste0("Data", season, " : ", field_names), - season = season, - field = field_names) - # }) #%>% - #{if (length(field_names) > 0) message("Done!")} - - return(CI) -} -# -# load_fields <- function(geojson_path) { -# field_geojson <- st_read(geojson_path) %>% -# select(pivot, pivot_quadrant) %>% -# vect() -# return(field_geojson) -# } -# -# load_harvest_data <- function(havest_data_path){ -# harvest_data <- readRDS(havest_data_path) -# return(harvest_data) -# } -# -# load_rasters <- function(raster_path, dates) { -# raster_files <- list.files(raster_path, full.names = TRUE, pattern = ".tif") -# filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>% -# compact() %>% -# flatten_chr() -# -# return(filtered_files) -# } -# -# mask_and_set_names <- function(filtered_files, fields) { -# rasters_masked <- map(filtered_files, mask_raster, fields = fields) %>% set_names(filtered_files) -# rasters_masked[sapply(rasters_masked, is.null)] <- NULL -# rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract)) -# -# return(rasters_masked) -# } -# -# calculate_total_pix_area <- function(filtered_files, fields_geojson) { -# # total_pix_area <- rast(filtered_files[1]) %>% -# # subset(1) %>% -# # crop(fields_geojson, mask = TRUE)%>% -# # global(.data, fun = "notNA") -# total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(fields_geojson, mask = TRUE) %>% freq(., usenames = TRUE) -# -# return(total_pix_area) -# } -# -# cloud_layer_extract <- function(rasters_masked){ -# cloud_layer_rast <- map(rasters_masked, function(spatraster) { -# spatraster[[5]] -# }) %>% rast() -# -# return(cloud_layer_rast) -# } -# -# calculate_cloud_coverage <- function(cloud_layer_rast, total_pix_area) { -# cloud_perc_list <- freq(cloud_layer_rast, usenames = TRUE) %>% -# mutate(cloud_perc = (100 -((count/total_pix_area$count)*100)), -# cloud_thres_5perc = as.integer(cloud_perc < 5), -# cloud_thres_40perc = as.integer(cloud_perc < 40)) %>% -# rename(Date = layer) %>% select(-value, -count) -# -# cloud_index_5perc <- which(cloud_perc_list$cloud_thres_5perc == max(cloud_perc_list$cloud_thres_5perc)) -# cloud_index_40perc <- which(cloud_perc_list$cloud_thres_40perc == max(cloud_perc_list$cloud_thres_40perc)) -# -# return(list(cloud_perc_list = cloud_perc_list, cloud_index_5perc = cloud_index_5perc, cloud_index_40perc = cloud_index_40perc)) -# } -# -# process_cloud_coverage <- function(cloud_coverage, rasters_masked) { -# if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) > 1) { -# message("More than 1 raster without clouds (<5%), max mosaic made ") -# -# cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_5perc] -# rsrc <- sprc(cloudy_rasters_list) -# x <- mosaic(rsrc) -# names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") -# -# } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) == 1) { -# message("Only 1 raster without clouds (<5%)") -# -# x <- rast(rasters_masked[cloud_coverage$cloud_index_5perc[1]]) -# names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") -# -# } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) > 1) { -# message("More than 1 image contains clouds, composite made of <40% cloud cover images") -# -# cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_40perc] -# rsrc <- sprc(cloudy_rasters_list) -# x <- mosaic(rsrc) -# names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") -# -# } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) == 0) { -# message("No cloud free images available") -# x <- rast(rasters_masked[1]) -# x[x] <- NA -# names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") -# } -# -# return(x) -# } - -# extract_rasters_daily_func <- function(daily_vals_dir, filtered_files, fields_geojson) { -# extracted_files <- walk(filtered_files, extract_rasters_daily, field_geojson = fields_geojson, quadrants = TRUE, daily_vals_dir) -# return(extracted_files) -# } - -# CI_load <- function(daily_vals_dir, grouping_variable){ -# extracted_values <- list.files(here(daily_vals_dir), full.names = TRUE) -# -# field_CI_values <- extracted_values %>% -# map_dfr(readRDS) %>% -# group_by(.data[[grouping_variable]]) %>% -# summarise(across(everything(), ~ first(na.omit(.)))) -# return(field_CI_values) -# } -# -# CI_long <- function(field_CI_values, pivot_long_cols){ -# field_CI_long <- field_CI_values %>% -# gather("Date", value, -pivot_long_cols) %>% -# mutate(Date = right(Date, 8), -# Date = ymd(Date) -# ) %>% -# drop_na(c("value","Date")) %>% -# mutate(value = as.numeric(value))%>% -# filter_all(all_vars(!is.infinite(.)))%>% -# rename(field = pivot_quadrant) -# -# return(field_CI_long) -# } -# -# process_year_data <- function(year, harvest_data, field_CI_long) { -# pivots_dates_year <- harvest_data %>% na.omit() %>% filter(year == year) -# pivot_select_model_year <- unique(pivots_dates_year$field) -# -# data <- map_dfr(pivot_select_model_year, ~ extract_CI_data(.x, harvest_data, field_CI_long, season = year)) -# -# return(data) -# } - - - -#functions for CI_data_prep -# create_mask_and_crop <- function(file, pivot_sf_q) { -# # file <- filtered_files[5] -# message("starting ", file) -# loaded_raster <- rast(file) -# names(loaded_raster) <- c("Red", "Green", "Blue", "NIR") -# # names(CI) <- c("green","nir") -# message("raster loaded") -# -# # CI <- CI[[2]] / CI[[1]] - 1 -# CI <- loaded_raster$NIR / loaded_raster$Green - 1 -# -# loaded_raster$CI <- CI -# # CI <- CI$nir/CI$green-1 -# message("CI calculated") -# loaded_raster <- terra::crop(loaded_raster, pivot_sf_q, mask = TRUE) #%>% CI_func() -# -# loaded_raster[loaded_raster == 0] <- NA -# # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI") -# # v_crop$CI <- v_crop$CI - 1 -# new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) -# writeRaster(loaded_raster, new_file, overwrite = TRUE) -# -# vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt")) -# terra::vrt(new_file, vrt_file, overwrite = TRUE) -# -# # v_crop <- mask_raster(v, pivot_sf_q) -# return(loaded_raster) -# } - - -# extract_rasters_daily <- function(file, field_geojson, quadrants = TRUE, save_dir) { -# # x <- rast(filtered_files[1])%>% CI_func(drop_layers = TRUE) -# # date <- date_extract(filtered_files[1]) -# # field_geojson <- sf::st_as_sf(pivot_sf_q) -# field_geojson <- sf::st_as_sf(field_geojson) -# x <- rast(file[1]) -# date <- date_extract(file) -# pivot_stats <- cbind(field_geojson, mean_CI = round(exactextractr::exact_extract(x$CI, field_geojson, fun = "mean"), 2)) %>% -# st_drop_geometry() %>% rename("{date}" := mean_CI) -# save_suffix <- if (quadrants){"quadrant"} else {"whole_field"} -# save_path <- here(save_dir, paste0("extracted_", date, "_", save_suffix, ".rds")) -# saveRDS(pivot_stats, save_path) -# } - diff --git a/r_app/interpolate_growth_model.R b/r_app/interpolate_growth_model.R index 0c046e9..b5f43eb 100644 --- a/r_app/interpolate_growth_model.R +++ b/r_app/interpolate_growth_model.R @@ -5,7 +5,7 @@ library(terra) library(tidyverse) library(lubridate) library(exactextractr) -library(readxl) +# library(readxl) # Vang alle command line argumenten op @@ -37,10 +37,13 @@ weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic") daily_vrt <- here(data_dir, "vrt") harvest_dir <- here(data_dir, "HarvestData") +source(here("r_app/parameters_project.R")) +source(here("r_app/ci_extraction_utils.R")) +# source(here("r_app/mosaic_creation_utils.R")) + source("parameters_project.R") source("ci_extraction_utils.R") - pivot_stats2 <- readRDS(here(cumulative_CI_vals_dir,"combined_CI_data.rds")) %>% ungroup() %>% group_by(field, sub_field) %>% @@ -61,76 +64,74 @@ pivot_stats_long <- pivot_stats2 %>% # sub_field = field) %>% unique() - -# #2021 -# pivot_select_model_Data_2021 <- harvesting_data %>% filter(Year == 2021) %>% pull(field) -# -# pivot_select_model_Data_2022 <- harvesting_data %>% filter(Year == 2022) %>% pull(field) - -# pivot_select_model_Data_2023 <- harvesting_data %>% filter(year == 2023) %>% filter(!is.na(season_start)) %>% pull(sub_field) - -pivot_select_model_Data_2024 <- harvesting_data %>% filter(year == 2024)%>% filter(!is.na(season_start)) %>% pull(sub_field) -message(pivot_select_model_Data_2024) -# pivots_dates_Data_2022 <- pivots_dates0 %>% filter(!is.na(season_end_2022)) -# pivot_select_model_Data_2022 <- unique(pivots_dates_Data_2022$pivot_quadrant ) -# -# pivots_dates_Data_2023 <- pivots_dates0 %>% filter(!is.na(season_start_2023)) -# pivot_select_model_Data_2023 <- unique(pivots_dates_Data_2023$pivot_quadrant) -# -# pivots_dates_Data_2024 <- pivots_dates0 %>% filter(!is.na(season_start_2024)) -# pivot_select_model_Data_2024 <- unique(pivots_dates_Data_2024$pivot_quadrant) +years <- harvesting_data %>% + filter(!is.na(season_start)) %>% + distinct(year) %>% + pull(year) extract_CI_data <- function(field_names, harvesting_data, field_CI_data, season) { - # field_names = "Nandi A1a" - # field_names = "1.13A" - # harvesting_data = harvesting_data - # field_CI_data = pivot_stats_long - # season= 2024 - + # Filter harvesting data for the given season and field names filtered_harvesting_data <- harvesting_data %>% - # na.omit() %>% filter(year == season, sub_field %in% field_names) - + + # Filter field CI data for the given field names filtered_field_CI_data <- field_CI_data %>% filter(sub_field %in% field_names) - if (nrow(filtered_field_CI_data) == 0) { - return(data.frame()) # Return an empty data frame if no data is found - } - - # CI <- map_df(field_names, ~ { + + # Return an empty data frame if no CI data is found + if (nrow(filtered_field_CI_data) == 0) { + return(data.frame()) + } + + # Create a linear interpolation function for the CI data ApproxFun <- approxfun(x = filtered_field_CI_data$Date, y = filtered_field_CI_data$value) Dates <- seq.Date(ymd(min(filtered_field_CI_data$Date)), ymd(max(filtered_field_CI_data$Date)), by = 1) LinearFit <- ApproxFun(Dates) - - CI <- data.frame(Date = Dates, FitData = LinearFit) %>% + + # Combine interpolated data with the original CI data + CI <- data.frame(Date = Dates, FitData = LinearFit) %>% left_join(., filtered_field_CI_data, by = "Date") %>% - filter(Date > filtered_harvesting_data$season_start & Date < filtered_harvesting_data$season_end) %>% + filter(Date > filtered_harvesting_data$season_start & Date < filtered_harvesting_data$season_end) + + # If CI is empty after filtering, return an empty dataframe + if (nrow(CI) == 0) { + return(data.frame()) + } + + # Add additional columns if data exists + CI <- CI %>% mutate(DOY = seq(1, n(), 1), model = paste0("Data", season, " : ", field_names), season = season, - sub_field = field_names) - # }) #%>% - + subField = field_names) + return(CI) } -## Extracting the correct CI values -# Data_2021 <- map(pivot_select_model_Data_2021, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2021)) %>% list_rbind() -# message('2021') -# Data_2022 <- map(pivot_select_model_Data_2022, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2022)) %>% list_rbind() -# message('2022') -# Data_2023 <- map(pivot_select_model_Data_2023, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2023)) %>% list_rbind() -# message('2023') -head(harvesting_data) -Data_2024 <- map(pivot_select_model_Data_2024, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = 2024)) %>% list_rbind() -message('2024') -head(Data_2024) -#CI_all <- rbind(Data_2023, Data_2024) -CI_all <- Data_2024 -message('CI_all created') -#CI_all <- Data_2023 +CI_all <- map_df(years, function(yr) { + # yr = 2021 + message(yr) + + # Get the fields harvested in this year + sub_fields <- harvesting_data %>% + filter(year == yr) %>% + filter(!is.na(season_start)) %>% + pull(sub_field) + + # Filter sub_fields to only include those with value data in pivot_stats_long + valid_sub_fields <- sub_fields %>% + keep(~ any(pivot_stats_long$sub_field == .x)) + + # Extract data for each valid field + map(valid_sub_fields, ~ extract_CI_data(.x, harvesting_data = harvesting_data, field_CI_data = pivot_stats_long, season = yr)) %>% + list_rbind() +}) + +CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitData), + cumulative_CI = cumsum(FitData)) + CI_all <- CI_all %>% group_by(model) %>% mutate(CI_per_day = FitData - lag(FitData), cumulative_CI = cumsum(FitData)) diff --git a/r_app/mosaic_creation.R b/r_app/mosaic_creation.R index 316ecf7..53de9e2 100644 --- a/r_app/mosaic_creation.R +++ b/r_app/mosaic_creation.R @@ -11,7 +11,7 @@ library(terra) library(tidyverse) library(lubridate) # library(exactextractr) -library(readxl) +# library(readxl) #funcion CI_prep package @@ -62,21 +62,24 @@ weekly_CI_mosaic <- here(laravel_storage_dir, "weekly_mosaic") daily_vrt <- here(data_dir, "vrt") harvest_dir <- here(data_dir, "HarvestData") +source(here("r_app/parameters_project.R")) +source(here("r_app/mosaic_creation_utils.R")) + source("parameters_project.R") source("mosaic_creation_utils.R") -dir.create(here(laravel_storage_dir)) -dir.create(here(reports_dir)) -dir.create(here(data_dir)) -dir.create(here(extracted_CI_dir)) -dir.create(here(daily_CI_vals_dir)) -dir.create(here(cumulative_CI_vals_dir)) -dir.create(here(weekly_CI_mosaic)) -dir.create(here(daily_vrt)) -dir.create(merged_final) -dir.create(harvest_dir) +# dir.create(here(laravel_storage_dir)) +# dir.create(here(reports_dir)) +# dir.create(here(data_dir)) +# dir.create(here(extracted_CI_dir)) +# dir.create(here(daily_CI_vals_dir)) +# dir.create(here(cumulative_CI_vals_dir)) +# dir.create(here(weekly_CI_mosaic)) +# dir.create(here(daily_vrt)) +# dir.create(merged_final) +# dir.create(harvest_dir) -# end_date <- lubridate::dmy("20-6-2024") +# end_date <- lubridate::dmy("28-8-2024") week <- week(end_date) @@ -101,78 +104,81 @@ vrt_list <- map(dates$days_filter, ~ vrt_files[grepl(pattern = .x, x = vrt_file compact() %>% flatten_chr() - raster_files_final <- list.files(merged_final,full.names = T, pattern = ".tif") -if (length(vrt_list) > 0 ){ +raster_files_final <- list.files(merged_final,full.names = T, pattern = ".tif") + +if (length(vrt_list) > 0) { print("vrt list made, preparing mosaic creation by counting cloud cover") - - total_pix_area <- rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>% + + total_pix_area <- + rast(vrt_list[1]) %>% terra::subset(1) %>% setValues(1) %>% crop(field_boundaries, mask = TRUE) %>% - global(., fun="notNA") #%>% - + global(., fun = "notNA") #%>% + layer_5_list <- purrr::map(vrt_list, function(vrt_list) { rast(vrt_list[1]) %>% terra::subset(1) }) %>% rast() - - missing_pixels_count <- layer_5_list %>% global(., fun="notNA") %>% + + missing_pixels_count <- + layer_5_list %>% global(., fun = "notNA") %>% mutate( total_pixels = total_pix_area$notNA, - missing_pixels_percentage = round(100 -((notNA/total_pix_area$notNA)*100)), + missing_pixels_percentage = round(100 - (( + notNA / total_pix_area$notNA + ) * 100)), thres_5perc = as.integer(missing_pixels_percentage < 5), thres_40perc = as.integer(missing_pixels_percentage < 45) ) - message('hier') - message(missing_pixels_count) - index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc) ) + + index_5perc <- which(missing_pixels_count$thres_5perc == max(missing_pixels_count$thres_5perc)) index_40perc <- which(missing_pixels_count$thres_40perc == max(missing_pixels_count$thres_40perc)) - + ## Create mosaic - - if(sum(missing_pixels_count$thres_5perc)>1){ + + if (sum(missing_pixels_count$thres_5perc) > 1) { message("More than 1 raster without clouds (<5%), max composite made") - + cloudy_rasters_list <- vrt_list[index_5perc] - + rsrc <- sprc(cloudy_rasters_list) x <- terra::mosaic(rsrc, fun = "max") names(x) <- c("Red", "Green", "Blue", "NIR", "CI") - - }else if(sum(missing_pixels_count$thres_5perc)==1){ + + } else if (sum(missing_pixels_count$thres_5perc) == 1) { message("Only 1 raster without clouds (<5%)") - + x <- rast(vrt_list[index_5perc[1]]) names(x) <- c("Red", "Green", "Blue", "NIR", "CI") - - }else if(sum(missing_pixels_count$thres_40perc)>1){ + + } else if (sum(missing_pixels_count$thres_40perc) > 1) { message("More than 1 image contains clouds, composite made of <40% cloud cover images") - + cloudy_rasters_list <- vrt_list[index_40perc] - + rsrc <- sprc(cloudy_rasters_list) x <- mosaic(rsrc, fun = "max") names(x) <- c("Red", "Green", "Blue", "NIR", "CI") - - }else if(sum(missing_pixels_count$thres_40perc)==1){ + + } else if (sum(missing_pixels_count$thres_40perc) == 1) { message("Only 1 image available but contains clouds") - + x <- rast(vrt_list[index_40perc[1]]) names(x) <- c("Red", "Green", "Blue", "NIR", "CI") - - }else if(sum(missing_pixels_count$thres_40perc)==0){ + + } else if (sum(missing_pixels_count$thres_40perc) == 0) { message("No cloud free images available, all images combined") - + rsrc <- sprc(vrt_list) x <- mosaic(rsrc, fun = "max") names(x) <- c("Red", "Green", "Blue", "NIR", "CI") - + } - + } else{ - message("No images available this week, empty mosaic created") - + x <- rast(raster_files_final[1]) %>% setValues(0) %>% crop(field_boundaries, mask = TRUE) - + names(x) <- c("Red", "Green", "Blue", "NIR", "CI") } plot(x$CI, main = paste("CI map ", dates$week)) diff --git a/r_app/mosaic_creation_utils.R b/r_app/mosaic_creation_utils.R index b07176f..68aadb8 100644 --- a/r_app/mosaic_creation_utils.R +++ b/r_app/mosaic_creation_utils.R @@ -11,207 +11,3 @@ date_list <- function(end_date, offset){ return(list("week" = week, "year" = year, "days_filter" = days_filter)) } - -# date_list <- function(weeks_in_the_paste){ -# week <- week(Sys.Date()- weeks(weeks_in_the_paste) ) -# year <- year(Sys.Date()- weeks(weeks_in_the_paste) ) -# days_filter <- Sys.Date() - weeks(weeks_in_the_paste) - days(0:6) -# -# return(c("week" = week, -# "year" = year, -# "days_filter" = list(days_filter))) -# -# } -create_mask_and_crop <- function(file, field_boundaries) { - # file <- filtered_files[5] - message("starting ", file) - loaded_raster <- rast(file) - names(loaded_raster) <- c("Red", "Green", "Blue", "NIR") - # names(CI) <- c("green","nir") - message("raster loaded") - - # CI <- CI[[2]] / CI[[1]] - 1 - CI <- loaded_raster$NIR / loaded_raster$Green - 1 - - loaded_raster$CI <- CI - # CI <- CI$nir/CI$green-1 - message("CI calculated") - loaded_raster <- terra::crop(loaded_raster, field_boundaries, mask = TRUE) #%>% CI_func() - - loaded_raster[loaded_raster == 0] <- NA - # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI") - # v_crop$CI <- v_crop$CI - 1 - new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) - writeRaster(loaded_raster, new_file, overwrite = TRUE) - - vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt")) - terra::vrt(new_file, vrt_file, overwrite = TRUE) - - # v_crop <- mask_raster(v, pivot_sf_q) - return(loaded_raster) -} - -CI_func <- function(x, drop_layers = FALSE){ - CI <- x[[4]]/x[[2]]-1 - add(x) <- CI - names(x) <- c("red", "green", "blue","nir", "cloud" ,"CI") - if(drop_layers == FALSE){ - return(x) - }else{ - return(x$CI) - } -} - -mask_raster <- function(raster, fields){ - # x <- rast(filtered_files[1]) - x <- rast(raster) - emtpy_or_full <- global(x, sum) - - if(emtpy_or_full[1,] >= 2000000){ - names(x) <- c("red", "green", "blue","nir", "cloud") - cloud <- x$cloud - - cloud[cloud == 0 ] <- NA - x_masked <- mask(x, cloud, inverse = T) %>% crop(.,fields, mask = TRUE ) - x_masked <- x_masked %>% CI_func() - - message(raster, " masked") - return(x_masked) - } -} - -date_extract <- function(file_path) { - str_extract(file_path, "\\d{4}-\\d{2}-\\d{2}") -} - - - -right = function(text, num_char) { - substr(text, nchar(text) - (num_char-1), nchar(text)) -} - -# load_fields <- function(geojson_path) { -# field_geojson <- st_read(geojson_path) %>% -# select(pivot, pivot_quadrant) %>% -# vect() -# return(field_geojson) -# } -# -# load_harvest_data <- function(havest_data_path){ -# harvest_data <- readRDS(havest_data_path) -# return(harvest_data) -# } - -# load_rasters <- function(raster_path, dates) { -# raster_files <- list.files(raster_path, full.names = TRUE, pattern = ".tif") -# filtered_files <- map(dates$days_filter, ~ raster_files[grepl(pattern = .x, x = raster_files)]) %>% -# compact() %>% -# flatten_chr() -# -# return(filtered_files) -# } - -# mask_and_set_names <- function(filtered_files, fields) { -# rasters_masked <- map(filtered_files, mask_raster, fields = fields) %>% set_names(filtered_files) -# rasters_masked[sapply(rasters_masked, is.null)] <- NULL -# rasters_masked <- setNames(rasters_masked, map_chr(names(rasters_masked), date_extract)) -# -# return(rasters_masked) -# } - -# calculate_total_pix_area <- function(filtered_files, fields_geojson) { -# # total_pix_area <- rast(filtered_files[1]) %>% -# # subset(1) %>% -# # crop(fields_geojson, mask = TRUE)%>% -# # global(.data, fun = "notNA") -# total_pix_area <- rast(filtered_files[1]) %>% subset(1) %>% crop(fields_geojson, mask = TRUE) %>% freq(., usenames = TRUE) -# -# return(total_pix_area) -# } - -# cloud_layer_extract <- function(rasters_masked){ -# cloud_layer_rast <- map(rasters_masked, function(spatraster) { -# spatraster[[5]] -# }) %>% rast() -# -# return(cloud_layer_rast) -# } - -# calculate_cloud_coverage <- function(cloud_layer_rast, total_pix_area) { -# cloud_perc_list <- freq(cloud_layer_rast, usenames = TRUE) %>% -# mutate(cloud_perc = (100 -((count/total_pix_area$count)*100)), -# cloud_thres_5perc = as.integer(cloud_perc < 5), -# cloud_thres_40perc = as.integer(cloud_perc < 40)) %>% -# rename(Date = layer) %>% select(-value, -count) -# -# cloud_index_5perc <- which(cloud_perc_list$cloud_thres_5perc == max(cloud_perc_list$cloud_thres_5perc)) -# cloud_index_40perc <- which(cloud_perc_list$cloud_thres_40perc == max(cloud_perc_list$cloud_thres_40perc)) -# -# return(list(cloud_perc_list = cloud_perc_list, cloud_index_5perc = cloud_index_5perc, cloud_index_40perc = cloud_index_40perc)) -# } - -process_cloud_coverage <- function(cloud_coverage, rasters_masked) { - if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) > 1) { - message("More than 1 raster without clouds (<5%), max mosaic made ") - - cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_5perc] - rsrc <- sprc(cloudy_rasters_list) - x <- mosaic(rsrc) - names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") - - } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_5perc) == 1) { - message("Only 1 raster without clouds (<5%)") - - x <- rast(rasters_masked[cloud_coverage$cloud_index_5perc[1]]) - names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") - - } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) > 1) { - message("More than 1 image contains clouds, composite made of <40% cloud cover images") - - cloudy_rasters_list <- rasters_masked[cloud_coverage$cloud_index_40perc] - rsrc <- sprc(cloudy_rasters_list) - x <- mosaic(rsrc) - names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") - - } else if (sum(cloud_coverage$cloud_perc_list$cloud_thres_40perc) == 0) { - message("No cloud free images available") - x <- rast(rasters_masked[1]) - x[x] <- NA - names(x) <- c("red", "green", "blue", "nir", "cloud", "CI") - } - - return(x) -} - - -#functions for CI_data_prep -# create_mask_and_crop <- function(file, pivot_sf_q) { -# # file <- filtered_files[5] -# message("starting ", file) -# loaded_raster <- rast(file) -# names(loaded_raster) <- c("Red", "Green", "Blue", "NIR") -# # names(CI) <- c("green","nir") -# message("raster loaded") -# -# # CI <- CI[[2]] / CI[[1]] - 1 -# CI <- loaded_raster$NIR / loaded_raster$Green - 1 -# -# loaded_raster$CI <- CI -# # CI <- CI$nir/CI$green-1 -# message("CI calculated") -# loaded_raster <- terra::crop(loaded_raster, pivot_sf_q, mask = TRUE) #%>% CI_func() -# -# loaded_raster[loaded_raster == 0] <- NA -# # names(v_crop) <- c("red", "green", "blue","nir", "cloud" ,"CI") -# # v_crop$CI <- v_crop$CI - 1 -# new_file <- here(merged_final, paste0(tools::file_path_sans_ext(basename(file)), ".tif")) -# writeRaster(loaded_raster, new_file, overwrite = TRUE) -# -# vrt_file <- here(daily_vrt, paste0(tools::file_path_sans_ext(basename(file)), ".vrt")) -# terra::vrt(new_file, vrt_file, overwrite = TRUE) -# -# # v_crop <- mask_raster(v, pivot_sf_q) -# return(loaded_raster) -# } - - diff --git a/r_app/parameters_project.R b/r_app/parameters_project.R index ebb8de9..78d00a1 100644 --- a/r_app/parameters_project.R +++ b/r_app/parameters_project.R @@ -1,111 +1,38 @@ library('readxl') #chemba -if(project_dir == "chemba1"){ - message("Yield data for Chemba") - - field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant)%>% - mutate(sub_area = case_when(pivot %in% c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", "1.14" , "1.16" , "1.17" , "1.18" , "6.1", "6.2", "DL1.1", "DL1.3") ~ "estate", - TRUE ~ "Cooperative")) - names(field_boundaries_sf) <- c("field", "sub_field", "geometry", "sub_area") - - field_boundaries <- field_boundaries_sf %>% vect() - names(field_boundaries) <- c("field", "sub_field", "sub_area") - - # field_boundaries_merged <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(pivot, pivot_quadrant) %>% group_by(pivot) %>% summarise() %>% vect() - - joined_spans <-st_read(here(data_dir, "span.geojson")) %>% st_transform(crs(field_boundaries_sf)) - names(joined_spans) <- c("field", "area", "radius", "spans", "span", "geometry") - - - pivots_dates0 <- readRDS(here(harvest_dir, "harvest_data_new")) %>% filter( - pivot %in% c("1.1", "1.2", "1.3", "1.4", "1.6", "1.7", "1.8", "1.9", "1.10", "1.11", "1.12", "1.13", - "1.14" , "1.16" , "1.17" , "1.18" ,"2.1", "2.2", "2.3" , "2.4", "2.5", "3.1", "3.2", "3.3", - "4.1", "4.2", "4.3", "4.4", "4.5", "4.6", "5.1" ,"5.2", "5.3", "5.4", "6.1", "6.2", "DL1.1", "DL1.3") +field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson")) + +names(field_boundaries_sf) <- c("field", "sub_field", "geometry") +field_boundaries <- field_boundaries_sf %>% terra::vect() + +harvesting_data <- read_excel(here(data_dir, "harvest.xlsx")) %>% + dplyr::select( + c( + "field", + "sub_field", + "year", + "season_start", + "season_end", + "age", + "sub_area", + "tonnage_ha" + ) + ) %>% + mutate( + field = as.character(field), + sub_field = as.character(sub_field), + year = as.numeric(year), + season_start = as.Date(season_start), + season_end = as.Date(season_end), + age = as.numeric(age), + sub_area = as.character(sub_area), + tonnage_ha = as.numeric(tonnage_ha) + ) %>% + mutate( + season_end = case_when(season_end > Sys.Date() ~ Sys.Date(), + TRUE ~ season_end), + age = round(as.numeric(season_end - season_start) / 7, 0) ) - - harvesting_data <- pivots_dates0 %>% - dplyr::select(c("pivot_quadrant", "pivot", "season_start_2021", "season_end_2021", "season_start_2022", - "season_end_2022", "season_start_2023", "season_end_2023", "season_start_2024", "season_end_2024", "Age")) %>% - pivot_longer(cols = starts_with("season"), names_to = "Year", values_to = "value") %>% - separate(Year, into = c("name", "Year"), sep = "(?<=season_start|season_end)\\_", remove = FALSE) %>% - mutate(name = str_to_title(name)) %>% - pivot_wider(names_from = name, values_from = value) %>% - rename(Field = pivot, - subField = pivot_quadrant) - - - - - - - - -} else if (project_dir == "xinavane"){ - library(readxl) - message("Yield data for Xinavane") - field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson")) %>% dplyr::select(-Pivot) - names(field_boundaries_sf) <- c("field", "sub_field", "geometry") - - field_boundaries <- field_boundaries_sf %>% vect() - - joined_spans <- field_boundaries_sf %>% dplyr::select(Field) - - pivots_dates0 <- read_excel(here(harvest_dir, "Yield data.xlsx"), - skip = 3, - col_types = c("numeric", "text", "skip", "text", "numeric", "numeric", "numeric", "numeric", "date", - "numeric", "skip", "numeric")) %>% - rename( - year = 1, - field = 2, - sub_area = 3, - hectares = 4, - tons = 5, - tcha = 6, - tchy = 7, - season_end = 8, - age = 9, - ratoon = 10 - ) %>% - mutate(Season_end = ymd(season_end), - Season_start = as.Date(season_end - (duration(months = age))), - subField = Field) #don't forget to add 2022 as a year for the 'curent' season - - pivots_dates0 <- pivots_dates0 %>% - mutate(year = year + 6) - - # Add 6 years to Season_end column - pivots_dates0 <- pivots_dates0 %>% - mutate(season_end = season_end + years(6)) - - # Add 6 years to Season_start column - pivots_dates0 <- pivots_dates0 %>% - mutate(season_start = season_start + years(6)) - - - harvesting_data <- pivots_dates0 %>% dplyr::select(c("field","sub_field", "year", "season_start","season_end", "age" , "sub_area", "tcha")) - - -} else { - field_boundaries_sf <- st_read(here(data_dir, "pivot.geojson")) - head(field_boundaries_sf) - names(field_boundaries_sf) <- c("field", "sub_field", "geometry") - field_boundaries <- field_boundaries_sf %>% terra::vect() - harvesting_data <- read_excel(here(data_dir, "harvest.xlsx")) %>% - dplyr::select(c("field", "sub_field", "year", "season_start", "season_end", "age", "sub_area", "tonnage_ha")) %>% - mutate(field = as.character(field), - sub_field = as.character(sub_field), - year = as.numeric(year), - season_start = as.Date(season_start), - season_end = as.Date(season_end), - age = as.numeric(age), - sub_area = as.character(sub_area), - tonnage_ha = as.numeric(tonnage_ha) - ) %>% - mutate(season_end = case_when( - season_end > Sys.Date() ~ Sys.Date(), - TRUE ~ season_end), - age = round(as.numeric(season_end - season_start)/7,0)) -}