-
Notifications
You must be signed in to change notification settings - Fork 0
/
tc_match_stns.r
118 lines (106 loc) · 5.1 KB
/
tc_match_stns.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
#Title: Program 'tc_match_stns.r'
#Author: Sanabria A., augusto.sanabria@ga.gov.au
#CreationDate: 2013-09-19
#Description: Program to match BoM stations with TC tracks.
#The station characteristics are read from BoM station database
#TC tracks are also provided by BoM.
#This script compares station location with TC tracks, if the TC track comes
#closer than 'distance' km to the station and the station was in operation at the time,
#write the station name and location; TC distance to the station, date and TC radius
#(MN_RADIUS_SF_WND) into an external file ("stn_TC_dist.txt").
#Keywords: BoM recording stations, BoM database of TC tracks
#Reference:
#SeeAlso:
#Version: 1.0
#Modified by:
#ModifiedDate:
#Modification:
#Required:
library(logging)
basicConfig()
addHandler(writeToFile, logger="", file = 'tc_match_stns.log')
loginfo(paste("Running ",sys.frame(1)$ofile, sep=''), logger='')
#Returns:
#void, produces external file.
#
#
#Read TC tracks:
source('loadData.r')
tc_fname <- "C:/WorkSpace/data/TC/Allstorms.ibtracs_wmo.v03r09.csv"
loginfo(paste("Source TC file: ", tc_fname, sep=""), logger='')
TC_datab <- read.csv(tc_fname, skip = 0, header = T)
OZ_TC_datab <- with(TC_datab,data.frame(serialno = Serial_Num, name = Name, UTC_time = ISO_time,
lon = Longitude, lat = Latitude,
pressure = Pres.WMO.))
#Read list of BoM stns and location:
datadir <- "C:/Workspace/data/raw/daily_max_wind_gust/"
stnFile <- paste(datadir, "DC02D_StnDet_999999999425050.txt", sep="")
stations <- readStationFile(stnFile, header = F, skip = 1)
loginfo(paste("Station file: ", stnFile, sep=""), logger="")
#Write all stations which are closer than 'distance' km to a TC track and are open at the time of the TC pass
#Write also stations not included in the two rules listed above with the reason (as a check)
t1 <- proc.time()
distance <- 200.0
keep_name <- c()
outdir <- "C:/Workspace/data/derived/tcobs/daily/"
outputFile <- paste(outdir, "stn_TC_dist.allstns_200.txt", sep = "")
loginfo(paste("Output will be stored in ", outputFile, sep= ''), logger='')
write('"Station", "BoM_ID", "State", "Lat", "Lon", "Name", "Serialnum", "dist (km)", "date passed", "TC pressure (hPa)"', file = outputFile, append = FALSE)
for (i in 1:length(stations$stnId)) {
stnID <- stations$stnId[i]
stn_name <- stations$stnName[i]
stn_state <- stations$stnState[i]
loginfo(paste("Station: ",stn_name, " (", stnID, ")",sep = ''), logger='')
stn_lat <- stations$stnLat[i]*(pi/180.0)
stn_lon <- stations$stnLon[i]*(pi/180.0)
nevents <- 0
stn_start <- stations$stnDataStart[i]
if (is.na(stations$stnDataEnd[i])) {
stn_close = strptime("2100-01-01", format = "%Y-%m-%d", tz = "GMT")
}else{
stn_close = stations$stnDataEnd[i]
}
for (k in 1:length(OZ_TC_datab$serialno)) {
lat2 <- OZ_TC_datab[k, 5]*(pi/180.0)
lon2 <- OZ_TC_datab[k, 4]*(pi/180.0)
if (is.na(lat2) | is.na(lon2) ) {
#print("**Warning: empty line (ignored)")
#print(paste("TC_name(previous) = ", OZ_TC_datab[k - 1, 1], sep = ""))
next
}
# Calculate distance of TC location from the station:
dd <- 2*asin(sqrt((sin((stn_lat - lat2)/2)) ^ 2 +
cos(stn_lat)*cos(lat2)*(sin((stn_lon - lon2)/2)) ^ 2))
dist_2stn <- 1.852*(180*60/pi)*dd
#print(paste("dist = ",dist_2stn,sep="") )
if (dist_2stn <= distance) {
keep_name <- as.character(OZ_TC_datab$name[k])
serialnum <- as.character(OZ_TC_datab$serialno[k])
TC_date = strptime(OZ_TC_datab[k, 3], format = "%Y-%m-%d %H:%M:%S", tz = "GMT")
if (TC_date >= stn_start & TC_date <= stn_close) {
#Station is close enough, record station, distance and time of passage:
write.table(paste(as.character(stn_name),
stnID, stn_state,
signif((180.0/pi)*stn_lat, 8),
signif((180.0/pi)*stn_lon, 8),
keep_name, serialnum, round(dist_2stn, 1),
OZ_TC_datab[k, 3], OZ_TC_datab[k, 6], sep = ","),
file = outputFile, append = TRUE, row.names = FALSE,
col.names = FALSE, quote = FALSE)
nevents = nevents + 1
}else{
#Station was NOT open at time of TC hit:
#write.table(paste(as.character(stn_name)," (",signif((180.0/pi)*stn_lat,8),",",signif((180.0/pi)*stn_lon,8),") ",keep_name," ",round(dist_2stn,1)," ",OZ_TC_datab[k,2]," No time match!",sep=""),
#file=outf00,append=TRUE,row.names=FALSE,col.names=FALSE,quote=FALSE)
}
}else{
#Station is farther than 'distance' km, write it in to 'No_TC_hit' file:
#write.table(paste(stn_name," (",signif((180.0/pi)*stn_lat,8),",",signif((180.0/pi)*stn_lon,8),") ",keep_name," ",round(dist_2stn,1)," ",OZ_TC_datab[k,2]," Too far!",sep=""),
#file=outf00,append=TRUE,row.names=FALSE,col.names=FALSE)
}
}
loginfo(paste(nevents, "passed within ", distance, " km of ", stn_name, sep=''), logger = '')
}
etime <- proc.time() - t1
loginfo(paste("Proc. time = ", etime, sep = ""), logger='')
loginfo(paste("Finished ", sys.frame(1)$ofile, sep=''), logger='')