setwd("~/Rclass/")설정된 작업 디렉토리를 확인하는 방법은 ‘getwd()’ 입니다.
getwd()## [1] "D:/OneDrive - hufs.ac.kr/Workspace_hufs/Rclass/Ehwa_240227"위와 같이 “D:/Rclass/”가 작업 디렉토리로 설정된 것을 볼 수 있습니다.
다음으로는 분석에 필요한 라이브러리를 로드해줍니다.
library("lubridate")##
## 다음의 패키지를 부착합니다: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary("plyr")
library("openair")rm(list=ls())
# library("epitools")
library("lubridate")
library("plyr")
setwd("D:/HUFS/Project/2023_NIER_intensive/")
source("../../../Dropbox/Rfunction/f_func.R")
solang2 = function(DOY1, HR1, lat, lon, tdf){
pi = 3.14159265358979323846
timelon= 15.0 * tdf #utc+4h
lat1 = lat * pi/180
lon1 = lon * pi/180
A = 2 * pi * DOY1 / 365
declination = 0.006918 - 0.399912*cos(A)+0.070257*sin(A) -0.006758*cos(2*A) + 0.000907*sin(2*A) - 0.002697*cos(3*A) + 0.001480 * sin(3*A)
localtime = HR1 + (lon-timelon)/15
B = 2 * pi * DOY1/365
et = (0.000075 + 0.001868 * cos(B) -0.032077*sin(B) -0.014615 * cos(2*B) -0.040849 * sin(2*B) ) * 12 / pi
t = 15 * (pi/180) * ( localtime + et - 12 )
coszenith = sin(declination)*sin(lat1) + cos(declination)*cos(lat1)*cos(t)
elevation = 0.5*pi - acos(coszenith)
sinkai = cos(declination)*sin(t)/sin(0.5*pi-elevation)
coskai = (-cos(lat1)*sin(declination)+sin(lat1) * cos(declination)*cos(t))/sin(0.5*pi-elevation)
kai = asin( sinkai )
if( coskai < 0 ) {
kai1 = pi - kai
} else if( coskai > 0 & sinkai < 0 ) {
kai1 = 2*pi+kai
} else {
kai1 = kai
}
kai1 = kai1+pi
if(kai1>2*pi) kai1 = kai1 - 2*pi
azimuth = kai1
return(c(elevation, azimuth))
}
day_of_year = function( yy, mm, dd){
md =c(31,28,31,30,31,30,31,31,30,31,30,31)
if( yy %% 4 == 0 ){
if( yy %% 100 == 0 ){
if( yy %% 400 == 0 ){
md[2]=29
} else {
md[2]=28
}
} else {
md[2]=29
}
} else {
md[2]=28
}
doy = 0
for (k in 1:(mm-1)) doy = doy+md[k]
doy = doy + dd
if (mm == 1) doy = dd
return (doy)
}
opt = read.csv(file = "DATA/OUT/hn_aod.csv")
ID = which(substr(opt$Time,12,20) == "")
opt$Time[ID] = paste(opt$Time[ID], "00:00:00")
opt$Time = as.POSIXct(opt$Time, format = "%Y-%m-%d %H:%M:%S") #- hours(9)
fact = 3.1415/180
ln = nrow(opt)
for (ln in 1:nrow(opt)){
opt$DOY[ln] = day_of_year(year(opt$Time)[ln], month(opt$Time)[ln], day(opt$Time)[ln])
opt$HR[ln] = hour(opt$Time[ln]) + minute(opt$Time[ln])/60 + second(opt$Time[ln])/3600
opt$DOY[ln] = opt$DOY[ln] + opt$HR[ln] / 24
opt$SZA[ln] = 90.0 - solang2(trunc(opt$DOY)[ln], opt$HR[ln], 35.2266035, 126.8489257, 0)[1]/fact
}
plot(opt$SZA ~ opt$Time)
opt$SZA = round(opt$SZA,2)
opt$DOY1 = trunc(opt$DOY)
opt$HR = hour(opt$Time)
write.csv(file = "DATA/OUT/hn_aod1.csv", opt[opt$SZA <= 80, c("SZA", "DOY1", "HR", "AOD", "S450", "SAE")], row.names = F)
###################
rm(list=ls())
# library("epitools")
library("lubridate")
library("plyr")
setwd("/mnt/d/HUFS/Project/2023_NIER_intensive/")
source("../../../Dropbox/Rfunction/f_func.R")
opt = read.csv(file = "DATA/OUT/hn_aod.csv")
opt$Time = as.POSIXct(opt$Time)
plot(opt$SZA ~ opt$Time)
ln = 1
for (ln in 1:nrow(opt)){
descript = paste0("
&INPUT
isat=0, !user defined wavelength range with no filter function
wlinf=0.55, !starting wavelength of 0.25 um
wlsup=0.55, !ending wavelength of 4.0 um
wlinc=-0.01, !wavelength increment is dL/L=.01
idatm=6, !atmosphere: 6 is US standard
nf=3, !use the more accurate MODTRAN3 solar spectrum
sza =", opt$SZA[ln],", !solar zenith angle (not used if zero)
iday=", trunc(opt$DOY[ln]),", !day of the year
time=", hour(opt$Time[ln]), ", !time in hours GMT
alat=35.2266035, !latitude on Earth
alon=126.8489257, !longitude on Earth
isalb=0,
albcon=0.3,
uw=1.418, !precipitable water vapor (-1 means use atmosphere); for US standard atmosphere, uw=1.418(g/cm2)
uo3=0.349, !total ozone (-1 means use atmosphere); US standard atmos. uo3=0.349 (atm-cm)
xco2=360, !carbon dioxide concentration (ppmv)
iaer = 5,
tbaer=", opt$AOD[ln], ",
wlbaer=.440,
qbaer=0.109,
gbaer=0.7,
wbaer=", opt$S450[ln], ",
abaer=", opt$SAE[ln], ",
zout=0,100, !get output at the surface and top of atmosphere
iout=10, !output mode for TOA and surface spectral fluxes
nstr=4, !number of discrete ordinates in both hemispheres in DISORT
/
")
write.table(file = "/mnt/d/Model/sbdart/INPUT", descript, row.names = F, col.names = F, quote = FALSE, append = F)
system("/mnt/d/Model/sbdart/./sbdart.e >> HN_sbdart.out")
}
opt = read.csv(file = "DATA/OUT/hn_aod.csv")
opt$Time = as.POSIXct(opt$Time)
aero = read.table(file = "hn_direct.out")
clear = read.table(file = "hn_clear2.out")
Ftoadn_aer = aero[,4]
Ftoaup_aer = aero[,5]
Fbotdn_aer = aero[,7]
Fbotup_aer = aero[,8]
Ftoadn_clr = clear[,4]
Ftoaup_clr = clear[,5]
Fbotdn_clr = clear[,7]
Fbotup_clr = clear[,8]
Ftoant_aer = Ftoadn_aer - Ftoaup_aer # TOA net
Fbotnt_aer = Fbotdn_aer - Fbotup_aer # bottom net
Ftoant_clr = Ftoadn_clr - Ftoaup_clr # TOA net
Fbotnt_clr = Fbotdn_clr - Fbotup_clr # bottom net
opt$tfc = Ftoant_aer - Ftoant_clr #TOA forcing
opt$bfc = Fbotnt_aer - Fbotnt_clr#Bottom forcing
boxplot(opt$bfc ~ month(opt$Time))
boxplot(opt$tfc ~ month(opt$Time))
mean(opt$bfc)
mean(opt$tfc)
d.opt = aggregate(opt[,c("tfc", "bfc")], by = list(format(opt$Time, "%Y-%m-%d")), mean, na.rm = T)
d.opt$N = aggregate(opt[,c("tfc", "bfc")], by = list(format(opt$Time, "%Y-%m-%d")), length)[,2]
d.opt = d.opt[d.opt$N >= 3, ]
d.opt$Group.1 = as.POSIXct(d.opt$Group.1)
m.opt = aggregate(d.opt[,c("tfc", "bfc")], by = list(format(d.opt$Group.1, "%Y-%m-01")), mean, na.rm = T)
m.opt$N = aggregate(d.opt[,c("tfc", "bfc")], by = list(format(d.opt$Group.1, "%Y-%m-01")), length)[,2]
m.opt = m.opt[m.opt$N >= 5, ]
m.opt$Group.1 = as.POSIXct(m.opt$Group.1)
mm.opt = aggregate(m.opt[,c("tfc", "bfc")], by = list(format(m.opt$Group.1, "%m")), mean, na.rm = T)