Contour Plot of pollutant concentration field

Description

The function contourPlot2 generates a contour plot of a given quantity, such as the ground concentration of an airborne pollutant or odour, defined on a regular grid.

Usage

contourPlot2(
  data,
  x = "x",
  y = "y",
  z = "z",
  domain = NULL,
  background = NULL,
  underlayer = NULL,
  overlayer = NULL,
  legend = NULL,
  levels = NULL,
  size = 0,
  fill = TRUE,
  tile = FALSE,
  transparency = 0.75,
  colors = NULL,
  mask = NULL,
  inverse_mask = FALSE,
  bare = FALSE
)

Arguments

Details

This is a convenience function to plot contour levels of a scalar quantity such as pollutants computed by a dispersion model, with ggplot2 version >= 3.3.0.

Data are required to be on a regular grid, typically (but not necessarily) in UTM coordinates. The input dataframe has to be in long format, i.e. one line per value to be plotted. The names of the columns corresponding to x, y and z can be specified in the input parameters.

If tile = TRUE data are shown as they are, without any graphical interpolation required for contour plots. This is helpful when you want to visualise the raw data. Since version 2.4.0, when tile = TRUE the intervals include the lowest bound and exclude the highest bound: ⁠[min, max)⁠. Note: In previous version it was the opposite.

underlayer and overlayer layers are ggplot2 objects to be shown at different levels of the vertical stack of the plot. These are useful to show topographical information related to the plot, such as sources or receptors locations.

When a shp file is given to the mask argument the plot is drawn only inside the polygon. In order to avoid boundary artifacts due to reduced resolution, original data are resampled to higher resolution (currently set to 10 times the original one.) Ifinverse_mask is set to TRUE, the plot is drawn outside the polygon. The mask feature is based on the terra::mask() function. The CRS of the shp file is applied to the data in the data.frame. Please, keep in mind this feature is still experimental.

Value

A ggplot2 object.

Examples

# Load example data in long format
data(volcano)
volcano3d <- reshape2::melt(volcano)
names(volcano3d) <- c("x", "y", "z")
# Contour plot with default options
v <- contourPlot2(volcano3d)
v

# Set levels, and properly format the legend title:
contourPlot2(
    volcano3d,
    levels = c(-Inf, seq(100, 200, 20), Inf),
    legend = expression("PM"[10] ~ "[" * mu * "g m"^-3 * "]")
)

# Sometimes, instead of a contour plot it is better to plot the original
# raster data, without any interpolation:
contourPlot2(
    volcano3d,
    levels = c(-Inf, seq(100, 200, 20), Inf),
    tile = TRUE
)

# Since contourPlot2 returns a `ggplot2` object, you can add instructions as:
library(ggplot2)
v +
    ggtitle("Example volcano data") +
    labs(x = NULL, y = NULL)

Download basemap from Italian National Geoportal

Description

Download the aerial orthophoto of the requested domain from the Italian National Geoportal.

Usage

downloadBasemap(
  file = file,
  xSW = NA,
  ySW = NA,
  xExt = NA,
  yExt = NA,
  crs = 32,
  width = 1024,
  height = 1024,
  units = "px",
  res = 72
)

Arguments

Details

The domain is specified by the South-West point coordinates, and its extension in the x and y directions. The Coordinate Reference System (CRS) is in UTM 32 or 33.

Note that, even if the downloading is successful the file might be empty due to some weird behaviour of the remote server from the PCN.

Value

The output is a tiff encoded with GeoTIFF metadata at the path provided. No value is returned.

Examples

## Not run: 
# Download a basemap of a domain with SW coordinates (410000, 5000500)
# in the UTM32 CRS and extension 5000m in both directions.

downloadBasemap(
    file = "./basemap.tif",
    xSW = 410000, ySW = 5000500, xExt = 5000, yExt = 5000
)

