This package provides tools for interfacing with input and output files used by the Stochastic Time Inverted Lagrangian Transport (STILT) model.
If R is using the intel fortran compiler, additional compilation flags must be included. Run the following command to check if the ifort compiler is in use, and add PKG_FCFLAGS = -assume byterecl to ~/.R/Makevars if necessary.
R CMD config FC | grep -q ifort && mkdir -p ~/.R && (echo "PKG_FCFLAGS = -assume byterecl" >> ~/.R/Makevars && echo "Modified compiler options at ~/.R/Makevars")
R CMD config FC | grep -q gfortran && mkdir -p ~/.R && (echo "PKG_FCFLAGS = -frecord-marker=8" >> ~/.R/Makevars && echo "Modified compiler options at ~/.R/Makevars")Package can be installed using devtools.
if (!require('devtools')) install.packages('devtools')
devtools::install_github('uataq/stiltread')NOTE: the rgdal package was removed from the CRAN repository as of October 2023.
To install rgdal nonetheless, it can be installed from source:
url <- "https://cran.r-project.org/src/contrib/Archive/rgdal/rgdal_1.6-7.tar.gz"
install.packages(url, type="source", repos=NULL)Then, load the package and install the dependencies.
library(stiltread)
install_stiltread()To read a raster layer from a HRRR data file,
library(stiltread)
shgt <- read_met('/path/to/arl/20150617.00z.hrrra',
var = 'shgt',
yy = 15,
mm = 6,
dd = 17,
hh = 3,
lvl = 0)
shgt
# class : RasterLayer
# dimensions : 1059, 1799, 1905141 (nrow, ncol, ncell)
# resolution : 3005.169, 2991.952 (x, y)
# extent : -2700991, 2705308, -1580359, 1588118 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=lcc +lon_0=-97.5 +lat_0=38.5 +lat_1=38.5 +lat_2=38.5 +ellps=WGS84 +no_defs
# data source : in memory
# names : layer
# values : -80.036, 4175.964 (min, max)which returns a raster::RasterLayer containing data, grid coordinates, and projection information using the standard PROJ4 format.
Since wind is grid relative, use read_met_wind instead of read_met to access wind fields. This is called in the same way as read_met, only no var argument is required. This returns a rasterStack with u and v layers that have been rotated into +proj=longlat.
library(stiltread)
uv <- read_met_wind('/path/to/arl/20150617.00z.hrrra',
yy = 15,
mm = 6,
dd = 17,
hh = 3,
lvl = 0)
uv
# class : RasterStack
# dimensions : 72, 78, 5616, 2 (nrow, ncol, ncell, nlayers)
# resolution : 0.0345, 0.0261 (x, y)
# extent : -113.4537, -110.7627, 39.81678, 41.69598 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84
# names : u, v
# min values : -6.965451, -10.460653
# max values : 5.198425, 5.289214The u and v layers can be accessed by name using
uv$u
# class : RasterLayer
# dimensions : 72, 78, 5616 (nrow, ncol, ncell)
# resolution : 0.0345, 0.0261 (x, y)
# extent : -113.4537, -110.7627, 39.81678, 41.69598 (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=longlat +ellps=WGS84
# data source : in memory
# names : u
# values : -6.965451, 5.198425 (min, max)Since read_met returns a raster::RasterLayer with projection metadata, we can leverage standard mapping libraries in R.
plot(shgt, main = 'Elevation')
For an interactive leaflet map,
library(leaflet)
leaflet() %>%
addProviderTiles('CartoDB.Positron') %>%
addRasterImage(shgt, opacity = 0.5, colors = viridis::viridis(16))
- Returned raster layers are limited to files that use Lambert Conformal Conic for the map projection, which is typical for midlatitudes.
- Upper levels of pressure variables are often recorded as a difference from the base level. This calculation needs to be performed by the user.