pressure_altitude (R)

#install.packages("readxl")
#install.packages("ggpubr")
#install.packages("ggplot2")
#install.packages("mapdata")
#install.packages("maptools")
#install.packages("plotly")

library(readxl)
library(ggplot2)
library(ggpubr)
library(plyr)
library(plotly)
library(dplyr)
library(reshape2)

# stations_near_airports.xlsx, Airport: basic information of each airport
 Airportinfo = read_excel("Airport1.xlsx", sheet = "Airport")
 APname=Airportinfo$IATA

######################## Pressure Altitude ########
PA = read_excel("new.xlsx", sheet = "PressureAltitude") # unit:m;1 m = 3.28 ft.
PA78 = PA[PA$Month=="7" | PA$Month=="8",]  # July & August
PA78p <- PA78[,1:31]
names(PA78p)[1] <- "Time(day)"
mup <- PA78p %>% group_by(`Time(day)`) %>% summarise_all(funs(mean)) %>% data.frame()
tmelt <- melt(PA78p)

cPA <- ggplot(tmelt, aes(value)) +
  # Histogram  plot
  geom_histogram(aes( color = `Time(day)`,fill = `Time(day)`),  
                 position = "dodge") +  
  scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  scale_fill_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  facet_wrap(~variable,ncol = 4,scales="free")+
  labs(x="Pressure altitude (m)", y = "Count")+theme(axis.title= element_text(size=40,vjust=0.5, hjust=0.5))+
  theme(title= element_text(size=30, vjust=0.5, hjust=0.5))+theme(axis.text= element_text(size=28,  vjust=0.5, hjust=0.5))+
  theme(legend.text = element_text(size =26))+theme(strip.text.x  =  element_text(size=32))+theme(legend.key.size=unit(1.2,'cm'))+
  theme(plot.margin =unit(c(4, 4, 1, 1), "lines"))
cPA + theme(legend.justification=c(0.8,0), legend.position=c(0.8,0))
ggsave("Cpamulti-panel.pdf", width = 27, height =35)
###########################

############# Temperature ##########
T = read_excel("new.xlsx", sheet = "Temperature") # unit:m;1 m = 3.28 ft.

Twinter = T[T$Season=="Winter",] # winter season
Tsummer = T[T$Season=="Summer",] # summer season
T78 = T[T$Month=="7" | T$Month=="8",]  # July & August

################ Temperature multi-panel plot ###############
T78p <- T78[,1:31]
names(T78p)[1] <- "Time(day)"
mut <- T78p %>% group_by(`Time(day)`) %>% summarise_all(funs(mean)) %>% data.frame()
tmelt <- melt(T78p)

cT <- ggplot(tmelt, aes(value)) +
  # Histogram  plot
  geom_histogram(aes( color = `Time(day)`,fill = `Time(day)`),  
                 position = "dodge") +  
  scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  scale_fill_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  facet_wrap(~variable, ncol = 4, scales="free")+
  labs(x="Temperature ( degree C)", y = "Count")+theme(axis.title= element_text(size=40,vjust=0.5, hjust=0.5))+
  theme(title= element_text(size=30, vjust=0.5, hjust=0.5))+theme(axis.text= element_text(size=29,  vjust=0.5, hjust=0.5))+
  theme(legend.text = element_text(size =26))+theme(strip.text.x  =  element_text(size=32))+theme(legend.key.size=unit(1.2,'cm'))+
  theme(plot.margin =unit(c(4, 4, 1, 1), "lines"))
cT + theme(legend.justification=c(0.8,0), legend.position=c(0.8,0))
ggsave("Ctmulti-panel.pdf", width = 25, height = 35)
###################################################
################ Takeoff factor ##############
TF = read_excel("new.xlsx", sheet = "TakeoffFactor") # unit:m;1 m = 3.28 ft.
TFwinter = TF[TF$Season=="Winter",] # winter season
TFsummer = TF[TF$Season=="Summer",] # summer season
TF78 = TF[TF$Month=="7" | TF$Month=="8",]  # July & August

################ Takeoff distance factor multi-panel plot ###############
TF78p <- TF78[,1:31]
names(TF78p)[1] <- "Time(day)"
mutf <- TF78p %>% group_by(`Time(day)`) %>% summarise_all(funs(mean)) %>% data.frame()
tfmelt <- melt(TF78p)

