From f31ff4579d907e9526c6ea041aea9ee30ca6628d Mon Sep 17 00:00:00 2001 From: Beni Stocker Date: Thu, 12 Sep 2024 22:28:04 +0200 Subject: [PATCH] filled gaps in canopy and reference height --- data-raw/02_compile_final_site_list.R | 183 ++++++++++++++++++++++++-- data/fdk_site_info.rda | Bin 9196 -> 10240 bytes 2 files changed, 173 insertions(+), 10 deletions(-) diff --git a/data-raw/02_compile_final_site_list.R b/data-raw/02_compile_final_site_list.R index fb0dafe..6b1fa61 100644 --- a/data-raw/02_compile_final_site_list.R +++ b/data-raw/02_compile_final_site_list.R @@ -349,7 +349,8 @@ df <- df |> ## root zone water storage capacity--------------------------------------------- # using the map from Stocker et al., 2023, obtainable from Zenodo at https://doi.org/10.5281/zenodo.5515246 -whc <- raster("/data/archive/whc_stocker_2023/data/zroot_cwdx80_forcing.nc") +# whc <- raster("/data/archive/whc_stocker_2023/data/zroot_cwdx80_forcing.nc") # on geco server +whc <- raster("~/data/mct_data/zroot_cwdx80_forcing.nc") whc_v <- raster::extract(whc, loc) # append to original data frame @@ -369,7 +370,8 @@ df <- df |> ## Get still missing elevation data from ETOPO1--------------------------------- # file is too large to add it to this repo. -etopo <- raster("/data/archive/etopo_NA_NA/data/ETOPO1_Bed_g_geotiff.tif") +# etopo <- raster("/data/archive/etopo_NA_NA/data/ETOPO1_Bed_g_geotiff.tif") # on geco server +etopo <- raster("~/data/etopo/ETOPO1_Bed_g_geotiff.tif") # on geco server etopo_v <- raster::extract(etopo, loc) # add etopo1 column @@ -408,7 +410,7 @@ df <- df |> # complement with AmeriFlux meta info left_join( ameriflux |> - select( + dplyr::select( sitename, canopy_height_amf = canopy_height, reference_height_amf = reference_height @@ -416,14 +418,14 @@ df <- df |> join_by(sitename) ) |> mutate( - ifelse(is.na(canopy_height), canopy_height_amf, canopy_height), - ifelse(is.na(reference_height), reference_height_amf, reference_height) + canopy_height = ifelse(is.na(canopy_height), canopy_height_amf, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_amf, reference_height) ) |> # complement with FLUXNET2015 meta info left_join( siteinfo_fluxnet2015 |> - select( + dplyr::select( sitename = id, canopy_height_flx = canopy_height, reference_height_flx = reference_height @@ -431,10 +433,173 @@ df <- df |> join_by(sitename) ) |> mutate( - ifelse(is.na(canopy_height), canopy_height_flx, canopy_height), - ifelse(is.na(reference_height), reference_height_flx, reference_height) + canopy_height= ifelse(is.na(canopy_height), canopy_height_flx, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_flx, reference_height) + ) |> + + # complement with Plumber data + left_join( + plumber |> + dplyr::select( + sitename, + canopy_height_plb = canopy_height, + reference_height_plb = reference_height + ), + join_by(sitename) + ) |> + mutate( + canopy_height= ifelse(is.na(canopy_height), canopy_height_plb, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_plb, reference_height) ) +# get mean ratio of measurement height over vegetation height (for later) +vec_ratio <- df |> + mutate(ratio = reference_height/canopy_height) |> + ungroup() |> + pull(ratio) + +ratio_q10 <- quantile(vec_ratio, probs = 0.25, na.rm = TRUE) +ratio_q50 <- quantile(vec_ratio, probs = 0.5, na.rm = TRUE) +ratio_q90 <- quantile(vec_ratio, probs = 0.75, na.rm = TRUE) + +# fill by mean by vegetation type and major climate zone (biome) +df_fill <- df |> + mutate(biome = stringr::str_sub(koeppen_code, start = 1, end = 1)) |> + group_by(biome, classid) |> + summarise( + canopy_height = mean(canopy_height, na.rm = TRUE), + reference_height = mean(reference_height, na.rm = TRUE) + ) + +# fill remaining only by veg type +df_fill2 <- df |> + group_by(classid) |> + summarise( + canopy_height2 = mean(canopy_height, na.rm = TRUE), + reference_height2 = mean(reference_height, na.rm = TRUE) + ) + +# fill remaining only by climate zone +df_fill3 <- df |> + group_by(koeppen_code) |> + summarise( + canopy_height3 = mean(canopy_height, na.rm = TRUE), + reference_height3 = mean(reference_height, na.rm = TRUE) + ) + +# df_fill <- df_fill |> +# left_join( +# df_fill2, +# join_by(classid) +# ) |> +# mutate( +# canopy_height = ifelse(is.na(canopy_height), canopy_height2, canopy_height), +# reference_height = ifelse(is.na(reference_height), reference_height2, reference_height) +# ) |> +# dplyr::select( +# -reference_height2, +# -canopy_height2 +# ) + +# some are still missing: take ENF for DNF +df_fill4 <- df_fill |> + filter( + classid == "ENF" + ) |> + rename( + canopy_height4 = canopy_height, + reference_height4 = reference_height + ) + +# fill missing +df <- df |> + + # fill with average by biome and vegetation type + mutate(biome = stringr::str_sub(koeppen_code, start = 1, end = 1)) |> + left_join( + df_fill |> + rename( + canopy_height_fill = canopy_height, + reference_height_fill = reference_height + ), + join_by(biome, classid) + ) |> + mutate( + canopy_height = ifelse(is.na(canopy_height), canopy_height_fill, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_fill, reference_height) + ) |> + dplyr::select( + -reference_height_fill, + -canopy_height_fill + ) |> + + # fill with average by vegetation type + left_join( + df_fill2 |> + rename( + canopy_height_fill = canopy_height2, + reference_height_fill = reference_height2 + ), + join_by(classid) + ) |> + mutate( + canopy_height = ifelse(is.na(canopy_height), canopy_height_fill, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_fill, reference_height) + ) |> + dplyr::select( + -reference_height_fill, + -canopy_height_fill + ) |> + + # fill remaining with ENF for DNF + left_join( + df_fill4 |> + mutate(classid = "DNF") |> + rename( + canopy_height_fill = canopy_height4, + reference_height_fill = reference_height4 + ), + join_by(biome, classid) + ) |> + mutate( + canopy_height = ifelse(is.na(canopy_height), canopy_height_fill, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_fill, reference_height) + ) |> + dplyr::select( + -reference_height_fill, + -canopy_height_fill + ) |> + + # fill with average by koeppen climate + left_join( + df_fill3 |> + rename( + canopy_height_fill = canopy_height3, + reference_height_fill = reference_height3 + ), + join_by(koeppen_code) + ) |> + mutate( + canopy_height = ifelse(is.na(canopy_height), canopy_height_fill, canopy_height), + reference_height = ifelse(is.na(reference_height), reference_height_fill, reference_height) + ) |> + dplyr::select( + -reference_height_fill, + -canopy_height_fill + ) + +# correct those where ratio of canopy and vegetation height is in 10% quantiles +df <- df |> + mutate(ratio = reference_height/canopy_height) |> + mutate(reference_height = ifelse(ratio > ratio_q90, canopy_height * ratio_q50, reference_height)) |> + mutate(reference_height = ifelse(ratio < ratio_q10, canopy_height * ratio_q50, reference_height)) + + +# df |> +# mutate(ratio = reference_height/canopy_height) |> +# ggplot(aes(ratio, ..count..)) + +# geom_histogram() + # ## get IGBP from MODIS data----------------------------------------------------- # # download IGBP class values (MODISTools) # igbp <- apply(df, 1, function(x){ @@ -509,8 +674,6 @@ message("IGBP: found ", print(igbp_by_source |> dplyr::filter(IGBP_fdk != IGBP_eu | IGBP_fdk != IGBP_flx2015)) - - # save the data, retaining only key columns (note: valuable information also in # other columns!) fdk_site_info <- df |> diff --git a/data/fdk_site_info.rda b/data/fdk_site_info.rda index d7b1c6bbe7e1c10856b5f491de8f4888f7fd7a1e..40546f06fea9b7867182c42046d2e654f2d392be 100644 GIT binary patch delta 4206 zcmV-!5RvcfM}Sa}83eiBC&IBEr5Aro84%EUlyvpj8&VXX>`X?i|L`C2{I7?(H9fZ9 z!Nu%S{%elyJp6}cZZMC1bJA*rpb$3@K%|p0+H^~vsd4#N_qKqqCa(~C3Ph@6N%vc1 z#Uyu-0?A<-38CFv!mQ!weE5aVWB?AqKcsFd16b6p3naC6oT+i{;)|#uH>7`RfV0n0 zj!OQEFn`oc&-0JWD`yI}s)CjvzWqy8yFfPb=(ajgDvg}KT$SR>yU;YSlyxlQn6WDX zq_K3h_$cmuGE@EJ%j;&C`y?*3dEc3R4u5OHlAwx5i48cJVx4NRyXL9sMjMdr0ESEz z1lz+R-Q-sH^pJfBqh_Bfvax@NBeNyB$4T4r)gdE$$_aw2?l!GTCkPX(4#mTMK-67A z>Nu7tsvggWIA;hG%sB$z zvUGUMK^`aGJdSs@JFI`3MNAaQo74QkPl(|nGTdN=MclQyWKAMb^vn=EWZ*opPU3Qv zn7T%MLEgIIDHe*YCF&7*gq7gRgu!k)%wdSC{XZHfBvR>yeXxy+WtVwB+)QE9IJRO) z6A0)+C=r1GjXiaz#IIfW^{&()oJ#QR8f~tCO~^@L+uIM*mb`yWOG2r5S|YPYg;tmq zom6i7VgxYsr}QgXKHW-O4geuujezmz47)$lk{+>v^yTSrKc4DMU;qy-AFEfcXWDyQr6YpLLs>#r;P|yS~pJjas{x)xh0X9 zsTr2728LH?l;?lxaql)Q=4T|(9T_0NtC3?8XkgP5J9b_$h;e@QVF|bgepyqnyF*c% z^5KZ+E0}fPWzu^3=0oncae#iDGAP-oZx5|)U!?cdJ1GLln9jCuCxzumLtSGxnW4Lc zfS6X&kG$4NYkRs?pt6Z#ad`b(f(x#{-5gc-Chcpo-nf6O0Ld!hy62qhKYHarv`AQj z231#8hso7qN}R+~2ZKNJg>y&!hD$G+J!A<8edEYP?$__78p=ON!Do#m<>+vB4fECR%QrES(e11=NNP-nTqef#_W!c3rJf6py+|)|e=ZQ{fss7pBW8mJ zsDdHr67GLYjkVG4Ke!b~+eZ&BUlr6y^BuVO&~%L$!=&=J$2VVJ+4j$|kkG+oaHG*Q zC7s8^khCa;w;>5_2W*bEg6Bevpu3`%BWmuq&jKY|*;m-lTY;>p_BUJZ&k%AN=uRtk z@L9Z~j-3mc7HTCo5}5@bmgR9Z#$3p{g$PbC+pT}Kt&-pi8Y@`n%OW|U-_>@r8Oh0{ z5g93Vws}gFX%@icmUoZAExZ*QFvIQxHnR(zEtTfySx@`6;cUgsm)+&i1}>?U>+BK& z;#T4hdVNxcv?7nt2dRAT&)s?laO~;fN4iD+y1~poTufMaag8 z7&?Cf#WP{@v?mn~@v`9E#BSJxAr1av5n2j@4?e=I=0v*kE7CfU^m709LoCu5$7GKA z3dK)t9Q8XS6g&=ENVhi$rdNE`?iC&r>wAgD##k8zw~w{V&6CL#_^6_a=T z!~!oRt$_La3qkEK%jWePp<0UYRUuxSKP!Jfd&Bays`jJ2it6EOdQ-iolumqZ;@z5QAzdU5Vokx}stlww3Hqi@NkywXx?=D+v)%Yue29SwZ zO57V+J4Hr!_d&(|fr$i1eR9{a?D2MoXQzZPyQZbL)k?s7nM&ZV?hiYqb%6?RMhB6- z3Ka4f0{$Tc`53xFl=AJNAU|BZ*)D%B1z@TTCZ@yK)iwIKuXJJ#jt_MC%A{bmrxH1A zh6$J}5eshUxh&*G*0BGGn&yv%wX1jMR&}E6beX(X%W~;Wd-gT%gMqEAu!5XC=3{)LRSt?tUV!A9Y1sp@~D)c|uiF)ekKQMoy$q+aE z^o;?WHB@AN1S~@e1?jgnuE46NoI5P6&mU#8|~5g&X{>gf4| zuv13XLh5rra_=QBs z7$5q093h@vatueTw)pi@0KqDhX3EW6iM9@7g$cEN z9lw3;1pAMXR%pCzc2tUU$3t;DKR9{dR0<_aFYE^< z_9KtYr6X@IybLAJRw#r|f6cqRK@}`H!dOnPt3x<@dBSp6*FCl^2V&wyMP*;VT|_L{ zS!bSH59$yRoP9wbssMkXeZ&4=ll`a4xQ@pUy1~e$(`y~)2Dr~P;wc)nWcZw!wPDE8 zL52mqP*x*j)U5FK3fHirzM6%0vp;v$ zWmcUj?B1`mgKU=WYQUbNi1A&3iH^6bjN}8~YfXV2Y@N;}o+*E&1?vQE4Hf+BTSy!Q zwC$=_Q&RQ(>kr3s-5$$v`pg1WJsy?@irc88XW#jwXxyuwg;4h?^p4#1iuugB!I(okQA_XgcsFW`zuHyl!+AtTlI{%4 zN1a!0j$1uU7!iMxq}M>Nm$8HJyOp=Ar((z_W`GEL3PWYKEIV$UY_B-t+l2W;l%0GV+`@lDyfN}`qh;2VH2jNYq>3IA z;Uz@+bx{@lo?L!Xjx!y14AAD5ji)_FWOWJrP&J^1SOj%hDRH9#Fg#|sJ$wk9Bdrgi zl#!Fyf9Oit4^+)2FkHKzO66LBXbCz9)|B4tS54m;qEHBZ5JZZ!BRRLs`@)f)NFw1- z@`P&yGpT=on$78}9;cwti{U~XtIIA+D!78C?nQIgy>QvqX!tzhh}^1bwqXX}7e11* z!>8?QU#2hReNIponBDvd8_@cAtX!1}BB4&F^tz=gng1}F>VsV?{c+`|Uhwhj8wIu+ z4%xLr4H5&y>IdO0D8Sh;cIS_(!zjn(SZz@a%ea4w>+*a^Tm}YfN^C zSZU&c4apiID-r=&%x_#s)n{%c20yNxDZ6!?ZKoznUpp`4wZ zkVbbX5_7OMQ4@E^o=|bu`MQbP;UlleA0*^ZDBc$huBfbS*S-og*i*)U(-hz?b%>$C zU_O79s!~(`9rZhXcAnT>N#90pQd^j!)s5udf-s@r8PyUeKG@H5p&9y?(-fZM5`I7} zv1K3@$YOsPD~LoXtZmfLCaBtc5&M{pHgM;Gb-|*QfAzZgEDE~mi%mu}Zv(^TX1M5H zh^o-UwFHm}>fLepcgqOFXFL)%h@k-B_~)L1z!oMVczC71XzVssA-Zc~5T zc^ux2*;%2cn&?Zmj&LP#rvs~k|DS%p7j-v~gEhf(J{r8&Iw0{aAvD#{TLQ#kV@DW2 ziOY%i_k{@Kr0C|@j7vZvri&sW9GWJ#8pciGnyQodC^8eAvMs(q5*4hJ>7!bF^J`Ak zFAR<|HcTe+b>c6l0^-RSk+n03;Mjj-DQm~16xO#Vk1+J?ifA9Sl`~V_bfNq^P(G7I z%RpwG56g=n?{G=f)OZU}Py@{z0yE1i2{=!a_`0ML0)d`gC|S?qQff-jQ35b}bLGel zfe0wcYr03#P^{Fj=UkURXBosA;rruvaP~>2%BEco4^KF>b!bWfaWv+#!`XjNxRa-< zR3T+*1ftqaT<-9MqY2};j4)Z_AhTpmeXLA1#OV$L94kKmjtX(te7fr#b1p+j*Nw7O>-ak>H8kA@M{lu@(8)zbO&06IDsxvG7ai9Z-+BZIEP+CJJ|i zB6xRRi1QzxV6-WTDjkPtu}FV*S7N-S&ZhV!HK}kZh2a+nubEklmxU7 zmURS4YeiR3Ftpt<&q)sSH3Nolz}owFzTpY@!MXqafbT)R-jcK@@bvJM!FY1oD|?hV zl8w6a-9DL+=9=x=0m*+yF;qAOWE4s4e((a~C%0`Cb)+_D+Y$pfFa0RTO$?J?|E?|$ zD;Jv}?Bx5ZVphW~Ck*>Z)~Tgr0(QC!hnUI7>HuMEpCE|T5WQBR<-ANJXR$w~O^s5~ z9NkOqz4rJDa4U{9!k%E;Td-VWD8dK7l!nRgI23F)T+myZ4YPlRs+u5mjbrem-Gnz9 zt@*s;Ww*cwWZ%-C*jREKg3|)GoKkynJJCqRnhw3n;ETR({ z(j`*VRnWAd*LX20P+VS$+mN8Jm%9IOA?&p3FpC;S_fS@l4Otk+`rx*VeG?Ox)7P0S($6B3 zu{f@r9I+v5CI>7pW0q{^b4MP{mt+6{0LP^SA^-u}Pv7$b0LUWHQ9dvYivj=u00045 ET4G=u{r~^~ delta 3154 zcmV-Y46XBkQ0zyL83e02BdxI=r5ArydN88yyA$z+K%-DpjZHHM$F%wv3styYD3Y9U zOV18em)Kb%O*OGPR5#v_JJh!7!nLEVs?SQdh-L`2k)h9ex253)wnh(`Ecy@CA;E)w z6-%L{a)2*K26%E006`sRQ%UBFwkj0s1%D6CR<*wEnsRlBc7M1juMlz}U)+Dg5!OAx zlE~lAUU2`?^H%9N5*S6hp=Dz9-S-|Px{u{)vs4fzX?NEKmY z%SF1ny7BR(R$jzvGRrZmEeSvGb<$cIFKZnkZFpq_UUB*7^~xNJaE8|1!9l`7y+B^@ zD}l3~_kXqZB??Cxx8L9?Hpay@p=9@S&uCchda)m@~yMaSP*v%(;qGKb|VL$Q=W?H>p6 z@6qg$D_#owmF5HbLF3_bag&K5miMP`xO9OzA$lRg`9undc?gy6ZQSw zTZLMWn{=`9Fpr#Du>%ic4~U}z8OU2enCmZj&Nj2321Cc&pO=+>NBw#y@fr-SLCY|` z+z8rvpm;00uj?g93L?ZUij|kprbiM+RQY!p!=2n({-ycrl5~=h2!7O~B$1QeQBR`d zrK(N9WlA9NO;OeC4rYJ+keabNIZ<(m$Fp##uTor#+QsWYEX-8LCM{}8)&+OOMAICy z^RE5HVA20(ME*#&o3Z#4i9nB6EhY}xriRXr+lQ5}`pVZ`j8|9sc!G$NK5mhQG*=g2 zc1_|M`QRNnrvNWtJ&D;Gl@oA^+ahTd2_%_b6p1)3gRG&cWCMTOayp1#@f->5EZf%` z@;JT4@v3(m9~W{sM6xSM&n49xGPP|b9gAGgZ|yQGa$*Soz+SWc7T$`iTxtK8=}<5> z#30RfSj=lGH9IBXiwzN@A5)g>LSHb{a?b*(OQv-_3=KVL2+ zcvW+}-^bNmLkfQxWsFT`kXLA?b45_*4cw<8&d(1@@<>3&gyj^V3;DvI=vvE_RX(}RBv1q9 z;;}1OwuWta3C`au-yP3M1V{m1nAVqa10~i}tv`JU={z!JUzlkbUjZB;Bt}kVl9n@R z&WXQB(pDzFCdF$K@s-m%dEZ*o zQ^tSAq-db|b^web9>tQR@`N7saiPN6g5pLXrQk&`_O1N*XrK_ZJB|h$w6aJ>bE*6F zphDED@4*UgvW);Yz59-|UcN&37^_;n?kN*}89->|n@-CNJ)%d?8|_?0;c^z9x9sRH zh`cV^>SO*fBAilHuHUk$+C;mNyAX@f@W6kZPPQ1vH$`v@*~8ZXeDdJnnj;+m+v31` z9_4iVwV#lw?QJQ{cbiO<*VJa)Q{)w&o=yQ%^O#nf!2cNr;0_}!fa8wA13)eF&uHfJ zdUZ?lI@*s-?^}PHvO}RVd&YldBbQsI-}~ka-fqJaPh(`6FX}=uJM#&{9vs&?aZ7(( z>9Yrd$s_( z@z))QcUrix9gBSF8NiH4g*UP+kk+^XyMv31nVf&S~uu^%rDKV$Nz1Oaz&~F7%B66 zSWqrvLQWR}2y)n%(5Y>@Q7PC4l-bM&1N@3Ywg&}7f&Y;^GGt&U@^eJk0y?>jiIuba zvXX2}OfT82+Y?Pk5(W zxaAw7S$Sl12c0d&JlVt92i(Q3>kMK#+tGBX6^vMDgtnf{|3|<6u8idR8(hWoY(=;s zh9mP)eyMli-K|HGMqYJ*^qU7oiWH-Gw7mAZF}y1xU%Vz8{)?F6Km2r4j0b|^NP*h0 zWWIH*JPG7_~K(Jg^%;M3A_A*)$`9g25H-SvO9J>y}4t4w!qhJ??^ zwJv6s@@nX(A?rS4no>)g3yy?>d}u3l+xL(75wc_5=CUB5lW%2l02L^Mog%W=jQvwSki)S9;|Lj4m#v1>-mQz*cV zBr%>Urdx#>*r%{f?EX$4B+JD@Wp1kKWn8xX*pME>G-gT=@EZEsNc7FwCE*^eoMAh( z#@tj{D}w)E?M-KxCP9A*NDcTpRkI0RM?yRhhC>z3&s~3v{Thei=<4Pc+vS-~%5;Jk2Z;s!l{p@S;!92ad>UMqcCrlA_DV>Pshuw>7 zp+pIO;v;75SwFw>0vf%b$8a@d6-FH1jZ-9PXG*}i_N%X&8w5KZR>b$+m1i$*80V1_ z6+5h}7QlaI>pIe~VrcP!5b;&6*FG54UU?Q+d%zg5a$c6_jRepzc!-NNW~M+*_N93# zjf;LcW+rz>enX9nt?=H$&(at4{LoK?g(|`TMf0_6wyV~$HZLt1|1O)QK&#Bg6^;(mOkf4;J1w?7MoDA0|3XHP#8KPuAB03E6+z#kUN!%j+|xS!46vHL-J zOQL^mDA%3GVeGzSdH>^fEO9f^CW1!JHcV8cBhQ4qBE1ywEJ`bsn$8_F+*vX7l+&fh zz{ARk)JK~asI`&`OWV!7*}fq`n?>;)o3S16(eJw5`i;jy(2rl#R_<}fo)^swd&(bb zS5CgYtL5zuLn0Jpo80yZoXGKWkHq`jaMMCL9kQrSV7{j(!HO3d)d6K-^Aj?yAW)pW sl<4HYZs6YE*Z=?k0R9i48UO*sN4nMm0DW3sp*}DTivj=u00045S|9#2@Bjb+