How to get daily data

[QuocPhamBao]

Hello
How can I modify this line to get the daily data from NCEP?
t2m <- t2m.NCEP(lon=c(0,360),lat=c(35,90))
If there is no way to get the daily data from NCEP, can we get daily data from ERA5 or some other datasets?

Best,

[brasmus]

Hi. The data that come with esd are only sample data and typically monthly or annual (to keep its size down). But you can download daily data with ERA5.CDS() after you have registered and installed the right python tools, and then read the the daily data using retrieve. I think you also can download daily data from e.g. http://apdrc.soest.hawaii.edu/datadoc/ncep2_daily.php and then read it using retrieve.

[QuocPhamBao]

Hello,
Thank you for your answer, I would like to ask is there any way to get daily data in some specific dates.
For example, I only want to get date in these days: "1966-01-08", "1966-01-10","1966-01-13",
"1966-01-19", for further analysis.
I used below lines:

rr <- retrieve('hgt.1966.nc', lev = 500,
it = as.Date(c('1966-01-08', '1966-01-19'))) # this worked well and it extract daily data from 1966-01-08 to 1966-01-19.

Then, I want to extract some specific dates in that period for further analysis:

rr = subset(rr, it = as.Date(c("1966-01-08", "1966-01-10","1966-01-13",
"1966-01-19"))) # but this is not working.
The "subset" is not work as I expected.
Please give me suggestion.

Best

[brasmus]

Here is the way to extract only one day:

rr <- retrieve('hgt.1966.nc', lev = 500)
rr1 <- subset(rr, it = c("1966-01-08", "1966-01-08"))
rr2 <- subset(rr, it = c("1966-01-10", "1966-01-10"))
rr3 <- subset(rr, it = c("1966-01-13", "1966-01-13"))
rr123 <- c(rr1,rr2,rr3)

You don't need to include as.Date in the subset-call. It takes the string as a date...

[kajsamp]

I'm guessing that the index of your object 'rr' is of a different class (POSIXt?) than your input 'it'. Is that right? Could you check the following, please?

class(index(rr))
index(rr)[1:10]

I'm working on updating the function subset to make it more flexible with regards to time formats, but here's a workaround in the meantime:

x <- subset(rr, it=as.Date(index(rr)) %in% as.Date(c("1966-01-08", "1966-01-10","1966-01-13",
"1966-01-19")))

[QuocPhamBao]

Hello
Thank you for your response.
Please see the outcome as below:

class(index(rr))
[1] "POSIXct" "POSIXt"
index(rr)[1:10]
[1] "1966-01-08 UTC" "1966-01-09 UTC" "1966-01-10 UTC" "1966-01-11 UTC"
[5] "1966-01-12 UTC" "1966-01-13 UTC" "1966-01-14 UTC" "1966-01-15 UTC"
[9] "1966-01-16 UTC" "1966-01-17 UTC"

[QuocPhamBao]

Hello
There is an error when I applied EOF

r1 <- retrieve('hgt.1966.nc', lev = 500, it = as.Date(c('1966-01-08',
'1966-01-17')),
lon=c(-180,180),lat=c(35,90))

r2 <- retrieve('hgt.1967.nc', lev = 500, it = as.Date(c('1967-12-19','1967-12-19')),
lon=c(-180,180),lat=c(35,90))

x1 <- subset(r1, it=as.Date(index(r1)) %in% as.Date(c("1966-01-08", "1966-01-09","1966-01-10", "1966-01-12","1966-01-13",
"1966-01-17", "1966-01-19")))
x2 <- subset(r2, it=as.Date(index(r2)) %in% as.Date(c("1967-12-19")))

x = c(x1, x2)

e_o_f = EOF(x)

Error in EOF.default(x) : object 'eof' not found

I guess the problem maybe from: x = c(x1, x2) because if I only use x1: e_o_f = EOF(x1)
It works well.

[kajsamp]

The problem comes from the way you combine x1 and x2. Try using cbind(x1, x2) instead.

[QuocPhamBao]

Yes, I also tried that way, but there is another error:

Error in svd(y, nu = min(c(ny, npca)), nv = min(c(ny, npca))) :
a dimension is zero
Error in svd(t(y), nu = min(c(ny, npca)), nv = min(c(ny, npca))) :
a dimension is zero
Error: $ operator is invalid for atomic vectors

[kajsamp]

Sorry, it should be rbind(x1, x2), not cbind.

You may still have trouble with the EOF analysis with such a small input object, containing only 8 time steps. What are you trying to figure out with the EOF analysis?

[QuocPhamBao]

Thank you very much.
It worked !
I am intending to collect geopotential height, sea level pressure and wind speed at 500hPa level.
I need to collect for 96 events (96 specific dates) from 1966-2013.

There is a limitation of daily NCEP data that I need to download for each year (1966, 1967,...2013) and extract data in specific date of each year. Then combine them (using rbind), finally I use EOF to see the pattern of the first principle component of these variables (geopotential height, sea level pressure and wind speed) in 96 events.
This is really time consuming way !
I hope I can download the data of the period in single *.nc file.
But I cannot find any website that I can download NCEP data for the whole period (1966-2013).
So, I need to download year by year and read one by one into R.

Is it possible for me to use PCA for that matrix? or PCA is only work for data in specific station?
Maybe EOF and PCA is quite similar and it both can apply for my purpose ?

Best

[kajsamp]

The PCA and EOF methods are implemented in esd in similar ways, using single value decomposition (svd), but EOF is for field objects and PCA for station objects. For your purposes, you should use the EOF function.

Regarding the many NCEP files, you could look into using cdo to extract data for the relevant dates and combine them into one netcdf file. Otherwise if you do this in R, I suggest saving the data in an rda-file so that you only have to do it once.

Sounds like you have an interesting project going on. Good luck and let us know if you have any further issues or questions about the esd-package.

[QuocPhamBao]

Thank you !
May I ask one more question.
In the example 2.6, I think it will plot the figure of PC1 (first principle component), could we plot the figures of other PC, let's say PC2, PC3,PC4 and PC5?
Example 2.6.
#Get NCEP 2m air temperature for the selected spatial window defined by lon and lat
t2m <- t2m.NCEP(lon=c(-30,30),lat=c(40,70))

Computes the EOFs

X <- EOF(t2m)

Plot the result

plot(X)

[brasmus]

Sure! plot(X,ip=2) to show the 2nd component (ip=3 shows the third, etc). We use a standard notation for arguments, where is refers to 'index space'. it is 'index time', and ip is index pattern.

[QuocPhamBao]

Thank you very much for your answer.
It helped me a lot.
Currently, can we create maps of EOFs in azimuthal projection ?

[QuocPhamBao]

Hello !
I would like to ask if possible to extract the values of 20 leading EOFs to see its values?
For the figure of ght (please see the attachment), is it possible to extract the data (longitude, latitude and its values) in that figure in order to draw in other application?
Best regards

[brasmus]

Yes, you see the values of the EOFs by addressing them attr(x,'pattern') if x <- EOF(...). Alternatively, if you want just one pattern, you can use

data("slp.ERA5")
eof <- EOF(slp.ERA5)
m1 <- subset(eof,ip=1)
map(m1) -> y
y

Many of the results are stored as attibutes - try names(attributes(eof)). You can also use str(eof) to show the structure.

[QuocPhamBao]

Thank you for great support. I understand how to extract the value of eof.

In the below figure, is there any way for us to know how many percentage that component 1 accounts for, component 2 accounts for, .....component 20 accounts for.
I mean if I want to take the values and plot in another software, how I can see the values.

Best