library(maps); library(mapdata); library(maptools); library(chron); library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
library(fields)
## Loading required package: spam
## Loading required package: grid
## Spam version 1.0-1 (2014-09-09) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## backsolve, forwardsolve
library(raster)
library(rgdal)
## rgdal: version: 0.9-2, (SVN revision 526)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
## Path to GDAL shared files: /usr/share/gdal/1.11
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
library(magic)
## Loading required package: abind
##
## Attaching package: 'magic'
##
## Następujący obiekt został zakryty z 'package:raster':
##
## shift
library(modiscloud)
## Loading required package: date
## Loading required package: sfsmisc
library(openair)
## Loading required package: lazyeval
## Loading required package: dplyr
##
## Attaching package: 'dplyr'
##
## Następujący obiekt został zakryty z 'package:sfsmisc':
##
## last
##
## Następujące obiekty zostały zakryte z 'package:raster':
##
## intersect, select, union
##
## Następujące obiekty zostały zakryte z 'package:stats':
##
## filter, lag
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## intersect, setdiff, setequal, union
library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
##
## Attaching package: 'plyr'
##
## Następujące obiekty zostały zakryte z 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
##
## Następujący obiekt został zakryty z 'package:fields':
##
## ozone
library(dplyr)
#punkty_kasi <- readOGR(dsn="/home/bartosz/Pulpit/kasia/", layer="feno_1_1_centroidy_wgs84")
#pl=readOGR(dsn="/home/bartosz/Dokumenty/projekty/feno2012/gis/shp/", layer="atpol_duze")
#pl_woj=readOGR(dsn="/home/bartosz/Dokumenty/projekty/isok/shp", layer="POL_adm1")
#atpol=readOGR(dsn="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/gis/", layer="atpol_3_3_wgs84")
#rzeki=readOGR(dsn="/home/bartosz/Dokumenty/projekty/isok/shp", layer="rzeki")
#pl_granice=readOGR(dsn="/home/bartosz/Dokumenty/projekty/feno2012/gis/shp/", layer="pl_granice")
#
# #setwd("/home/bartosz/Pulpit/modis/przyciete/aqua/")
# setwd("/media/bartosz/toshiba/przyciete")
# #setwd("/media/bartosz/solarcity/modis/przyciete/")
#
# wynik <- NULL
# pliki=dir(pattern=".tif")
# daty <- unique(substr(pliki,1,8))
#
#
# for (i in 51:length(daty)){
#
# evi_tmp <- paste0(daty[i],".250m_16_days_EVI.tif")
# ndvi_tmp <- paste0(daty[i],".250m_16_days_NDVI.tif")
# doy_tmp <- paste0(daty[i],".250m_16_days_composite_day_of_the_year.tif")
# rel_tmp <- paste0(daty[i],".250m_16_days_pixel_reliability.tif")
#
#
# evi_raster <- raster(evi_tmp)
# evi_raster[which(evi_raster[]==-3000)]=NA
# ndvi_raster <- raster(ndvi_tmp)
# ndvi_raster[which(ndvi_raster[]==-3000)]=NA
# doy_raster <- raster(doy_tmp)
# doy_raster[which(doy_raster[]==-1)]=NA
# rel_raster <- raster(rel_tmp)
# rel_raster[which(rel_raster[]==255)]=NA
#
#
# rok <- as.numeric(substr(daty[i],1,4))
# miesiac <- as.numeric(substr(daty[i],5,6))
#
# tmp <- cbind(as.data.frame(punkty_kasi)[,2:3],extract(ndvi_raster,punkty_kasi), extract(evi_raster,punkty_kasi), extract(doy_raster,punkty_kasi), rok, extract(rel_raster,punkty_kasi))
# names(tmp) <- c("feno","meteo","ndvi","evi","doy","year", "reliability")
#
# if(miesiac==12) tmp$doy <- ifelse(tmp$doy<300,tmp$doy+365,tmp$doy)
#
#
#
# wynik=rbind(wynik,tmp)
# print(i)
#
# if(i%%50==0) saveRDS(wynik,file="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/NDVI_kasia.rds", compress=T)
# }
# saveRDS(wynik,file="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/NDVI_kasia.rds", compress=T)
wynik <- readRDS(file="/home/bartosz/Dokumenty/publikacje własne/2015/aerobiologia/NDVI_kasia.rds")
library(dplyr)
library(hydroTSM)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## Następujące obiekty zostały zakryte z 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: xts
##
## Attaching package: 'xts'
##
## Następujące obiekty zostały zakryte z 'package:dplyr':
##
## first, last
##
## Następujący obiekt został zakryty z 'package:sfsmisc':
##
## last
##
##
## Attaching package: 'hydroTSM'
##
## Następujący obiekt został zakryty z 'package:raster':
##
## extract
wynik$data <- as.Date(paste0(wynik$year,"-01-01"))+wynik$doy-1
wynik$ndvi <- ifelse(wynik$reliability<=2,yes = wynik$ndvi,no=NA)
wynik$evi <- ifelse(wynik$reliability<=2,yes = wynik$evi,no=NA)
time_full=as.data.frame(seq(from=as.Date("2006-12-01"),to=as.Date("2015-01-31"),by="day"))
names(time_full) <- "data"
time <- seq(from=as.Date("2006-12-01"),to=as.Date("2015-01-31"),by="day")
for (i in as.character(unique(wynik$feno))) {
tmp <- subset(wynik,feno==i)
tmp <- join(time_full,tmp,by="data")
proste <- as.data.frame(approx(tmp$data,tmp$ndvi, xout=time))
splajn <- as.data.frame(spline(tmp$data,tmp$ndvi, xout=time,method = "natural"))
tmp <- cbind.data.frame(proste,splajn[,2])
tmp <- subset(tmp,x>=as.Date("2007-01-01") & x<=as.Date("2014-12-31"))
names(tmp)=c("date","ndvi_simpl_1x1","ndvi_splajn_1x1")
# myplot <- timeVariation(tmp,pollutant = c("ndvi_splajn_1x1","ndvi_simpl_1x1"), normalise=F, labels=i)
# plot(myplot, subset="month", xlab="miesiac", legend=F)
#print(i)
barplot(monthlyfunction(tmp[,1:2], mean)/10000, main=i)
}
#
#
#
# tmp2 <- join(time_full,tmp,by="data")
#
# bialy <- subset(wynik,meteo=="Białystok")
# bialy2 <- subset(wynik,meteo=="Białystok")
# bialy$ndvi <- ifelse(bialy$reliability<=2,yes = bialy$ndvi,no=NA)
# bialy$evi <- ifelse(bialy$reliability<=2,yes = bialy$evi,no=NA)
# plot(bialy$data,bialy$ndvi/10000, type='l',col='red')
# lines(bialy2$data,bialy$ndvi/10000, type='p',col='black',cex=0.4)
# #lines(bialy$data,bialy$evi/10000, type='l', col="red")
#
# library(ggplot2)
# qplot(bialy$data,bialy$evi/10000, type='l', col="red")
#
# feno <- readRDS("~/Pulpit/feno.rds")
#