setwd("~/Downloads/hurricane")
require(plyr)
## Loading required package: plyr
# http://vis.computer.org/vis2004contest/data.html
# http://www.vets.ucar.edu/vg/isabeldata/readme.html Dimensions: 500 x 500
# x 100 Physical Scale: 2139km (east-west) x 2004km (north-south) x 19.8km
# (vertical) Note: The east-west distance of 2139km is measured at the
# southern latitude; the distance is only 1741km at the northern latitude
# due to the curvature of the earth. Physical Location: Longitude (x):
# 83W to 62W Latitude (y): 23.7N to 41.7N Height (z): 0.035km to 19.835km
# Rescale x y z values to actual global coordinates
rescale <- function(x, newmin, newmax) {
(x - min(x))/diff(range(x)) * (newmax - newmin) + newmin
}
read_dat <- function(filename, varname) {
require(reshape2)
xres <- 500
yres <- 500
zres <- 100
# Map the 1:500 to the x coordinates, and we need to reverse for some
# reason
xcoord <- rev(rescale(1:xres, -62, -83))
ycoord <- rev(rescale(1:yres, 23.7, 41.7))
zcoord <- rescale(1:zres, 0.035, 19.835)
# xcoord <- xres:1 ycoord <- yres:1 zcoord <- 1:zres
fh <- file(filename, "rb")
dat <- array(readBin(fh, what = "double", n = xres * yres * zres, size = 4,
endian = "big"), dim = c(xres, yres, zres), dimnames = list(y = ycoord,
x = xcoord, z = zcoord))
close(fh)
# Take a subset, skipping every 5
xseq <- seq(1, xres, by = 5)
yseq <- seq(1, yres, by = 5)
zseq <- seq(1, zres, by = 10)
dat <- dat[yseq, xseq, zseq]
# Convert 1e35 values to NA, accounting for FP error
dat[abs(dat/1e+35 - 1) < 1e-06] <- NA
melt(dat, value.name = varname)
}
# Get each x y z data
vx <- read_dat("Uf40.bin", "vx")
## Loading required package: reshape2
vy <- read_dat("Vf40.bin", "vy")
vz <- read_dat("Wf40.bin", "vz")
t <- read_dat("TCf40.bin", "t")
# Merge them into a single data frame
isabel <- join(join(join(vx, vy), vz), t)
## Joining by: y, x, z
## Joining by: y, x, z
## Joining by: y, x, z
# Calculate speed
isabel <- transform(isabel, speed = sqrt(vx^2 + vy^2))
# Take a slice Keeps 1 out of every 'by' values in x every_n(3:10) # 3 5 7
# 9 every_n(3:10, 3) # 3 6 9
every_n <- function(x, by = 2) {
x <- sort(x)
x[seq(1, length(x), by = by)]
}
keepx <- every_n(unique(isabel$x), 4)
keepy <- every_n(unique(isabel$y), 4)
keepz <- every_n(unique(isabel$z), 2)
isabel1 <- subset(isabel, x %in% keepx & y %in% keepy & z %in% keepz)
library(grid)
library(ggplot2)
# Speed
ggplot(subset(isabel, z == min(z)), aes(x = x, y = y, fill = speed)) + geom_tile()
# Speed faceted
ggplot(isabel, aes(x = x, y = y, fill = speed)) + geom_raster() + facet_wrap(~z)
# Temp
ggplot(subset(isabel, z == min(z)), aes(x = x, y = y, fill = t)) + geom_tile()
# Temp faceted
ggplot(subset(isabel, z < 8), aes(x = x, y = y, fill = t)) + geom_raster() +
facet_wrap(~z)
# Vector field
ggplot(subset(isabel1, z == min(z)), aes(x = x, y = y)) + geom_segment(aes(xend = x +
vx/50, yend = y + vy/50), arrow = arrow(length = unit(0.1, "cm")), size = 0.25)
# Faceted vector field
ggplot(isabel1, aes(x = x, y = y)) + geom_segment(aes(xend = x + vx/40, yend = y +
vy/40), arrow = arrow(length = unit(0.075, "cm")), size = 0.25) + facet_wrap(~z)
# Map with speed
require(maps)
## Loading required package: maps
usa <- map_data("usa")
ggplot(subset(isabel, z == 0.035), aes(x = x, y = y)) + geom_raster(aes(fill = speed),
alpha = 0.5) + geom_path(aes(x = long, y = lat), data = usa) + scale_x_continuous(breaks = seq(-100,
60, by = 2)) + scale_y_continuous(breaks = seq(24, 50, by = 2))
ggplot(subset(isabel1, z == min(z)), aes(x = x, y = y)) + geom_segment(aes(xend = x +
vx/50, yend = y + vy/50), arrow = arrow(length = unit(0.1, "cm")), size = 0.25) +
geom_path(aes(x = long, y = lat), data = usa)
# ================ test
dat <- array(1:12, dim = c(2, 2, 3), dimnames = list(x = 1:2, y = 1:2, z = 1:3))
d2 <- melt(dat)
d2
## x y z value
## 1 1 1 1 1
## 2 2 1 1 2
## 3 1 2 1 3
## 4 2 2 1 4
## 5 1 1 2 5
## 6 2 1 2 6
## 7 1 2 2 7
## 8 2 2 2 8
## 9 1 1 3 9
## 10 2 1 3 10
## 11 1 2 3 11
## 12 2 2 3 12
ggplot(subset(d2, z == 2), aes(x = x, y = y, fill = value)) + geom_tile()
