phenofit shinyapp released

eco-hydro · Dec 17, 2018 · fd12ee7 · fd12ee7
1 parent a5aeaac
commit fd12ee7
Show file tree

Hide file tree

Showing 14 changed files with 5,774 additions and 5,586 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,15 +1,13 @@
 Package: phenofit
 Type: Package
 Title: Extract Remote Sensing Vegetation Phenology 
-Version: 0.1.7
+Version: 0.1.8
 Authors@R: c(
     person("Dongdong", "Kong", role = c("aut", "cre"), email = "kongdd.sysu@gmail.com"),
     person("Yongqiang", "Zhang", role = c("aut"), email = "yongqiang.zhang2014@gmail.com"),
-    person("Qiang", "Zhang", role = c("aut"), email = "zhangq68@bnu.edu.cn"),
-    person("Xihui", "Gu", role = c("aut"), email = "guxh@cug.edu.cn"),
     person("Mingzhong", "Xiao", role = c("aut"), email = "xmingzh@mail2.sysu.edu.cn"),
-    person("Jianjian", "Cui", role = c("aut"), email = "cuijj6@mail2.sysu.edu.cn"),
-    person("Wenlin", "Huang", role = c("aut"), email = "smile.room@qq.com"))
+    person("Xihui", "Gu", role = c("aut"), email = "guxh@cug.edu.cn"),
+    person("Jianjian", "Cui", role = c("aut"), email = "cuijj6@mail2.sysu.edu.cn"))
 Description: 
     The merits of 'TIMESAT' and 'phenopix' are adopted. Besides, a simple and 
     growing season dividing method and a practical snow elimination method 

