davidcarslaw / openair

Tools for air quality data analysis
https://davidcarslaw.github.io/openair/
GNU General Public License v2.0
299 stars 113 forks source link

Black trajectory functions for other sites #35

Closed lyggd closed 8 years ago

lyggd commented 8 years ago

Hi David, I find black trajectory functions is very helpful but it is only works for limited sites. I am working in the US and will have the Hysplist results and pollutant data for some sites. I would like to check with you if there are functions in open air package can help me to do the PSCF analysis ? Happy new year ! Thank you very much. Yi

JaysonAP commented 8 years ago

Do you have the trajectories loaded into R, or you just have the trajectories direct from Hysplit?

I believe you will need to process trajectories within OpenAir so they are read properly in R. Look into the getMet and procTraj functions, and read Appendix D in the OpenAir manual. This should provide you with all the steps necessary. This process will allow you to create back trajectories for any lat/lon you would like.

Once you have your trajectories in OpenAir, or if you already do, then Chapter 26 of the OpenAir manual will be of benefit. PSCF analysis is rain with the format:

trajLevel(traj, pollutant = "pm2.5", statistic = "pscf", col = "increment")

lyggd commented 8 years ago

Thank you JaysonAP. I make some progress based on your suggestions which is very exciting. Best Regards, Yi

JaysonAP commented 8 years ago

Not a problem. I have been through the process for a number of different locations. Be sure to ask if you have any additional questions.

lyggd commented 8 years ago

Hi JaysonAP, In Appendix D of Open air mannal, ProcTraj. function has one line code :
x <- "c:\\users\\david\\hysplit4\\exec\\hyts_std"

Should I change it based on my setting ? For example: x <- "c:\\hysplit4\\exec\\hyts_std"

I saved read.files, add.met, procTraj function in myFunction.r . Met data had been successfully downloaded through getMet For example: I would like to run 24hr analysis for one year (2014) at JLG Supersite site in AZ.

Then my code is

