#zabawa z pakietem hydroTSM
#całość z
#https://cran.r-project.org/web/packages/hydroTSM/vignettes/hydroTSM_Vignette-knitr.pdf
#
#tu jest instrukcja do pakietu
#https://cran.r-project.org/web/packages/hydroTSM/hydroTSM.pdf

install.packages("hydroTSM")
library(hydroTSM)

#wczytanie konkretnych danych które są "na stanie" pakietu
data(SanMartinoPPts)

#jako dane pakiet stosuje obiekty typu zoo

x.Date <- as.Date("2003-02-01") + c(1, 3, 7, 9, 14) - 1
x <- zoo(rnorm(5), x.Date)
plot(x)
time(x)


dane=data(SanMartinoPPts)
plot(dane)
plot(SanMartinoPPts)

#wycinamy troszkę danych
x <- window(SanMartinoPPts, start=as.Date("1985-01-01"))

#i robimy miesięczne sumy
m <- daily2monthly(x, FUN=sum)
m <- daily2monthly(x, FUN=mean)
plot(m)
(m <- daily2monthly(x, FUN=sum))

#sprawdźmy ile tam jest dat
dates <- time(x)
plot(dates)

#oraz ile jest pełnych lat w tych danych
( nyears <- yip(from=start(x), to=end(x), out.type="nmbr" ) )

smry(x)

hydroplot(x, var.type="Precipitation", main="at San Martino",
          pfreq = "dm", from="1987-01-01")


#odpowiedni R tabble, ilośc dni z danymi
dwi(x)
dwi(x, out.unit="mpy")



#mozaika, miesięczna suma opadóW
M <- matrix(m, ncol=12, byrow=TRUE)
colnames(M) <- month.abb
rownames(M) <- unique(format(time(m), "%Y"))

require(lattice) #widzimy że do czegoś bardziej zaawansowanego pakiet wymaga zewnętrznej pomocy
print(matrixplot(M, ColorRamp="Precipitation",
                 main="Monthly precipitation at San Martino st., [mm/month]"))


#działanie na wartościach rocznych
daily2annual(x, FUN=sum, na.rm=TRUE)
mean( daily2annual(x, FUN=sum, na.rm=TRUE) )
#miesięczna mediana, pamiętamy że działamy na sumie opadów dziennych
monthlyfunction(m, FUN=median, na.rm=TRUE)

#automatyczne zrobienie nazw miesięcy 
cmonth <- format(time(m), "%b")

#robienie boxplotu z opadem, w sumie robimy go standardowym poleceniem
#szeregujemy nazwy miesięcy
months <- factor(cmonth, levels=unique(cmonth), ordered=TRUE)
boxplot( coredata(m) ~ months, col="lightblue", main="Monthly Precipitation",
         ylab="Precipitation, [mm]", xlab="Month")


#sezonowo
seasonalfunction(x, FUN=sum, na.rm=TRUE) / nyears
#ciekawostka, w Polinezji Francuskiej wiosna i jesień trwają krócej
seasonalfunction(x, FUN=sum, na.rm=TRUE,type="FrenchPolynesia") / nyears
DJF <- dm2seasonal(x, season="DJF", FUN=sum) 
hydroplot(x, pfreq="seasonal", FUN=sum, stype="default")





data(SanMartinoPPts)
x <- window(SanMartinoPPts, start=as.Date("1988-01-01"))
hydroplot(x, ptype="ts", pfreq="o", var.unit="mm")

( R10mm <- length( x[x>10] ) )
wet.index <- which(x >= 1)
( PRwn95 <- quantile(x[wet.index], probs=0.95, na.rm=TRUE) )
(very.wet.index <- which(x >= PRwn95))
( R95p <- sum(x[very.wet.index]) )


#suma 5-cio dniówek
x.5max <- rollapply(data=x, width=5, FUN=sum, fill=NA, partial= TRUE,
                    align="center")
hydroplot(x.5max, ptype="ts+boxplot", pfreq="o", var.unit="mm")

(x.5max.annual <- daily2annual(x.5max, FUN=max, na.rm=TRUE))


#klimatogram - bardzo lubię to słowo

data(MaquehueTemuco)
plot(MaquehueTemuco$pcp)

#serie danych opad, temperatura min, temperatura max
pcp <- MaquehueTemuco[, 1]
tmx <- MaquehueTemuco[, 2]
tmn <- MaquehueTemuco[, 3]

m <- climograph(pcp=pcp, tmx=tmx, tmn=tmn, na.rm=TRUE)



#teraz troszkę spatial
############
## Loading the monthly time series of precipitation within the Ebro River basin.
data(EbroPPtsMonthly)
## Loading the gis data
data(EbroPPgis)
## Putting the measurements of the first row of 'EbroPPtsMonthly' into their
## corresponding spatial location given by 'x.gis'
require(sp)
x.spt <- gists2spt(x.ts=EbroPPtsMonthly[1,], x.gis=EbroPPgis, X="EAST_ED50",
                   Y="NORTH_ED50", na.rm=FALSE, sname="ID")
## Plotting the measured values (only the first row of 'EbroPPtsMonthly') at their
## corresponding spatial location
spplot(x.spt, zcol="value")




hydropairs(EbroPPtsMonthly[,2:4])

install.packages("rgdal")

require(rgdal)
data(EbroDEM1000m)
dem <- EbroDEM1000m
# Plotting the DEM

require(sp)
spplot(dem, scales=list(draw=TRUE, y=list(rot=90)))
# Computing and plotting the hypsometric curve
hypsometric(dem)
# If the raster file has more than 1 band, and the elevation data are in a
# band different from 1:
dem$ELEVATION <- EbroDEM1000m$band1 # dummy example
hypsometric(dem, band= 2)
hypsometric(dem, band= "ELEVATION") # same as before, but user-friendly