cTF <- ggplot(tfmelt, aes(value)) +
  # Histogram  plot
  geom_histogram(aes( color = `Time(day)`,fill = `Time(day)`),  
                 position = "dodge") +  
  scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  scale_fill_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  facet_wrap(~variable, ncol =4,scales="free")+
  labs(x="Takeoff factor", y = "Count")+theme(axis.title= element_text(size=36,vjust=-1, hjust=0.5))+
  theme(title= element_text(size=30, vjust=0.5, hjust=0.5))+theme(axis.text= element_text(size=29,  vjust=0.5, hjust=0.5))+
  theme(legend.text = element_text(size =26,vjust=0.5, hjust=0.5))+theme(strip.text.x  =  element_text(size=32))+theme(legend.key.size=unit(1.2,'cm'))+
  theme(plot.margin =unit(c(4, 4, 6, 1), "lines"))+theme(legend.title = element_text(size =34,vjust=0.5, hjust=0.5))
cTF + theme(legend.justification=c(0.8,0), legend.position=c(0.8,0))
#tfmelt$Region[ tfmelt$variable %in% north ] <- "North"
#tfmelt$Region[ tfmelt$variable %in% northeast ] <- "Northeast"
#tfmelt$Region[ tfmelt$variable %in% east ] <- "East"
#tfmelt$Region[ tfmelt$variable %in% southcentral ] <- "Southcentral"
#tfmelt$Region[ tfmelt$variable %in% southwest ] <- "Southwest"
#tfmelt$Region[ tfmelt$variable %in% northwest ] <- "Northwest"
ggsave("Tfmulti-panel.pdf", width = 25, height = 35)

####################################################
############### Climb factor ############
CF = read_excel("new.xlsx", sheet = "ClimbFactor") # unit:m;1 m = 3.28 ft.

CFwinter = CF[CF$Season=="Winter",] # winter season
CFsummer = CF[CF$Season=="Summer",] # summer season
CF78 = CF[CF$Month=="7" | CF$Month=="8",]  # July & August

################ Climb factor multi-panel plot ###############
CF78p <- CF78[,1:31]
names(CF78p)[1] <- "Time(day)"
mucf <- CF78p %>% group_by(`Time(day)`) %>% summarise_all(funs(mean)) %>% data.frame()
cfmelt <- melt(CF78p)

cCF <- ggplot(cfmelt, aes(value)) +
  # Histogram  plot
  geom_histogram(aes( color = `Time(day)`,fill = `Time(day)`),  
                 position = "dodge") +  
  scale_color_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  scale_fill_manual(values = c("#999999", "#E69F00", "#56B4E9")) +
  facet_wrap(~variable, ncol =4,scales="free")+
  labs(x="Climb factor", y = "Count")+theme(axis.title= element_text(size=40,vjust=0.5, hjust=0.5))+
  theme(title= element_text(size=30, vjust=0.5, hjust=0.5))+theme(axis.text= element_text(size=29,  vjust=0.5, hjust=0.5))+
  theme(legend.text = element_text(size =26))+theme(strip.text.x  =  element_text(size=32))+theme(legend.key.size=unit(1.2,'cm'))+
  theme(plot.margin =unit(c(4, 4, 6, 1), "lines"))+theme(legend.title = element_text(size =34,vjust=0.5, hjust=0.5))
cCF + theme(legend.justification=c(0.8,0), legend.position=c(0.8,0))
ggsave("Cfmulti-panel.pdf", width = 25, height = 35)

#####################################################
# pressure altitude
m_PA = aggregate(PA78[, 2:31], list(PA78$`Time(day)`), FUN=mean, na.rm=TRUE)
md_PA = aggregate(PA78[, 2:31], list(PA78$`Time(day)`), FUN=median, na.rm=TRUE)
mx_PA = aggregate(PA78[, 2:31], list(PA78$`Time(day)`), FUN=max, na.rm=TRUE)
mn_PA = aggregate(PA78[, 2:31], list(PA78$`Time(day)`), FUN=min, na.rm=TRUE)
# temperature
m_t = aggregate(T78[, 2:31], list(T78$`Time(day)`), FUN=mean, na.rm=TRUE)
md_t = aggregate(T78[, 2:31], list(T78$`Time(day)`), FUN=median, na.rm=TRUE)
mx_t= aggregate(T78[, 2:31], list(T78$`Time(day)`), FUN=max, na.rm=TRUE)
mn_t = aggregate(T78[, 2:31], list(T78$`Time(day)`), FUN=min, na.rm=TRUE)
# takeoff factor
m_takeoff = aggregate(TF78[, 2:31], list(TF78$`Time(day)`), FUN=mean, na.rm=TRUE)
md_takeoff = aggregate(TF78[, 2:31], list(TF78$`Time(day)`), FUN=median, na.rm=TRUE)
mx_takeoff = aggregate(TF78[, 2:31], list(TF78$`Time(day)`), FUN=max, na.rm=TRUE)
mn_takeoff = aggregate(TF78[, 2:31], list(TF78$`Time(day)`), FUN=min, na.rm=TRUE)
# climb rate factor
m_climb = aggregate(CF78[, 2:31], list(CF78$`Time(day)`), FUN=mean, na.rm=TRUE)
md_climb = aggregate(CF78[, 2:31], list(CF78$`Time(day)`), FUN=median, na.rm=TRUE)
mx_climb = aggregate(CF78[, 2:31], list(CF78$`Time(day)`), FUN=max, na.rm=TRUE)
mn_climb = aggregate(CF78[, 2:31], list(CF78$`Time(day)`), FUN=min, na.rm=TRUE)