# Download a basemap of a domain with SW coordinates (410000, 5000500)
# in the UTM32 CRS and extension 5000m in both directions.
# The file has to be 2048 x 2048 pixels.

downloadBasemap(
    file = "./basemap.tif",
    xSW = 410000, ySW = 5000500, xExt = 5000, yExt = 5000,
    width = 2048, height = 2048
)

# Download a basemap of a domain with SW coordinates (410000, 5000500)
# in the UTM32 CRS and extension 5000m in both directions.
# The file has to be 10cm x 10cm with a resolution of 150 dpi.

downloadBasemap(
    file = "./basemap.tif",
    xSW = 410000, ySW = 5000500, xExt = 5000, yExt = 5000,
    width = 10, height = 10, units = "cm", res = 144
)

## End(Not run)

ADSO/BIN data import function

Description

Import data from ADSO/BIN binary file. It requires an active Python installation with the arinfopy library.

Usage

importADSOBIN(
  file = file.choose(),
  variable = NULL,
  slice = 1,
  deadline = 1,
  k = 1,
  kz = 1,
  dx = 0,
  dy = 0,
  destaggering = FALSE,
  raster.object = FALSE,
  verbose = FALSE
)

Arguments

Details

The importADSIOBIN() function imports data from an ADSO/BIN binary file. It relies on the ‘arinfopy’ (version >= 2.2.0) python library. For more information on the library see the GitHub repository.

For more information on the active python installation, check the documentation of reticulate.

Value

In standard use, importADSOBIN() return a data frame with ⁠(X, Y, Z)⁠ columns. Column Z contains the values of the requested variable. If the raster.object option is set, it returns a RasterLayer object.

See Also

importRaster(), importSurferGrd()

Examples

## Not run: 
# Read ground level (slice = 1) value of variable M001S001.
pm10 <- importADSOBIN(
    file = "average_2018.bin",
    variable = "M001S001",
    slice = 1
)

# Read deadline 12 of the second vertical level of temperature:
temperature <- importADSOBIN(
    file = "swift_surfpro_01-10_01_2018",
    variable = "TEMPK",
    slice = 2,
    deadline = 12
)

# Read varibale M001S001 at ground level, at given date and time,
# and print basic information:
nox <- importADSOBIN(
    file = "conc_01-10_07_2018",
    variable = "M001S001",
    slice = 1,
    deadline = "2018/07/02 12:00",
    verbose = TRUE
)

## End(Not run)

Import generic raster file

Description

A function to import the first layer of a generic raster file.

Usage

importRaster(
  file = file.choose(),
  k = 1,
  kz = 1,
  dx = 0,
  dy = 0,
  destaggering = FALSE,
  variable = NULL,
  verbose = FALSE
)

Arguments

Details

This function is based on the terra package and it can import any format managed by it as NetCDF.

Destaggering applies a shift equal to half grid size in both horizontal directions. It is useful for importing data from the SPRAY air quality dispoersion model and it is not applied by default.

An optional summary output can be printed out by setting the verbose parameter to TRUE.

Value

A data.frame with x, y and z columns for the grid cells coordiantes and the variable value.

See Also

importADSOBIN(), importSurferGrd()

Examples

## Not run: 
# Import binary (netcdf) file and convert coordinates from km to m,
# without destaggering. Variable name is "NOx".
mydata <- importRaster(
    file = "/path_to_file/filename.nc",
    variable = "NOx",
    k = 1000,
    destaggering = FALSE
)

# Import binary (netcdf) file and convert coordinates from km to m,
# with shift of 100 m in both directions:
mydata <- importRaster(
    file = "/path_to_file/filename.nc",
    variable = "pm10",
    k = 1000,
    dx = 100,
    dy = 100
)

## End(Not run)

Import Grid file

Description

A function to import data from Surfer text grid file.

Usage

importSurferGrd(fname, k = 1000, destaggering = FALSE)

Arguments

Details

Surfer grd file is imported and an array of x, y, z columns is returned X and y coordinates can be converted from km to m (default k=1000) and vice versa. Destaggering is applied by default.