diff --git a/README.Rmd b/README.Rmd
@@ -38,7 +38,13 @@ You can install phenofit from github with:
 devtools::install_github("kongdd/phenofit")
 ```
 
-## Example
+Or run with `shiny`:
+```{r, eval = FALSE}
+shiny::runGitHub("phenofit", "kongdd", subdir = "inst/shiny/phenofit")
+```
+![](man/Figure/phenofit_shiny.png)
+
+## R code Example
 
 Here, we illustrate how to use `phenofit` to extract vegetation phenology from 
 MOD13A1 in the sampled points. Regional analysis also can be conducted in the 
@@ -183,7 +189,7 @@ brks2 <- season_3y(INPUT, south = sp$lat[1] < 0,
 ## 2.3 Curve fitting
 ```{r curve fitting, fig.height=7, fig.align="center"}
 fit  <- curvefits(INPUT, brks2,
-                  methods = c("AG", "zhang", "beck", "elmore"), #,"klos",, 'Gu'
+                  methods = c("AG", "Zhang", "Beck", "Elmore"), #,"klos",, 'Gu'
                   debug = F, 
                   wFUN = wFUN,
                   nextent = 2, maxExtendMonth = 3, minExtendMonth = 1,

diff --git a/README.md b/README.md
@@ -1,27 +1,21 @@
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
+phenofit
+--------
 
-## phenofit
+[![Travis Build Status](https://travis-ci.org/kongdd/phenofit.svg?branch=master)](https://travis-ci.org/kongdd/phenofit) [![AppVeyor Build Status](https://ci.appveyor.com/api/projects/status/github/kongdd/phenofit?branch=master&svg=true)](https://ci.appveyor.com/project/kongdd/phenofit) [![codecov](https://codecov.io/gh/kongdd/phenofit/branch/master/graph/badge.svg)](https://codecov.io/gh/kongdd/phenofit)
 
-[![Travis Build
-Status](https://travis-ci.org/kongdd/phenofit.svg?branch=master)](https://travis-ci.org/kongdd/phenofit)
-[![AppVeyor Build
-Status](https://ci.appveyor.com/api/projects/status/github/kongdd/phenofit?branch=master&svg=true)](https://ci.appveyor.com/project/kongdd/phenofit)
-[![codecov](https://codecov.io/gh/kongdd/phenofit/branch/master/graph/badge.svg)](https://codecov.io/gh/kongdd/phenofit)
+A state-of-the-art **remote sensing vegetation phenology** extraction package: `phenofit`
 
-A state-of-the-art **remote sensing vegetation phenology** extraction
-package: `phenofit`
+-   `phenofit` combine merits of TIMESAT and phenopix
+-   A simple and stable growing season dividing methods was proposed
+-   Provide a practical snow elimination method, based on Whittaker
+-   7 curve fitting methods and 4 phenology extraction methods
+-   We add parameters boundary for every curve fitting methods according to their ecological meaning.
+-   `optimx` is used to select best optimization method for different curve fitting methods.
 
-  - `phenofit` combine merits of TIMESAT and phenopix
-  - A simple and stable growing season dividing methods was proposed
-  - Provide a practical snow elimination method, based on Whittaker
-  - 7 curve fitting methods and 4 phenology extraction methods
-  - We add parameters boundary for every curve fitting methods according
-    to their ecological meaning.
-  - `optimx` is used to select best optimization method for different
-    curve fitting methods.
-
-## Installation
+Installation
+------------
 
 You can install phenofit from github with:
 
@@ -30,20 +24,26 @@ You can install phenofit from github with:
 devtools::install_github("kongdd/phenofit")
 ```
 
-## Example
+Or run with `shiny`:
+
+``` r
+shiny::runGitHub("phenofit", "kongdd", subdir = "inst/shiny/phenofit")
+```
+
+![](man/Figure/phenofit_shiny.png)
+
+R code Example
+--------------
 
-Here, we illustrate how to use `phenofit` to extract vegetation
-phenology from MOD13A1 in the sampled points. Regional analysis also can
-be conducted in the similar way.
+Here, we illustrate how to use `phenofit` to extract vegetation phenology from MOD13A1 in the sampled points. Regional analysis also can be conducted in the similar way.
 
-## 1.1 Download MOD13A1 data
+1.1 Download MOD13A1 data
+-------------------------
 
-Upload point shapefile into GEE, clip MOD13A1 and download vegetation
-index data.
-[Here](https://code.earthengine.google.com/ee3ec39cf3061374dab435c358d008a3)
-is the corresponding GEE script.
+Upload point shapefile into GEE, clip MOD13A1 and download vegetation index data. [Here](https://code.earthengine.google.com/ee3ec39cf3061374dab435c358d008a3) is the corresponding GEE script.
 
-## 1.2 Initial weights for input data
+1.2 Initial weights for input data
+----------------------------------
 
 Load packages.
 
@@ -59,8 +59,7 @@ suppressMessages({
 })
 ```
 
-Set global parameters for
-`phenofit`
+Set global parameters for `phenofit`
 
 ``` r
 # lambda   <- 5    # non-parameter Whittaker, only suit for 16-day. Other time-scale
@@ -71,24 +70,20 @@ nptperyear <- 23   # How many points for a single year
 wFUN       <- wBisquare #wTSM #wBisquare # Weights updating function, could be one of `wTSM`, 'wBisquare', `wChen` and `wSELF`. 
 ```
 
-Read the point shapefile to get points coordinate information. Read
-Enhanced Vegetation Index (EVI) exported by `GEE`.
+Read the point shapefile to get points coordinate information. Read Enhanced Vegetation Index (EVI) exported by `GEE`.
 
-  - Add date according to composite day of the year (DayOfYear), other
-    than image date.
-  - Add weights according to `SummaryQA`.
+-   Add date according to composite day of the year (DayOfYear), other than image date.
+-   Add weights according to `SummaryQA`.
 
-For MOD13A1, Weights can by initialed by `SummaryQA` band (also suit for
-MOD13A2 and MOD13Q1). We have written a qc function for `SummaryQA`,
-`qc_summary`.
+For MOD13A1, Weights can by initialed by `SummaryQA` band (also suit for MOD13A2 and MOD13Q1). We have written a qc function for `SummaryQA`, `qc_summary`.
 
-| SummaryQA        | Pixel reliability summary QA                        | weight |
-| ---------------- | --------------------------------------------------- | ------ |
-| \-1 Fill/No data | Not processed                                       | `wmin` |
-| 0 Good data      | Use with confidence                                 | 1      |
-| 1 Marginal data  | Useful but look at detailed QA for more information | 0.5    |
-| 2 Snow/ice       | Pixel covered with snow/ice                         | `wmin` |
-| 3 Cloudy         | Pixel is cloudy                                     | `wmin` |
+| SummaryQA       | Pixel reliability summary QA                        | weight |
+|-----------------|-----------------------------------------------------|--------|
+| -1 Fill/No data | Not processed                                       | `wmin` |
+| 0 Good data     | Use with confidence                                 | 1      |
+| 1 Marginal data | Useful but look at detailed QA for more information | 0.5    |
+| 2 Snow/ice      | Pixel covered with snow/ice                         | `wmin` |
+| 3 Cloudy        | Pixel is cloudy                                     | `wmin` |
 
 ``` r
 data('MOD13A1')
@@ -117,11 +112,10 @@ df[, c("w", "QC_flag") := qc_summary(SummaryQA)]
 df <- df[, .(site, y = EVI/1e4, t, date, w, QC_flag)]
 ```
 
-Add one year in head and tail, for growing season dividing. For example,
-the input data period is 20000218 ~ 20171219. After adding one year in
-head and tail, it becomes 19990101 ~ 20181219.
+Add one year in head and tail, for growing season dividing. For example, the input data period is 20000218 ~ 20171219. After adding one year in head and tail, it becomes 19990101 ~ 20181219.
 
-## 2.1 load data
+2.1 load data
+-------------
 
 ``` r
 sites        <- unique(df$site)
@@ -139,15 +133,15 @@ titlestr   <- with(sp, sprintf('[%03d,%s] %s, lat = %5.2f, lon = %6.2f',
 file_pdf   <- sprintf('Figure/%s_[%03d]_%s.pdf', prefix_fig, sp$ID[1], sp$site[1])
 ```
 
-If need night temperature (Tn) to constrain ungrowing season backgroud
-value, NA values in Tn should be filled.
+If need night temperature (Tn) to constrain ungrowing season backgroud value, NA values in Tn should be filled.
 
 ``` r
 d$Tn %<>% zoo::na.approx(maxgap = 4)
 plot(d$Tn, type = "l"); abline(a = 5, b = 0, col = "red")
 ```
 
-## 2.1 Check input data
+2.1 Check input data
+--------------------
 
 ``` r
 dnew  <- add_HeadTail(d, south, nptperyear = 23) # add additional one year in head and tail
@@ -157,18 +151,14 @@ INPUT <- check_input(dnew$t, dnew$y, dnew$w,
 INPUT$y0 <- dnew$y 
 ```
 
-## 2.2 Divide growing seasons
+2.2 Divide growing seasons
+--------------------------
 
-Simply treating calendar year as a complete growing season will induce a
-considerable error for phenology extraction. A simple growing season
-dividing method was proposed in `phenofit`.
+Simply treating calendar year as a complete growing season will induce a considerable error for phenology extraction. A simple growing season dividing method was proposed in `phenofit`.
 
-The growing season dividing method rely on heavily in Whittaker
-smoother.
+The growing season dividing method rely on heavily in Whittaker smoother.
 
-Procedures of initial weight, growing season dividing and curve fitting
-are separated. Phenology extraction and curve fitting are combined
-together.
+Procedures of initial weight, growing season dividing and curve fitting are separated. Phenology extraction and curve fitting are combined together.
 
 ``` r
 par(mar = c(3, 2, 2, 1), mgp = c(3, 0.6, 0))
@@ -208,11 +198,12 @@ brks2 <- season_3y(INPUT, south = sp$lat[1] < 0,
 
 <img src="man/Figure/divide growing season-1.svg" style="display: block; margin: auto;" />
 
-## 2.3 Curve fitting
+2.3 Curve fitting
+-----------------
 
 ``` r
 fit  <- curvefits(INPUT, brks2,
-                  methods = c("AG", "zhang", "beck", "elmore"), #,"klos",, 'Gu'
+                  methods = c("AG", "Zhang", "Beck", "Elmore"), #,"klos",, 'Gu'
                   debug = F, 
                   wFUN = wFUN,
                   nextent = 2, maxExtendMonth = 3, minExtendMonth = 1,
@@ -252,24 +243,24 @@ print(params$AG)
 # # A tibble: 18 x 8
 #    flag      t0    mn    mx    rsp    a3    rau    a5
 #    <fct>  <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl>
-#  1 2000_1  571. 0.168 0.406 0.0256  4.59 0.0173  5.08
-#  2 2001_1  925. 0.174 0.405 0.0215  4.48 0.0168  6   
-#  3 2002_1 1299. 0.189 0.494 0.0310  2.25 0.0178  3.57
-#  4 2003_1 1637. 0.168 0.432 0.0263  2    0.0124  3.71
-#  5 2004_1 2021. 0.181 0.449 0.0395  2    0.0181  4.79
-#  6 2005_1 2415. 0.180 0.462 0.0132  6    0.0332  2   
-#  7 2006_1 2743. 0.175 0.436 0.0198  2.08 0.0135  3.78
-#  8 2007_1 3114. 0.166 0.483 0.0207  2    0.0151  2.89
-#  9 2008_1 3479. 0.175 0.501 0.0265  2    0.0152  3.90
-# 10 2009_1 3887. 0.174 0.485 0.0135  5.25 0.0319  2   
-# 11 2010_1 4203. 0.181 0.493 0.0234  2    0.0147  2   
-# 12 2011_1 4568. 0.184 0.465 0.0310  2    0.0127  4.80
-# 13 2012_1 4925. 0.170 0.511 0.0391  2.57 0.0115  6   
-# 14 2013_1 5317. 0.171 0.489 0.0165  6    0.0153  2.20
-# 15 2014_1 5666. 0.201 0.500 0.0342  2    0.0137  3.57
-# 16 2015_1 6041. 0.217 0.504 0.0231  4.40 0.0210  5.76
-# 17 2016_1 6385. 0.195 0.485 0.0407  2    0.0126  4.76
-# 18 2017_1 6766. 0.175 0.445 0.0210  6    0.0126  3.70
+#  1 2000_1  584. 0.158 0.410 0.0189  5.50 0.0220  3.06
+#  2 2001_1  926. 0.171 0.402 0.0208  4.76 0.0171  6   
+#  3 2002_1 1300. 0.186 0.501 0.0296  2.27 0.0181  2.85
+#  4 2003_1 1639. 0.162 0.447 0.0251  2    0.0128  3.46
+#  5 2004_1 2022. 0.176 0.449 0.0392  2    0.0178  4.55
+#  6 2005_1 2415. 0.175 0.461 0.0131  6    0.0319  2   
+#  7 2006_1 2742. 0.162 0.437 0.0192  2    0.0129  3.65
+#  8 2007_1 3116. 0.156 0.480 0.0189  2    0.0146  2.72
+#  9 2008_1 3478. 0.169 0.502 0.0265  2    0.0148  3.57
+# 10 2009_1 3887. 0.167 0.479 0.0132  5.25 0.0308  2   
+# 11 2010_1 4201. 0.174 0.482 0.0238  2    0.0121  2   
+# 12 2011_1 4569. 0.180 0.465 0.0301  2    0.0131  4.97
+# 13 2012_1 4939. 0.160 0.509 0.0242  3.91 0.0136  4.98
+# 14 2013_1 5290. 0.188 0.482 0.0386  6    0.0112  4.04
+# 15 2014_1 5667. 0.201 0.501 0.0341  2    0.0137  3.51
+# 16 2015_1 6053. 0.211 0.503 0.0178  6.00 0.0289  4.44
+# 17 2016_1 6385. 0.192 0.484 0.0398  2    0.0126  4.78
+# 18 2017_1 6765. 0.170 0.445 0.0213  5.88 0.0123  3.71
 
 ## Get GOF information
 stat  <- ldply(fit$fits, function(fits_meth){
@@ -278,27 +269,27 @@ stat  <- ldply(fit$fits, function(fits_meth){
 fit$stat <- stat
 print(head(stat))
 #   meth   flag       RMSE       NSE         R       pvalue  n
-# 1   AG 2000_1 0.10118415 0.5059028 0.8994179 5.361953e-09 23
-# 2   AG 2001_1 0.08849003 0.5841794 0.9040939 3.322258e-09 23
-# 3   AG 2002_1 0.08954279 0.6682757 0.9018838 1.765036e-09 24
-# 4   AG 2003_1 0.06113231 0.7855540 0.9429517 1.691407e-11 23
-# 5   AG 2004_1 0.06255657 0.7342777 0.9313366 1.124177e-10 23
-# 6   AG 2005_1 0.08182111 0.7225756 0.9266332 1.638453e-09 21
+# 1   AG 2000_1 0.10118884 0.5058570 0.9009524 4.594455e-09 23
+# 2   AG 2001_1 0.09080665 0.5421617 0.9087421 1.372861e-10 26
+# 3   AG 2002_1 0.08280808 0.7139161 0.9045688 7.588185e-09 22
+# 4   AG 2003_1 0.05988494 0.7940177 0.9398957 8.653382e-11 22
+# 5   AG 2004_1 0.06100681 0.7512312 0.9302713 3.670010e-10 22
+# 6   AG 2005_1 0.08166482 0.7178998 0.9281228 4.926381e-10 22
 
 print(fit$fits$AG$`2002_1`$ws)
 # $iter1
-#  [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
-#  [8] 0.2000000 0.2000000 0.5000000 1.0000000 0.9507082 0.9452160 0.1000000
-# [15] 0.5000000 0.6913968 1.0000000 1.0000000 0.1000000 0.2000000 0.1000000
+#  [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.5000000
+#  [8] 1.0000000 0.9192712 0.9979039 0.1000000 0.5000000 0.8826912 1.0000000
+# [15] 1.0000000 0.1000000 0.2000000 0.1000000 0.2000000 0.2000000 0.2000000
 # [22] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
-# [29] 0.2000000 0.2000000 0.2000000 1.0000000 0.8705899
+# [29] 1.0000000 0.7849687
 # 
 # $iter2
-#  [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
-#  [8] 0.2000000 0.2000000 0.5000000 1.0000000 0.9507082 0.9149649 0.2000000
-# [15] 0.5000000 0.6790401 1.0000000 0.9928448 0.2000000 0.2000000 0.2000000
+#  [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.5000000
+#  [8] 1.0000000 0.9192712 0.9739109 0.2000000 0.5000000 0.8736643 1.0000000
+# [15] 0.9935016 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
 # [22] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
-# [29] 0.2000000 0.2000000 0.2000000 1.0000000 0.8705899
+# [29] 1.0000000 0.7849687
 ## visualization
 # svg("Figure1_phenofit_curve_fitting.svg", 11, 7)
 # Cairo::CairoPDF(file_pdf, 11, 6) #
@@ -310,7 +301,8 @@ grid::grid.newpage(); grid::grid.draw(g)# plot to check the curve fitting
 
 <img src="man/Figure/curve fitting-1.svg" style="display: block; margin: auto;" />
 
-## 2.4 Extract phenology.
+2.4 Extract phenology.
+----------------------
 
 ``` r
 # pheno: list(p_date, p_doy)
@@ -350,12 +342,12 @@ head(pheno$doy$AG)
 # # A tibble: 6 x 21
 #   flag  origin     TRS1.sos TRS1.eos TRS2.sos TRS2.eos TRS5.sos TRS5.eos
 #   <fct> <date>        <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
-# 1 2000~ 2000-01-01      163      278      167      273      175      264
-# 2 2001~ 2001-01-01      141      267      146      263      155      254
-# 3 2002~ 2002-01-01      159      278      166      271      180      258
-# 4 2003~ 2003-01-01      124      279      133      270      149      252
-# 5 2004~ 2004-01-01      160      266      167      260      178      251
-# 6 2005~ 2005-01-01      140      274      145      266      156      252
+# 1 2000~ 2000-01-01      163      279      167      273      174      262
+# 2 2001~ 2001-01-01      142      267      146      263      155      254
+# 3 2002~ 2002-01-01      158      284      166      274      180      257
+# 4 2003~ 2003-01-01      125      279      134      270      150      251
+# 5 2004~ 2004-01-01      162      266      168      262      178      251
+# 6 2005~ 2005-01-01      140      274      144      266      156      253
 # # ... with 13 more variables: TRS6.sos <dbl>, TRS6.eos <dbl>,
 # #   DER.sos <dbl>, DER.pop <dbl>, DER.eos <dbl>, GU.UD <dbl>, GU.SD <dbl>,
 # #   GU.DD <dbl>, GU.RD <dbl>, ZHANG.Greenup <dbl>, ZHANG.Maturity <dbl>,

diff --git a/inst/shiny/phenofit/data/phenoflux115.rda b/inst/shiny/phenofit/data/phenoflux115.rda