> data1 = read.csv("NDVI_veg.csv", header = T)
> data2 = read.csv('NDVI_all.csv',header = T)
> tell_type <- function(x){
+   b = 0
+   if(x == 'PA'){
+     b = 2
+   }else if(x == 'SA'){
+     b = 4
+   }else if(x == 'SM'){
+     b = 6   
+   }else{
+     b = 8
+   }
+ }
>
> data1$type = sapply(data1$types,tell_type)
> labels = unique(data1$types)
> df1 = data.frame(NDVI = data2$NDVI,season = data2$season,type = NA)
> df2 = data.frame(NDVI = data1$ndvi,season = tolower(data1$season),type = data1$type )
> df = rbind(df1,df2)
>
> p =ggplot()+
+   geom_density(data = df,aes(x = NDVI,y= ..density..,fill = season),alpha = 0.5)+
+   geom_point(data = df,aes(x = NDVI,y = type,color = factor(type)),size = 4)+
+   scale_y_continuous(
+     breaks = c(2,4,6,8),
+     sec.axis = sec_axis( ~ . +0, breaks = c(2,4,6,8),labels = labels,
+                          name = "Type")
+   )+
+   scale_color_manual(values = c('red','yellow','blue','green'))+
+   theme_bw()+
+   theme(
+     axis.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
+     axis.title = element_text(face = 'bold',color = 'black',size = 14,hjust = 0.5),
+     legend.text = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
+     legend.title = element_text(face = 'bold',color = 'black',size = 12,hjust = 0.5),
+     legend.position = 'bottom',
+     legend.direction = 'horizontal'
+   )+
+   ylab("NDVI")+
+   xlab("Density")
> png('plot1.png',
+     height = 20,
+     width = 20,
+     units = 'cm',
+     res = 800)
> print(p)
Warning message:
Removed 33009 rows containing missing values (geom_point).
> dev.off()
png
  2