if (FALSE) {
library(RNCEP)
######################################################
## In the first example, we use the internal R function
## 'mean' to calculate average temperatures per
## month and year
## First gather temperature data from a spatial extent
## for January and February from 2000-2001.
wx.extent <- NCEP.gather(variable='air', level=850,
months.minmax=c(1,2), years.minmax=c(2000,2001),
lat.southnorth=c(50, 55), lon.westeast=c(0, 5),
reanalysis2=FALSE, return.units=TRUE)
## Now calculate the average temperature per month and year ##
wx.ag <- NCEP.aggregate(wx.data=wx.extent, YEARS=TRUE, MONTHS=TRUE,
DAYS=FALSE, HOURS=FALSE, fxn='mean')
## Notice that aggregated time components have been replaced
## with "XX" or "XXXX" ##
dimnames(wx.ag)[[3]][1]
#######################################################
## In the second example, we create our own function to
## calculate the percentage of observations at each grid
## point with a temperature less than -5 degrees Celsius ##
## First create the function ##
## Note: temperature is in degrees Kelvin in NCEP database ##
COLD <- function(x){
return(length(which(x < 273.15-5))/(length(x) - sum(is.na(x))))
}
## Now calculate the percentage of occuence of temperatures
## less than -5 degrees Celsius per month and year ##
wx.cold <- NCEP.aggregate(wx.data=wx.extent, YEARS=TRUE, MONTHS=TRUE,
DAYS=FALSE, HOURS=FALSE, fxn='COLD')
##########################################################
## As explained in the Details section above,
## calculating some variables requires sequential aggregations ##
## Here we calculate the monthly mean of daily average
## relative humidity.
## First gather relative humidity data from near the surface
## for a spatial extent for October through November
## from 2001-2002.
wx.extent <- NCEP.gather(variable='rhum.sig995', level='surface',
months.minmax=c(10,11), years.minmax=c(2001,2002),
lat.southnorth=c(50, 55), lon.westeast=c(0, 5),
reanalysis2=FALSE, return.units=FALSE)
## First calculate maximum daily relative humidity ##
wx.ag <- NCEP.aggregate(wx.data=wx.extent,
HOURS=FALSE, fxn='max')
## Then calculate the monthly average of those daily maximums ##
wx.ag2 <- NCEP.aggregate(wx.data=wx.ag,
DAYS=FALSE, fxn='mean')
##########################################################
## Data that have been aggregated may then be visualized
## or exported in various formats ##
###########################################
## Visualize the aggregated temperatures ##
NCEP.vis.area(wx.ag2, title.args=
list(main='Monthly average of daily maximum relative humidity
\n October 2000'), image.plot.args=
list(legend.args=list(text='%',
cex=1.15, padj=-1, adj=-.25)))
###########################################################
## Export a layer of the data to a format that can then be
## imported into ArcGIS ##
## Convert the first layer of the aggregated array
## to a data.frame ##
wx.df <- NCEP.array2df(wx.ag[,,1])
## Specify that the data.frame is a spatial object
library(sp)
coordinates(wx.df) <- ~longitude+latitude
gridded(wx.df) <- TRUE
proj4string(wx.df) <- CRS('+proj=longlat + datum=WGS84')
## Save the data in .asc format
write.asciigrid(wx.df, fname='wx.asc')
## Note: Data will be written to your working directory ##
## The resulting .asc file can be imported into ArcMap
## using ArcMap's "ASCII to Raster" tool in the "Conversion Tools"
## section of the ArcToolbox. ##
}
Run the code above in your browser using DataLab