-2

I want to plot single scatter plot showing a different color for two station months rainfall. for example A station Jan rainfall value in red, Feb in yellow and B station Jan rainfall value in blue, Feb in green and so on which appears on the legend. Also, I want to include a smooth line for the both stations data that also appear on the legend like the red smooth line for A station and blue for station B. In this link, you can find the both stations CSV data: https://drive.google.com/file/d/0B3fQ9_46L-O0TjJwYmF6UThNSGs/view?usp=sharing https://drive.google.com/file/d/0B3fQ9_46L-O0ZXVYb3lzZDBZaHM/view?usp=sharing

Below is the code I tried but could not succeed.

#reading csv file of ramoili station of rautahat[Scatterplot of two stations][1]
ram = read.csv('preci_ramoili.csv',header=TRUE, stringsAsFactors=FALSE)
#reading CSV file of gaur station of rautahat
gaur= read.csv('preci_Gaur.csv',header=TRUE, stringsAsFactors=FALSE)
#gaur rainfall
rain <- data.frame(index(agg),stack(as.data.frame(coredata(agg))))
rain
head(rain)
tail(rain)
names(rain)[1] <- "Year"
names(rain)[2] <- "Rainfall"
names(rain)[3] <- "Month"
#ramoili rainfall
rain1<-data.frame(index(core),stack(as.data.frame(coredata(core))))
rain1
head(rain1)
names(rain1)[1] <- "Year"
names(rain1)[2] <- "Rainfall"
names(rain1)[3] <- "Month"
head(rain1)
#ramoili premonsoon rainfall
rain1_pre<-data.frame(index(core[,3:5]),stack(as.data.frame(coredata(core[,3:5]))))
head(rain_pre)
tail(rain1_pre)
names(rain1_pre)[1] <- "Year"
names(rain1_pre)[2] <- "Rainfall"
names(rain1_pre)[3] <- "Month"
#ggplot of two stations gaur and ramoili yearly rainfall of rautahat in same plot
p9 <- ggplot(rain, aes(x =Year, y=Rainfall, size=Rainfall)) + geom_point(shape = 21,color = "#000000", fill = "#40b8d0") + 
  geom_smooth(aes(fill="Gaur"), colour="darkblue", size=1)

p10 <- p9 + geom_point(data=rain1, aes(x =Year, y=Rainfall, color=Month )) + 
  geom_smooth(data=rain1, aes(fill="Ramoili"), colour="red", size=1)+ 
  ggtitle(" Yearly rainfall at two stations of Rautahat")+
  scale_fill_manual(name="Stations", values=c("blue", "red"))
print(p10)
  • The data helps, but your question and code is missing references to `agg` and `core`. It also helps to list the packages your are using in your sample code. – Jake Kaupp Dec 20 '16 at 15:27

1 Answers1

-2

In the absence of complete sample data, and using the sample data you provided, I've illustrated an approach.

I'm uncertain as to why you want the different month colors for each station, and I think that difference would be illustrated better using facets. If not, I would still recommend keeping the month colors consistent and remove the facet grid.

You would need to modify the axes, titles, scales, etc. to fit your liking.

library(dplyr)
library(ggplot2)
library(purrr)
library(tidyr)
library(trend)

gaur <- list.files("~/Desktop", pattern = "Gaur", full.names = TRUE)  %>% 
  read.csv() %>% 
  mutate(station = "Gaur")

ramoili <- list.files("~/Desktop", pattern = "ramoili", full.names = TRUE)  %>% 
  read.csv() %>% 
  mutate(station = "Ramoli")

plot_data <- bind_rows(gaur, ramoili) %>% 
  gather(month, rainfall, -Year, -station) 


ggplot(plot_data, aes(x = Year, y = rainfall)) +
  geom_line(aes(color = month)) +
  geom_point(aes(color = month), show.legend = FALSE) +
  geom_smooth(aes(fill = station), size = 0.1) +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  scale_fill_manual(name = "Stations", values = c("blue", "red")) +
  facet_grid(month ~ station) +
  theme_minimal()

enter image description here

And for models:

models <- bind_rows(gaur, ramoili) %>% 
      select(-Year) %>% 
      nest(-station) %>% 
      mutate(ts_data = map(data, ~ts(.x, frequency = 1, start = c(1984,1)))) %>% 
      mutate(mk_model = map(ts_data, mk.test),
             sens_slope = map(ts_data, sens.slope))

> models$mk_model[1]
[[1]]
Mann-Kendall Test

two-sided homogeinity test
H0: S = 0 (no trend)
HA: S != 0 (monotonic trend)

Statistics for total series
      S     varS    Z    tau     pvalue
1 -1637 39765.67 -8.2 -0.275 2.2289e-16

> models$sens_slope[1]
[[1]]

Sen's slope and intercept


slope:  0
95 percent confidence intervall for slope
0 0

intercept: 14.9
nr. of observations: 384
Jake Kaupp
  • 7,892
  • 2
  • 26
  • 36
  • I am trying to find the trend of these two stations. Also, test its trend using Mann Kendall and magnitude using sen's slope. Please, could you give me some direction how I should proceed for that on this two datasets. – Pawan Thapa Dec 21 '16 at 11:49
  • I'd look into the `kendall` or `trend` packages. Remeber to accept/up vote helpful answers. – Jake Kaupp Dec 21 '16 at 14:10
  • I added an approach to the modelling as well, using the `trend`, `purrr` and `tidyr` packages. – Jake Kaupp Dec 22 '16 at 00:59
  • Thank you again@Jake Kaupp. While running the model I got an error which is: – Pawan Thapa Dec 23 '16 at 03:20
  • > models <- bind_rows(gaur, ramoili) %>% + select(-Year) %>% + nest(-station) %>% + mutate(ts_data = map(data, ~ts(.x, frequency = 1, start = c(1984,1)))) %>% + mutate(mk_model = map(ts_data, mk.test), + sens_slope = map(ts_data, sens.slope)) Error: missing value where TRUE/FALSE needed In addition: Warning messages: 1: In data.matrix(data) : NAs introduced by coercion 2: In data.matrix(data) : NAs introduced by coercion – Pawan Thapa Dec 23 '16 at 03:23