/
tc_TC_wnds_stns_daily.r
134 lines (124 loc) · 5.34 KB
/
tc_TC_wnds_stns_daily.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
#Title: Program 'tc_TC_wnds_stns.r'
#Author: Sanabria A., augusto.sanabria@ga.gov.au
#CreationDate: 2013-10-04
#Description: Program to extract cyclonic winds from BoM stations
#located in the cyclonic region of Australia.
#It reads a dataset of cyclonic tracks affecting a BoM recording station
#and extracts corresponding values from the station wind speeds.
#Dataset includes max wind speed of each cyclone coming closer than 600 km to
#the station.
#Output is written to external file using station ID as identifier.
#Keywords: BoM recording stations, BoM database of TC tracks, observed cyclonic winds
#Reference:
#SeeAlso:
#Version: 1.1
#Modified by: AS
#ModifiedDate: 2014-01-10
#Modification: The cyclonic wind dataset now includes only max wind speeds per cyclone.
#Required:
#sdir <<- "/home/asanabri/R/my_rsoftw/"
indir <- "C:/WorkSpace/tcobs/"
setwd(indir)
#source(paste(sdir,"get_TC_gusts.r",sep="") )
#source(paste(sdir,"get_TC_gusts_per_cyclone.r",sep="") )
source('get_TC_gusts_per_cyclone.r')
source('loadData.r')
#Returns:
#void, produces external file.
#
t1 <- proc.time()
#Read database of stations affected by TCs:
outdir <- "C:/Workspace/data/derived/tc/"
datadir <- "C:/Workspace/data/raw/daily_max_wind_gust/"
fname <- paste(outdir, "stn_TC_dist.allstns_200.txt", sep = "")
bom_stns_wTC <- read.table(fname, sep = ",", skip = 1, header = F)
#Read BoM datasets (half-hour) located in Dir 'datadir'
setwd(datadir)
all_bomd <- dir(pattern = "DC02D_Data")
#List stns not in bom data:
outf001 <- "stns_not_in_bomd.txt"
write.table("List of stations not in BoM-provided datasets",
file = paste(outdir, outf001, sep = ""), append = FALSE)
write.table("stn_ID stn_name stn_ST stn_lat stn_lon",
file = paste(outdir, outf001, sep = ""),
append = TRUE)
#Transform BoM stn_ID in stn name into numeric:
bom_substr <- c()
for (i in 1:length(all_bomd)) {
bom_substr[i] <- as.numeric(substring(all_bomd[i], 12, 17))
}
#Split the data into stations, extract wind gusts at the time of the TC track approaching the station and
#create a dataset of TC winds at that station:
grouped_data <- split(bom_stns_wTC, bom_stns_wTC$V2)
stn_char <- c()
cnt_stns <- 0
all_TC_gusts <- c()
for (i in 1:length(grouped_data)) {
# extract wind gust from BoM dataset
stn_id <- grouped_data[[i]][2]
stn_ST <- unique(grouped_data[[i]][3])
stn_name <- unique(grouped_data[[i]][1])
cycl_name <- grouped_data[[i]][6]
stn_lat <- grouped_data[[i]][4]
stn_lon <- grouped_data[[i]][5]
stn_char[i] <- paste(i, unique(stn_id),
stn_name$V1,stn_ST$V3,
" (",unique(stn_lat), unique(stn_lon),
")", sep = ",")
# open corresponding BoM dataset:
bdatas <- which(as.numeric(unique(stn_id)) == as.numeric(bom_substr))
# Test for errors (n.b. character(0) is numeric!!)
ZZ <- try(ifelse(bdatas >= 1, bdatas, NA))
# Write to external file list of stations not found in BoM-provided datasets:
if (!is.numeric(ZZ)) {
write.table(stn_char[i],
file = paste(outdir, outf001, sep = ""),
append = TRUE, row.names = FALSE,
col.names = FALSE, quote = FALSE)
next
}
# Write all stns found in BoM-provided datasets to R log file (to plot them in the OZ map):
str1 <- paste(unique(stn_id), ", ", stn_name$V1, ",",
unique(stn_lat), ",", unique(stn_lon),
",above, y" ,sep = "")
print(str1)
# TC_gusts <- get_TC_gusts(all_bomd[bdatas],unique(stn_ST),grouped_data[[i]][8],grouped_data[[i]][9])
#print("Split data into cyclones in order to calc. max speed per cyclone")
cycl_wsp <- cbind(grouped_data[[i]][8], grouped_data[[i]][9], grouped_data[[i]][6])
cycl_mx_wsp <- split(cycl_wsp, cycl_wsp$V6, drop = TRUE)
# loop tr' all cyclones:
for (j in 1:length(cycl_mx_wsp)) {
print(paste("Extracting TC gusts for ",as.character(cycl_mx_wsp[[j]][3][1,1])))
obs_gusts <- readDailyObs(all_bomd[bdatas], units = "m/s")
TC_gusts <- get_TC_gusts_per_cyclone(obs_gusts, unique(stn_ST),
cycl_mx_wsp[[j]][1], cycl_mx_wsp[[j]][2],
cycl_mx_wsp[[j]][3])
if (all(is.na(TC_gusts))) next
# Keep only max wind speed of this cyclone:
#print("Filter maximum gusts")
gusts = as.numeric(TC_gusts$gust)
mx_gwsp <- which(gusts == max(gusts, na.rm = TRUE))
if (length(mx_gwsp) > 1) {
mx_gwsp = mx_gwsp[1]
}
#print("Appending maximum gust to list of all records")
all_TC_gusts <- try(rbind(all_TC_gusts, TC_gusts[mx_gwsp,]))
# all_TC_gusts <- try(rbind(all_TC_gusts,TC_gusts ) ) #write all wnd speeds not only max
}
# Write to external file winds in BoM-provided datasets which match date/time in cyclone track,
# one dataset per station ID:
cnt_stns <- cnt_stns + 1
# write name with 6 digits (as in BoM station name convention):
new_str <- formatC(unique(grouped_data[[i]]$V2), width = 6, flag = "0")
out_main <- paste("bom_", new_str, ".csv",sep = "")
print(paste("Writing data to: ",out_main))
if (length(all_TC_gusts) > 0) {
write.csv(all_TC_gusts, file = paste(outdir, out_main, sep = "") )
all_TC_gusts <- c()
}
}
etime <- proc.time() - t1
print(paste("Stations in BoM-provided datasets hit by TC = ", cnt_stns, sep = "") )
print(paste("Proc. time = ", etime, sep = ""))
#
print("Normal program termination")