-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcyclone_paral.R
106 lines (90 loc) · 4.75 KB
/
cyclone_paral.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
#https://www.ncdc.noaa.gov/ibtracs/index.php?name=numbering
#https://www.ncdc.noaa.gov/ibtracs/index.php?name=wmo-var
#Time zone: UTC
rm(list = ls())
#Prepare ----
library(data.table)
cycl = fread("Allstorms.ibtracs_wmo.v03r10.csv", skip = 1)
#separate the row that describes unit
unit = cycl[1, ]
cycl = cycl[-1, ]
#convert columns to correct format
cycl[, c(9:12, 14:15)] = lapply(cycl[, c(9:12, 14:15)], as.numeric)
cycl$ISO_time = as.POSIXct(cycl$ISO_time, tz = "GMT")
cycl$Season_Name = paste0(cycl$Season, cycl$Name)
#Retrieve ----
library(geosphere, lib.loc = "R_library")
library(sp)
library(raster)
RetrCycl = function(year, x0, x1, y0, y1, resol, thres = 150, mapout = T, mapname = "map"){
cycl_filt = cycl[Season %in% year, ] #filter by year
#set up global coordinates
lon = rep(seq(x0 + resol / 2, x1 - resol / 2, by = resol), (y1 - y0) / resol)
lat = rep(seq(y1 - resol / 2, y0 + resol / 2, by = -resol), each = (x1 - x0) / resol)
coord = SpatialPointsDataFrame(matrix(c(lon, lat), ncol = 2), data.frame(ID = seq(1:length(lon))),
proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
#set up coordinates of typhoon events
coord_cycl = SpatialPointsDataFrame(as.matrix(cycl_filt[, list(Longitude, Latitude)]), data.frame(pt = seq(1:nrow(cycl_filt))),
proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
mat = distm(coord_cycl, coord) / 1000 #calculate distance
lst = split(t(mat), seq(nrow(t(mat))))
cycl_sel = lapply(lst, function(x) cycl_filt[which(x < thres), ]) #filter by distance
frequ = lapply(cycl_sel, function(x) length(unique(x$Season_Name)) / length(year)) #frequency (number of typhoons)
dur = lapply(cycl_sel, function(x) nrow(x) / length(year)) #duration (number of records)
wind = lapply(cycl_sel, function(x) mean(tapply(x$`Wind(WMO)`, x$Season_Name, max))) #maximum wind speed
press = lapply(cycl_sel, function(x) mean(tapply(x$`Pres(WMO)`, x$Season_Name, min))) #minimum atmospheric pressure
#store in rasters
ras.freq = raster(xmn = x0, xmx = x1, ymn = y0, ymx = y1, resolution = resol)
ras.dur = raster(xmn = x0, xmx = x1, ymn = y0, ymx = y1, resolution = resol)
ras.wind = raster(xmn = x0, xmx = x1, ymn = y0, ymx = y1, resolution = resol)
ras.press = raster(xmn = x0, xmx = x1, ymn = y0, ymx = y1, resolution = resol)
values(ras.freq) = unlist(frequ)
values(ras.dur) = unlist(dur)
values(ras.wind) = unlist(wind)
values(ras.press) = unlist(press)
dirct = "res_cycl/"
writeRaster(ras.freq, paste0(dirct, mapname, "_freq.grd"), format = "raster")
writeRaster(ras.dur, paste0(dirct, mapname, "_dur.grd"), format = "raster")
writeRaster(ras.wind, paste0(dirct, mapname, "_wind.grd"), format = "raster")
writeRaster(ras.press, paste0(dirct, mapname, "_press.grd"), format = "raster")
if(mapout){
library(maps)
#plot frequency
pdf(paste0(dirct, mapname, "_freq.pdf"))
plot(ras.freq, col = rev(heat.colors(300)))
map(database = "world", xlim = c(x0, x1), ylim = c(y0, y1), add = T)
abline(h = 0, v = 0) #equator and prime meridian
abline(h = c(-23.5, 23.5), lty = 2) #tropics
abline(h = 24, v = 121, lty = 3) #fushan
dev.off()
#plot duration
pdf(paste0(dirct, mapname, "_dur.pdf"))
plot(ras.dur, col = rev(heat.colors(300)))
map(database = "world", xlim = c(x0, x1), ylim = c(y0, y1), add = T)
abline(h = 0, v = 0) #equator and prime meridian
abline(h = c(-23.5, 23.5), lty = 2) #tropics
abline(h = 24, v = 121, lty = 3) #fushan
dev.off()
#plot wind speed
pdf(paste0(dirct, mapname, "_wind.pdf"))
plot(ras.wind, col = rev(heat.colors(300)))
map(database = "world", xlim = c(x0, x1), ylim = c(y0, y1), add = T)
abline(h = 0, v = 0) #equator and prime meridian
abline(h = c(-23.5, 23.5), lty = 2) #tropics
abline(h = 24, v = 121, lty = 3) #fushan
dev.off()
#plot atmospheric pressure
pdf(paste0(dirct, mapname, "_press.pdf"))
arg = list(at = seq(800, 1100, by = 50), labels = seq(800, 1100, by = 50))
plot(ras.press, col = heat.colors(300), breaks = 800:1100, axis.args = arg, zlim = c(800, 1100))
map(database = "world", xlim = c(x0, x1), ylim = c(y0, y1), add = T)
abline(h = 0, v = 0) #equator and prime meridian
abline(h = c(-23.5, 23.5), lty = 2) #tropics
abline(h = 24, v = 121, lty = 3) #fushan
dev.off()
}
}
#RetrCycl(year = 2001:2003, x0 = 0, x1 = 180, y0 = 0, y1 = 90, resol = 1, mapname = "NE_coarse")
RetrCycl(year = 2001:2003, x0 = -180, x1 = 180, y0 = -90, y1 = 90, resol = 1, mapname = "world_coarse")
RetrCycl(year = 2001:2003, x0 = 0, x1 = 180, y0 = 0, y1 = 90, resol = 0.5, mapname = "NE")
RetrCycl(year = 2001:2003, x0 = -180, x1 = 180, y0 = -90, y1 = 90, resol = 0.5, mapname = "world")