Value

A dataset with x, y and z columns is returned.

See Also

importRaster(), importADSOBIN()

Examples

## Not run: 
# Import Surfer Grd file and convert coordinates from km to m,
# with destaggering
mydata <- importSurferGrd("/path_to_file/filename.grd", k = 1000)

# Import Surfer Grd file and do not convert coordinates, without destaggering
mydata <- importSurferGrd(
    "path_to_file/filename.grd",
    k = 1,
    destaggering = FALSE
)

## End(Not run)

Plot hourly average radiation

Description

Plot a histogram with hourly average of solar radiation, together with hourly maxima for June and December.

Usage

plotAvgRad(
  mydata,
  date = "date",
  rad = "radg",
  ylabel = NULL,
  title = "",
  locale = NULL
)

Arguments

Value

A ggplot2 plot.

See Also

plotStabilityClass(), plotAvgTemp()

Examples

data(stMeteo)
plotAvgRad(stMeteo, date = "date", rad = "radg")

Plot average temperature

Description

plotAvgTemp builds a bar plot of time average temperature and two line plots with maximum and minimum temperature.

Usage

plotAvgTemp(
  mydata,
  date = "date",
  temp = "temp",
  avg.time = "1 month",
  ylabel = NULL,
  title = "",
  locale = NULL
)

Arguments

Value

A plot with average, min and max temperature in a given range of time.

See Also

plotStabilityClass(), plotAvgRad()

Examples

# Plot average monthly temperature and curves with monthly maximum and minimum
data(stMeteo)
str(stMeteo)
plotAvgTemp(stMeteo)
# Add a custom title
plotAvgTemp(stMeteo, title = "Monthly temperature")

# Override default locale
plotAvgTemp(stMeteo, avg.time = "1 month", locale = "it_IT")

Plot stability class

Description

Histogram plot of stability classes by season or hour.

Usage

plotStabilityClass(
  mydata,
  date = "date",
  sc = "sc",
  type = "season",
  locale = NULL
)

Arguments

Details

Numerical values of stability classes are mapped as: 1 = A, 2 = B, ..., 6 = F.

Value

A ggplot2 plot.

See Also

stabilityClass(), plotAvgRad(), plotAvgTemp()

Examples

data("stMeteo")

# Season plot of stability class pgt
plotStabilityClass(stMeteo, date = "date", sc = "pgt", type = "season")

# Hourly plot of stability class pgt
plotStabilityClass(stMeteo, date = "date", sc = "pgt", type = "hour")

# Override default locale
plotStabilityClass(
    stMeteo,
    date = "date",
    sc = "pgt",
    type = "season",
    locale = "it_IT"
)

Remove data outliers

Description

Remove data outliers based on the interquartile range.

Usage

removeOutliers(x, k = 1.5)

Arguments

Details

The interquartile range IQR is computed from input dataset as IQR = Q3 - Q1, where Q1 is 25th percentile and Q3 is the 75th percentile. Values larger than Q3 + k * IQR and smaller than Q1 - k * IQR are deemed as outliers and substituted with NA's.

The default value of k is 1.5.

Value

A numeric vector with the same length as input vector.

Examples

mydata <- c(-10 * runif(10), runif(10))
removeOutliers(mydata)

Compute Rolling Max

Description

The function computes the rolling maximum value along a time series.

Usage

rollingMax(mydata, length = 24)

Arguments

Details

It calculates the maximum over consecutive elements centered within a specified window.

For each index i, it considers a window of length points centered around i. When length is odd, the center falls exactly on i and the window extends equally to both sides. When length is even, the window extends one less point to the left than to the right and the rolling max is not exactly centered.

Values near the start of the series use windows with fewer than length data points if there are not enough preceding elements to form a full window. Similarly for values at the end.

Value

A numeric vector containing rolling maximum values, with same dimensions as mydata.

Examples

