0

I currently have a series of air mass back trajectory plots where are grid by the frequency of back trajectory crossings, produced using the OpenAir package for R. I would like to add to these plots points which indicate the location of fires. The points read from a .csv file where each unique fire has a latitude and longitude value.

I have read the OpenAir manual, but have not seen any solution presented. I have tried to add the points with the following code:

Week_One_Freq + points (FireData$longitude, FireData$latitude, col = "red", cex = .6)

But encounter the error: 

Error in Week_One_Freq + points(FireData$longitude, FireData$latitude,:  non-numeric argument to binary operator

I have tested to see if the issue is with the fire data. I have plotted the points as follows: 

library(rworldmap)
newmap <- getMap(resolution = "low")
plot(newmap, xlim = c(-80, 80), ylim = c(-40, 40), asp = 1)
points (FireData$longitude, FireData$latitude, col = "red", cex = .6)

And it works fine, correctly plotting fire location.

Here is my code for the frequency trajectory plots:

Week_One_Freq <- trajLevel(Week_1_Traj, statistic = "frequency",
                          smooth = TRUE,
                          map.cols = "grey40", map.alpha = 0.2,
                          projection = "mercator", parameters = NULL,
                          grid.col = "transparent",
                          key.pos = "right", 
                          orientation = c(90,0,0),
                          npoints = 20, origin = TRUE,
                          xlim= c(-80, 80), ylim = c(-40, 40))

So far, I can plot fire data points on a map, and I can produce a series of frequency plots for air mass back trajectories. I would like to combine the two, so that on one map there is both back trajectories and points indicating the location of fires.

EDIT: Example of trajectory data used to create frequency plots:

structure(list(X = 1:5, Receptor = c(1L, 1L, 1L, 1L, 
1L, 1L), Year = c(2017L, 2017L, 2017L, 2017L, 2017L, 2017L), 
Month = c(10L, 10L, 10L, 10L, 10L, 10L), Day = c(6L, 5L, 
5L, 5L, 5L, 5L), Hour = c(0L, 23L, 22L, 21L, 20L, 19L), hour.inc = 0:-5, 
lat = c(-3.162, -3.189, -2.937, -2.87, -2.809, -2.753), lon = c(-60.6, 
-60.439, -60.282, -60.13, -60.986, -60.851), height = c(500, 
502.2, 500.6, 495.2, 486.4, 474.3), pressure = c(944.2, 942.9, 
942.4, 942.3, 942.7, 943.5), date = structure(c(1L, 1L, 1L, 
1L, 1L, 1L), .Label = c("06/10/2017 00:00", "06/10/2017 03:00", 
"06/10/2017 06:00", "06/10/2017 09:00", "06/10/2017 12:00", 
"06/10/2017 15:00", "06/10/2017 18:00", "06/10/2017 21:00", 
"07/10/2017 00:00", "07/10/2017 03:00", "07/10/2017 06:00", 
"07/10/2017 09:00", "07/10/2017 12:00", "07/10/2017 15:00", 
"07/10/2017 18:00", "07/10/2017 21:00", "08/10/2017 00:00", 
"08/10/2017 03:00", "08/10/2017 06:00", "08/10/2017 09:00", 
"08/10/2017 12:00", "08/10/2017 15:00", "08/10/2017 18:00", 
"08/10/2017 21:00", "09/10/2017 00:00", "09/10/2017 03:00", 
"09/10/2017 06:00", "09/10/2017 09:00", "09/10/2017 12:00", 
"09/10/2017 15:00", "09/10/2017 18:00", "09/10/2017 21:00", 
"10/10/2017 00:00", "10/10/2017 03:00", "10/10/2017 06:00", 
"10/10/2017 09:00", "10/10/2017 12:00", "10/10/2017 15:00", 
"10/10/2017 18:00", "10/10/2017 21:00", "11/10/2017 00:00", 
"11/10/2017 03:00", "11/10/2017 06:00", "11/10/2017 09:00", 
"11/10/2017 12:00", "11/10/2017 15:00", "11/10/2017 18:00", 
"11/10/2017 21:00", "12/10/2017 00:00", "12/10/2017 03:00", 
"12/10/2017 06:00", "12/10/2017 09:00", "12/10/2017 12:00", 
"12/10/2017 15:00", "12/10/2017 18:00", "12/10/2017 21:00"
), class = "factor"), Date = structure(c(1L, 1L, 1L, 1L, 
1L, 1L), .Label = c("06/10/2017", "07/10/2017", "08/10/2017", 
"09/10/2017", "10/10/2017", "11/10/2017", "12/10/2017"), class = 
"factor")), row.names = c(NA, 
6L), class = "data.frame")

And example of fire data:

structure(list(acq_date = structure(c(6L, 6L, 6L, 6L, 6L, 6L), .Label = 
c("01/11/2017", "02/11/2017", "03/11/2017", "04/11/2017", "05/11/2017", 
"06/10/2017", "06/11/2017", "07/10/2017", "08/10/2017", "09/10/2017", 
"10/10/2017", "11/10/2017", "12/10/2017", "13/10/2017", "14/10/2017", 
"15/10/2017", "16/10/2017", "17/10/2017", "18/10/2017", "19/10/2017", 
"20/10/2017", "21/10/2017", "22/10/2017", "23/10/2017", "24/10/2017", 
"25/10/2017", "26/10/2017", "27/10/2017", "28/10/2017", "29/10/2017", 
"30/10/2017", "31/10/2017"), class = "factor"), latitude = c(-0.18499, 
-0.18777, -0.18832, -0.1932, -0.19654, -0.21077), longitude = c(-48.89484, 
-48.89885, -48.89536, -49.06121, -49.06173, -48.97051)), row.names = c(NA, 
6L), class = "data.frame")
R Ramsay
  • 130
  • 7
  • 1
    Please edit your question and add some data (using `dput()`) as to make it reproducible. – erc Jan 29 '19 at 13:32

1 Answers1

1

Someone has solved the same issue at https://github.com/davidcarslaw/openair/issues/68. You can add points to the last trellis object using the layer function in latticeExtra, but first you need to get your points onto the correct projection using mapproject. I've expanded on their example to include some additional information like the required libraries.

library(openair)
library(maps)
library(mapproj)
library(latticeExtra)

# import your data
Trajectory_Data <- read.table("Analysis.txt")
Fire_Data <- read.csv("Fire_Data.csv")

lon <- Fire_Data$longitude  
lat <- Fire_Data$latitude  
#convert points to projected lat/long
tmp <- mapproject(x=lon, y=lat, projection='lambert', parameters = c(51, 51), orientation = c(90, 0, 0))  
# I found the projection, parameters, and orientation by looking at the default values for arguments in the trajLevel function 

# Plot trajectories
trajLevel(subset(Trajectory_Data, lon > -70 & lon < -20 & lat > -15 & lat < 15),
          pollutant = "BC", statistic = "cwt", smooth = TRUE, col = "increment", border = NA)
# add the points. you can adjust the size of the points with cex
c <- trellis.last.object() + layer(lpoints(tmp$x, tmp$y, pch=16, cex=0.05, col="black"))
print(c)
cumulonimbus
  • 26
  • 2