<library("openair") source('M:/R/myFunction.r') for (i in 2014) { procTraj(lat = 33.5038, lon = -112.0957, year = i, name = "JLG Supersite", hours = 24, met = "M:/HYSPLIT/Mete/TrajData/", out = "M:/HYSPLIT/output/TrajProc/", hy.path = "c:/hysplit4/") }>`

The code ran 20 mins and stop with Error I upload the results in the attachment. Do you know how to solve this problem ? Should I do something for 'The system cannot find the path specified' ? capture

Thank you very much for your help. I appreciated very much. Best Regards, Yi

JaysonAP commented 8 years ago

Yes, you will need to change any file paths to match what is on your local machine. Double check your HYSPLIT path as well.

I always find it best to review any files before I run them to ensure that paths and information are correct.

Also, if you ever believe you will need trajectories longer than 24 hours you can always run your myFunction.r with hours=120 and then when you begin to do your PSCF analysis use a subset to view only 24 hours.

Example:

JLGTraj_24 <- subset(JLGTraj, hour.inc > -25)

lyggd commented 8 years ago

Hi JaysonAp, Thank you for your help. With your help, my Hysplit start running. My code is

library("openair")
library("plyr")
library("lattice")
library("latticeExtra")
library('maps')
library('mapproj')
library('mapdata')
source('M:/R/myFunction.r')
rm(list=ls())
for (i in 2014) {
  procTraj(lat = 33.5038, lon = -112.0957, year = i,
           name = "JLG Supersite", hours = 120,
           met = "M:/HYSPLIT/Mete/TrajData/", out = "M:/HYSPLIT/output/TrajProc/",
           hy.path = "c:/hysplit4/")
}

traj <- importTraj(site = "JLG Supersite", year = 2014, 
                   local ='M:/HYSPLIT/output/TrajProc/' )
JLGTraj_24 <- subset(traj, hour.inc > -25)
JLGTraj_24$day <- as.Date(JLGTraj_24$date)
trajPlot(selectByDate(JLGTraj_24, start = "15/8/2014", end ="21/8/2014"), lon = "lon", lat = "lat",
         map = TRUE ,map.fill = F, orientation = c(45,270,0), map.res = "state",
         grid.col = F, type = "day", col = "jet", lwd = 2, key.pos = "right",
         key.col = 1)

And I get the figure below rplot

But I still have some questions want to check with you: 1: What is orientation = c(45,270,0) ? I check the help file but still confused.

2: From the figure above, we find there are 8 tracks for each day. how to explain those 8 tracks ? Could you tell me where I can set the number of tack ?

  1. I want this map below showing the boundary of California, Nevada and Mexico. Could you tell me what should I do ?

I find some code (not include Mexico) but do not know how combine them together. map('state', region = c('Arizona', 'California', 'Nevada''),boundary = FALSE, lty = 2, add = TRUE)

rplot02

Again, I appreciate your patient, fast response and help very much. Best Regards, Yi

JaysonAP commented 8 years ago

1) OpenAir utilizes the maps package for trajectory plots, orientation is a function of that package where orientation = c(latitude, longitude, rotation) to set the center point and rotation of the map. For more on the maps package, see link: R Maps Package

2) The 8 tracks per day is the function of the HYSPLIT model that you used (I didn't see it in any previous posts) but I presume you used a NAM/WRF or NARR for your HYSPLIT which has a 3 hour temporal resolution. A 3 hour temporal resolution means there would be 8 trajectories per day. If you wanted to only use a specific time you could subset or cut data to a specific time (i.e. at 0600).

3) You may be plotting all of those boundaries, but your map is zoomed in to the extent to show only your trajectories. Also, map.res="state" or map.res="hires" will show all state boundaries by default. Try setting your map to a specific lat/lon and projection using xlim, ylim and projection, example:

trajPlot(selectByDate(JLGTraj_24, start = "15/8/2014", end ="21/8/2014"), lon = "lon", lat = "lat", map = TRUE ,map.fill = F, map.res = "state", projection="albers", xlim=c(-125,-110), ylim=c(30,45), grid.col = F, type = "day", col = "jet", lwd = 2, key.pos = "right", key.col = 1)

lyggd commented 8 years ago

Hi JaysonAP, Thank you very much for your information.

But I am still confuse on the 2nd question. I download monthly meteorological (.gbl) files from the NOAA website by using getMet function, which is NCEP/NCAR Reanalysis data.

1: What do you mean “has a 3 hour temporal resolution? A 3 hour temporal resolution means there would be 8 trajectories per day.”?

That is what I understand: take 2014-08-17 as example, there are 8 tracks on the map starting from JLG site. The first track start from JLG site at 00:00, 2014-08-17 and end at 00:00 on 2014-08-16. The second track start from JLG site at 3:00, 2014-08-17 and end at 3:00 on 2014-08-16, and so on. Is it correct?

2: I find if I choose different Met data, the Hysplit results will be different. See the figure below ( Right one is using GDAS and left one is using NCEPNCAR Reanalysis).

Do you know which Met data I should choose if I want to run Hysplit for O3 in AZ ?

capture

3: My final goal is doing some statistical analysis using the openair functions such as PSCF, CWT.

The next problem I have is how to pair Hysplit results with Hourly air quality data. For this purpose, what I should set the Hysplit before running?

It looks like we can import air quality data in British using importAURN function. Do you know if there is any way to directly download air quality data for US site? Otherwise, I can manual download data from AQS website https://aqs.epa.gov/api

Thank you very much for your help. i appreciate very much. Best Regards, Yi

lyggd commented 8 years ago

Hi Jayson, Hope you still following our conversion. I still have some questions need your help. Thank you very much. Best Regards, Yi

JaysonAP commented 8 years ago

Yes, NCAR Reanalysis data will have 8 'start times' each day because of the 3 hour temporal resolution. Different models can have different temporal resolution, so you can simply choose which is best for your research/study. Reanalysis data, specifically the NARR, is likely the best option.

To merge air quality data you will likely run a join command (see dplyr and tidyr packages) to combine the HYSPLIT data with air quality data. Be sure when you are merging that you merge using the by = "date" column to match up the air quality data with proper HYSPLIT time.

For importing US AQ data directly into R try using the raqdm package - please take time to read through the users guide, it will provide you with much of the information you need to operate the package successfully.

umeshdumka commented 8 years ago

Hi David and JaysonAP, I find black trajectory functions is very helpful. I would like to generate the trajectory for my location but after running the code it is giving the following error

Error in $<-.data.frame(*tmp*, "date", value = c(1420070400, 1420070400, : replacement has 493311 rows, data has 493268

Please let me know how to resolve this error.

In addition to this I want to generate the trajectory in hourly basis then please let me know where I have to change in the code "procTraj"

Thank you very much. @davidcarslaw @JaysonAP