# Compute rolling max over a 24-hour period on hourly time series data
data(stMeteo)
ws_24h <- rollingMax(mydata = stMeteo$ws, length = 24)

Meteorological dataset with hourly values

Description

A dataset containing 8760 hourly values of some meteorological variables corresponding to a full solar year.

Usage

stMeteo

Format

A data frame with 8760 rows and 7 variables:

date

date time in yyyy-mm-hh HH:MM:SS

ws

wind speed in m/s

wd

wind direction in deg.

temp

air temperature in C

radg

Global solar radiation in W/m^2

tcc

Total cloud cover in integers ranging from 0 to 8

pgt

Pasquill-Gifford-Turner stability class

Source

Self derived dataset.

Stability class.

Description

Computes stability class given net radiation, total cloud cover and wind speed.

Usage

stabilityClass(rad, tcc, ws, option = "iaea")

Arguments

Details

stabilityClass() computes stability class according to IAEA method based on net radiation, total cloud cover tcc and wind speed. Net radiation and wind are used by day; tcc and wind are used by night.

Three different alogorithms are implemented, selected by the option argument.

iaea option implements the *radiation-wind method recommended by the International Atomic Energy Agency (IAEA) and it is based on the net radiation during the day and cloud cover by night.

pasquill option is based on the original Pasquill formulation and lacks the "very weak" solar insolation present in the modified iaea version.

Eventually, the custom options is similar to iaea, with slightly different set of parameters for net radiation, wind speed and cloud cover.

Previously used option impact is the same as iaea and it is now deprecated.

Value

stabilityClass returns a numeric vector with Pasquill stability classes coded as: A = 1, B = 2, ... , F = 6 ranging from "very unstable" to "very stable".

See Also

plotStabilityClass()

Examples

# Compute stability class with custom algorithm
stMeteo$cst <- stabilityClass(
    rad = stMeteo$rad,
    tcc = stMeteo$tcc,
    ws = stMeteo$ws,
    option = "custom"
)

Vector field plot

Description

Simple function to plot a vector field given two components.

Usage

vectorField(
  data,
  scale = 1,
  everyx = 1,
  everyy = 1,
  size = 0.25,
  preview = TRUE
)

Arguments

Details

This function plots a vector field given a data.frame with coordinates (x, y) and corresponding velocity components (u, v). Vectors are coloured by magnitude (speed). The coordinates are assumed to be on a regular rectangular grid in the UTM reference system.

This function is heavily inspired by snippets of code in R Graphics Cookbook by Winston Chang (https://r-graphics.org/index.html).

Value

A ggplot2 object if preview = TRUE. A ggplot2 a plot, as a contourPlot2() and the vector field will be overlapped.

Examples

## Not run: 
metU <- importADSOBIN(
    "/path/to/meteofile",
    variable = 'U',
    slice = 2,
    k = 1000,
    verbose = TRUE
)
metU <- as.data.frame(metU)
metU <- metU %>%
    mutate(u = z, z = NULL)

metV <- importADSOBIN(
    "/path/to/meteofile",
    variable = 'V',
    slice = 2,
    k = 1000,
    verbose = TRUE
 )
metV <- as.data.frame(metV)
metV <- metV |>
    mutate(v = z, z = NULL)

met <- merge(metU, metV, by = c("x", "y"))

vectorField(
    met,
    everyx = 2,
    everyy = 2,
    scale = 10
) +
    coord_fixed(ratio = 1, xlim = c(0, 1000), ylim = c(0, 1000)) +
    scale_color_viridis_c()

# Overlap the vector field to a contour plot and set vector colours to black
met$ws <- sqrt(met$u^2 + met$v^2)
contourPlot2(met, z = "ws") +
    vectorField(
        met,
        everyx = 2,
        everyy = 2,
        scale = 10,
        preview = FALSE
    ) +
    scale_colour_gradient(low = "black", high = "black", guide = NULL)

## End(Not run)