Add a layer to TrajLevel function of Openair

[jordinagili]

I am using the TrajLevel function of Openair and I would like to add a layer with fire information (which I have in shp format) to the resulting map. Would you know if is that possible?

[jack-davison]

Hi Jordina,

{openairmaps} now has an (experimental) {ggplot2} implementation of trajLevel(). How are you with {ggplot2}? It will let you add on extra layers much more straightforwardly. The {sf} package will let you read in your shapefile, and then the geom_sf() function will let you add on your extra layer.

library(openairmaps) # openair maps
library(sf) # work with shape files
#> Warning: package 'sf' was built under R version 4.2.2
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(ggplot2) # ggplot2 extensions

# get a shapefile
# (here I mock one up)
st <- data.frame(lat = c(55, 56, 57), lng = c(0, 1, 2))
sf <- st_as_sf(st, coords = c("lng", "lat"), crs = st_crs(4326))

# create plot
plt <-
  trajLevelMapStatic(traj_data, pollutant = "pm2.5", statistic = "cwt") +
  scale_fill_gradientn(colours = openair::openColours())

# view
plt

# add our extra shapefile
plt +
  geom_sf(data = sf, inherit.aes = FALSE) +
  coord_sf(xlim = c(-55, 18), ylim = c(48, 82),
           crs = st_crs(4326), default_crs = st_crs(4326))
#> Coordinate system already present. Adding new coordinate system, which will
#> replace the existing one.

Alternatively, you could just run your trajLevel() function and pull out the data, and re-plot it however you'd like!

# or just pull out the data and plot however you'd like:
plt$data
#> # A tibble: 712 × 8
#> # Groups:   xgrid, ygrid [712]
#>    xgrid ygrid default                  N date                count pm2.5 .group
#>    <dbl> <dbl> <fct>                <int> <dttm>              <dbl> <dbl>  <int>
#>  1   -49    57 15 April 2010 to 21…     2 2010-04-17 03:00:00  15    0.75      1
#>  2   -49    58 15 April 2010 to 21…     1 2010-04-17 06:00:00  22    1.1       2
#>  3   -48    55 15 April 2010 to 21…     1 2010-04-17 00:00:00  14    0.7       3
#>  4   -48    57 15 April 2010 to 21…     2 2010-04-17 03:00:00  15    0.75      4
#>  5   -48    58 15 April 2010 to 21…     3 2010-04-17 06:00:00  20.7  1.03      5
#>  6   -48    61 15 April 2010 to 21…     1 2010-04-18 06:00:00  32    1.6       6
#>  7   -47    56 15 April 2010 to 21…     2 2010-04-17 00:00:00  14    0.7       7
#>  8   -47    57 15 April 2010 to 21…     1 2010-04-17 03:00:00  15    0.75      8
#>  9   -47    58 15 April 2010 to 21…     2 2010-04-17 06:00:00  20    1         9
#> 10   -47    59 15 April 2010 to 21…     1 2010-04-18 09:00:00  34    1.7      10
#> # … with 702 more rows

^{Created on 2023-03-12 with reprex v2.0.2}

Further, if you're happy with an interactive figure, you can use openairmaps::trajLevelMap() which is built in {leaflet}. That'll make it even easier to add extra layers on using shapefiles, with the extra utility of e.g. turning layers on and off.

If you have any more questions do shout, or feel free to shoot over your traj data and shapefile and we can take a look.

Cheers, Jack

[jordinagili]

Hi Jack, thank you very much for your detailed answer. I appreciate it.

The code I am using is the following:

trajLevel(subset(traj,lon > -10 & lon < -6 & lat > 40.5 & lat < 45), pollutant = "avg", statistic="cwt",col = c("skyblue", "yellow", "red"),border = "white", smooth = TRUE, lon.inc = 0.1, lat.inc = 0.1, map.res = "hires", projection = "lambert", grid.col = "darkgray", key.header = "CWT of PM2.5 (ug/m3)", key.position = "bottom",percentile="80", map.fill = FALSE, parameters = c(-10, 100))

And I would like to add to it the extra shp layer attached here effis.zip

Thanks,

Jordina

[jack-davison]

Hi again Jordina,

I don't believe you've provided traj so I've had to make some data up.

Is this the sort of thing you're after?

library(openairmaps)
library(openair)
library(ggplot2)
library(sf)
#> Warning: package 'sf' was built under R version 4.2.2
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

# read shapefile
sf <- read_sf("~/effis/effis4.shp")

# import traj
traj <- importTraj("elche")

# need a pollutant
traj$avg <- rnorm(nrow(traj), mean = 50, sd = 5)

trajLevelMapStatic(
  subset(traj, lon > -10 & lon < -6 & lat > 40.5 & lat < 45),
  pollutant = "avg",
  statistic = "cwt",
  lon.inc = 0.1,
  lat.inc = 0.1,
  percentile = "80",
  map.fill = NA
) +
  geom_sf(
    data = sf,
    inherit.aes = FALSE,
    fill = NA,
    color = "black"
  ) +
  scale_fill_gradientn(colors = c("skyblue", "yellow", "red")) +
  labs(fill = quickText("CWT of \nPM2.5 (ug/m3)")) +
  coord_sf(
    xlim = c(-10,-6),
    ylim = c(40.5, 45),
    crs = st_crs(4326),
    expand = FALSE
  )