library(sp)
library(maps)
library(mapdata)
#install.packages("rgeos")
library(maptools)

########### Location of Airports #############
world=map_data('world')

##########
library(grid)
library(gridExtra)

PA1=m_PA[1,2:31]
PA2=m_PA[2,2:31]
PA3=m_PA[3,2:31]

p00 <- data.frame(APname,Airportinfo$Airport_Longitude,Airportinfo$Airport_Latitude,t(PA1),t(PA2-PA1),t(PA3-PA2))
names(p00) <- c("Airport","Longitude","Latitude","1976-2005","from 2021-2050 to 1976-2005", "from 2071-2100 to 2021-2050")
# 1976-2005
mymap1 = ggplot(world) +
  geom_polygon(aes(x = long, y = lat, group = group), colour = "gray",fill = NA) +
  geom_point(data=p00, aes(x=Longitude, y=Latitude, colour=`1976-2005`)) +
  scale_colour_gradientn(colours = rainbow(4),name="m") + # terrain.colors(10)
  theme_bw()+ theme(axis.text = element_text(size=12),axis.title = element_text(size=12),title = element_text(size=12))+
  theme(legend.text = element_text(size=12))+
  labs(x="Longitude", y = "Latitude",subtitle = "Pressure altitude (1976-2005)") +
  #   geom_text(aes(x=Airport_Longitude,y=Airport_Latitude-1,label=IATA), data=Airportinfo) +
  coord_map(xlim = c(-180,180),ylim = c(0,55))
# mymap1
# 2021-2050
mymap2 = ggplot(world) +
  geom_polygon(aes(x = long, y = lat, group = group), colour = "gray",fill = NA) +
  geom_point(data=p00, aes(x=Longitude, y=Latitude, colour=`from 2021-2050 to 1976-2005`)) +
  scale_colour_gradientn(colours = rainbow(5),name="m")+#,limits = c(1.1, 3.2)) + # colours=rainbow(2)  #colours = heat.colors(10)
  theme_bw()+ theme(axis.text = element_text(size=12),axis.title = element_text(size=12),title = element_text(size=12))+
  theme(legend.text = element_text(size=12))+
  labs(x="Longitude", y = "Latitude",subtitle = "Pressure altitude change from 2021-2050 to 1976-2005") +
  #  geom_text(aes(x=Airport_Longitude,y=Airport_Latitude-1,label=IATA), data=Airportinfo) +
  coord_map(xlim = c(-180,180),ylim = c(0,55))
# mymap2
# 2071-2100
mymap3 = ggplot(world) +
  geom_polygon(aes(x = long, y = lat, group = group), colour = "gray",fill = NA) +
  geom_point(data=p00, aes(x=Longitude, y=Latitude, colour=`from 2071-2100 to 2021-2050`)) +
  scale_colour_gradientn(colours = rainbow(5),name="m")+#,limits = c(1.1, 3.2)) + # colours=rainbow(2)  
  theme_bw()+ theme(axis.text = element_text(size=12),axis.title = element_text(size=12),title = element_text(size=12))+
  theme(legend.text = element_text(size=12))+
  labs(x="Longitude", y = "Latitude",subtitle = "Pressure altitude change from 2071-2100 to 2021-2050") +
  #  geom_text(aes(x=Airport_Longitude,y=Airport_Latitude-1,label=IATA), data=Airportinfo) +
  coord_map(xlim = c(-180,180),ylim = c(0,55))
# mymap3

# install.packages('mapproj')
mymap=grid.arrange(mymap1, mymap2, mymap3, ncol = 1)
ggsave("PressureAltitudeChange.pdf",mymap,width = 11,height = 8)