Subset model data output

Description

Subset data structure returned by ReadGrib or DODSGrab by variables, levels, etc.

Usage

1
2
SubsetNOMADS(model.data, levels = NULL, variables = NULL, lon = NULL, 
lat = NULL, ensembles = NULL, forecast.date = NULL, model.run.date = NULL)

Arguments

model.data

Data structure from ReadGrib or DODSGrab

levels

Vector of levels to keep

variables

Vector of variables to keep

lon

Vector of longitudes of model nodes to keep

lat

Vector of latitudes of model nodes to keep

ensembles

Vector of ensemble runs to keep

forecast.date

Vector of forecast dates to keep

model.run.date

Vector of model run dates to keep

Value

model.data.sub

A subset of model.data.

Note

Multiple elements in each argument vector are obviously OR (i.e. variables are “tmpprs” OR “hgtprs”) but multiple subset vectors are AND. Thus it is simple to construct a model.data.sub with variables: tmpprs and hgptprs only from ensemble runs 3 and 4, for example.

Author(s)

Daniel C. Bowman daniel.bowman@unc.edu

See Also

ReadGrib, DODSGrab, ModelGrid, BuildProfile

Examples

 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
## Not run: 
#Plot winds from 20 GENS model runs

#Get the latest ensemble model run
model.urls <- GetDODSDates("gens")
latest.model <- tail(model.urls$url, 1)
model.runs <- GetDODSModelRuns(latest.model)
model.run <- tail(model.runs$model.run[grepl("all", 
    model.runs$model.run)], 1)

#Define region of interest: Chapel Hill, NC
lon <- -79.052104
lat <- 35.907553

lons <- seq(0, 359, by = 1)
lats <- seq(-90, 90, by = 1)

lon.diff <- abs(lon + 360 - lons)
lat.diff <- abs(lat - lats)

model.lon.ind <- which(lon.diff == min(lon.diff)) - 1
model.lat.ind <- which(lat.diff == min(lat.diff)) - 1

#Set up call to NOMADS
time <- c(0, 0)      #Analysis(?) model only
node.lon  <- c(model.lon.ind - 2, model.lon.ind + 2)  #Longitude grid
node.lat  <- c(model.lat.ind - 2, model.lat.ind + 2)   #Latitude grid
variables <- c("ugrdprs", "vgrdprs", "hgtprs") #Wind speeds, and geopotential height
levels    <- c(0, 25) #All available levels
ensembles <- c(0, 20)  #All available ensembles

model.data <- DODSGrab(latest.model, model.run, variables, time, 
    node.lon, node.lat, levels = levels, ensembles = ensembles)

#Plot winds
zonal.wind <- NULL
merid.wind <- NULL
height     <- NULL

for(k in ((ensembles[1]:ensembles[2] + 1))) {
    model.data.sub <- SubsetNOMADS(model.data, ensembles = c(k), 
    variables = c("hgtprs", "ugrdprs", "vgrdprs"))
    profile <- BuildProfile(model.data.sub, lon, lat)
    hgt     <- profile[[1]]$profile.data[, which(profile[[1]]$variables == "hgtprs"),]
    ugrd    <- profile[[1]]$profile.data[, which(profile[[1]]$variables == "ugrdprs"),]
    vgrd    <- profile[[1]]$profile.data[, which(profile[[1]]$variables == "vgrdprs"),]

   synth.hgt <- seq(min(hgt),
       max(hgt), length.out = 1000)
   ugrd.spline <- splinefun(hgt, ugrd, method = "natural")
   vgrd.spline <- splinefun(hgt, vgrd, method = "natural")
   zonal.wind[[k]] <- ugrd.spline(synth.hgt)
   merid.wind[[k]] <- vgrd.spline(synth.hgt)
   height[[k]] <- synth.hgt
}

PlotWindProfile(zonal.wind, merid.wind, height, lines = TRUE,
    points = FALSE, elev.circles = c(0, 15000, 30000), elev.labels = c(0, 15, 30),
    radial.lines = seq(45, 360, by = 45), colorbar = TRUE, invert = FALSE,
    point.cex = 2, pch = 19, lty = 1, lwd = 1,
    height.range = c(0, 30000), colorbar.label = "Wind Speed (m/s)")

## End(Not run)