#> Warning in min(x): no non-missing arguments to min; returning Inf
#> Warning in max(x): no non-missing arguments to max; returning -Inf
#> Coordinate system already present. Adding new coordinate system, which will
#> replace the existing one.

^{Created on 2023-03-12 with reprex v2.0.2}

[jordinagili]

Yes, it is! Thank you very much. I see that the smooth argument is not valid for ggplot2 as I am getting this error:

Error in ggplot2::coord_sf(xlim = xlim, ylim = ylim, default_crs = sf::st_crs(4326), : unused argument (smooth = TRUE)

Is not possible to smooth then?

Thanks,

Jordina

[jack-davison]

I'm afraid right now its not available - it acts as a wrapper around openair::trajLevel() at the moment and openair doesn't return the smoothed values it uses to plot. I'll have a word with David to see if the smoothed values can be returned.

I'll also confirm whether a shapefile can be added to the original trajLevel() function.

Thanks, Jack

[jordinagili]

Thanks, I'll appreciate it.

Thank you very much,

Jordina

[jack-davison]

Hi again Jordina,

smooth should now be an option for the static maps - do download the dev version of {openairmaps}. Note that it may be slower to draw as there are far more rectangles to map to the chosen coordinate system.

openairmaps::trajLevelMapStatic(
  openairmaps::traj_data,
  pollutant = "pm10",
  statistic = "sqtba",
  smooth = TRUE
) +
  ggplot2::scale_fill_gradientn(colors = c("skyblue", "yellow", "red"))

[jordinagili]

Great! Thank you very much,

Jordina

[jordinagili]

Hi again Jack,

Now I am trying to overlay the multisite SQTBA plot with my fire shapefile but I am getting this error:

Error in seq.default(h[1], h[2], length = n) : 'length.out' must be a non-negative number In addition: Warning messages: 1: In min(breaks) : no non-missing arguments to min; returning Inf 2: In max(breaks) : no non-missing arguments to max; returning -Inf

The trajLevel function works correctly with this multisite_test.zip, but the trajLevelMapStatic function does not. Here the code:

trajLevelMapStatic( subset(multisite_test, lon > -10 & lon < -6 & lat > 40.5 & lat < 45), pollutant = "avg", statistic = "sqtba", lon.inc = 0.1, lat.inc = 0.1, percentile = "80", map.fill = NA, smooth = TRUE, .combine = "site" )

Does the trajLevelMapStatic function works with multisite trajectory data?

Thank you very much for all your help,

Jordina

[jack-davison]

Ah, I think that is an oversight in the way I wrote the code! I'll have some time to give this a look tomorrow, I hope, and will keep you updated.

[jack-davison]

Okay, hopefully now we're away! Please download the dev version of openairmaps again and give this a go.

On my PC it takes about 40 seconds to run, so please be patient with it as there's a lot of computation going on!

library(openairmaps)
library(openair)
library(sf)
library(ggplot2)

multisite_test <- readr::read_rds("../../multisite_test.rds")

sf <- sf::read_sf("../../effis/effis4.shp")

trajLevelMapStatic(
  subset(multisite_test, lon > -10 & lon < -6 & lat > 40.5 & lat < 45),
  pollutant = "avg",
  statistic = "sqtba",
  lon.inc = 0.1,
  lat.inc = 0.1,
  percentile = "80",
  map.fill = NA,
  smooth = TRUE,
  .combine = "site",
  crs = NULL
) +
  geom_sf(
    data = sf,
    inherit.aes = FALSE,
    fill = NA,
    color = "black"
  ) +
  scale_fill_gradientn(colors = c("skyblue", "yellow", "red")) +
  labs(fill = quickText("CWT of \nPM2.5 (ug/m3)")) +
  coord_sf(
    xlim = c(-10,-6),
    ylim = c(41, 45),
    crs = st_crs(4326),
    expand = FALSE
  )

[jordinagili]

Thank you very much, Jack! I really appreciate it.

Jordina

[jordinagili]

Hi again Jack,

Would it be possible to add the shp layer to the trajPlot function?

selectByDate(traj_PA7, start = "15/07/2022", end = "15/07/2022", hour = 11:12 ) %>% trajPlot(pollutant = "PA7", col = "turbo", lwd = 2,
xlim = c(-10, -3), ylim = c(39, 45), map.res = "hires", projection = "lambert", grid.col = "darkgray", parameters = c(-10, 100))

Thank you!

Jordina

[jack-davison]

Hi again Jordina,

{openairmaps} also has a trajMapStatic() function which you can add to in the same way to the trajLevelMapStatic() function.

openairmaps::trajMapStatic(openairmaps::traj_data,
                           colour = "pm2.5") +
  ggplot2::scale_color_viridis_c(name = openair::quickText("PM2.5"),
                                 option = "turbo") +
  ggplot2::theme(panel.grid = ggplot2::element_line(color = "darkgrey"))

^{Created on 2023-05-26 with reprex v2.0.2}

[jordinagili]

Many thanks Jack!

Jordina