if (FALSE) {
library(RNCEP)
## Query the temperature from a particular pressure level ##
wx.extent1 <- NCEP.gather(variable='air', level=850,
months.minmax=c(9,10), years.minmax=c(1996,1997),
lat.southnorth=c(50,55), lon.westeast=c(5,10),
reanalysis2 = FALSE, return.units = TRUE)
## Query the temperature at 2 meters altitude with reference to
## the surface
wx.extent2 <- NCEP.gather(variable='air.sig995', level='surface',
months.minmax=c(2,3), years.minmax=c(2000,2001),
lat.southnorth=c(50,55), lon.westeast=c(0,5),
reanalysis2 = FALSE, return.units = TRUE)
## Query the temperature at 2 meters altitude with reference to
## a T62 Gaussian grid
wx.extent3 <- NCEP.gather(variable='air.2m', level='gaussian',
months.minmax=c(4,5), years.minmax=c(2006,2007),
lat.southnorth=c(32,35), lon.westeast=c(-35,-32),
reanalysis2 = FALSE, return.units = TRUE)
## Note that the dimnames of the data array indicate the
## latitudes, longitudes, and datetimes of the data. ##
dimnames(wx.extent1)
## Therefore, the latitudes, longitudes, and datetimes
## can be called. ##
dimnames(wx.extent1)[[1]] ## latitudes
dimnames(wx.extent1)[[2]] ## longitudes
dimnames(wx.extent1)[[3]] ## datetimes
#################################################################
#################################################################
#################################################################
## THERE ARE MANY OPTIONS FOR CREATING DIFFERENT R OBJECTS
## AND/OR FOR EXPORTING THESE WEATHER DATA.
## HERE ARE A FEW EXAMPLES
#########################################################
#########################################################
## The data array may be saved directly as an R object ##
save(wx.extent, file='wx_extent.Rdata')
## And then later recalled ##
load(file='wx_extent.Rdata')
################################################################
################################################################
## Another option is to create a raster object from the array ##
## For more info, see package raster
library(raster)
## Using the data from a query above ##
## First create a stacked raster object using the first
## layer (i.e. datetime) of the weather data array ##
## Notice the offset of 1.25 degrees (1/2 the spatial resolution)
## to describe the limits of the bounding box not the points
ras <- stack(raster(wx.extent1[,,1], crs="+proj=longlat +datum=WGS84",
xmn=min(as.numeric(dimnames(wx.extent1)[[2]])) - 1.25,
xmx=max(as.numeric(dimnames(wx.extent1)[[2]])) + 1.25,
ymn=min(as.numeric(dimnames(wx.extent1)[[1]])) - 1.25,
ymx=max(as.numeric(dimnames(wx.extent1)[[1]])) + 1.25))
## Then add each subsequent layer to the raster stack ##
for(i in 2:length(dimnames(wx.extent1)[[3]])){
ras <- addLayer(ras, raster(wx.extent1[,,i],
crs="+proj=longlat +datum=WGS84",
xmn=min(as.numeric(dimnames(wx.extent1)[[2]])) - 1.25,
xmx=max(as.numeric(dimnames(wx.extent1)[[2]])) + 1.25,
ymn=min(as.numeric(dimnames(wx.extent1)[[1]])) - 1.25,
ymx=max(as.numeric(dimnames(wx.extent1)[[1]])) + 1.25))
}
#######################################################
## Optionally, export a layer from the raster stack to
## a format that can be imported by Esri's ArcGIS products ##
## First, select a layer from the raster stack ##
ras1 <- raster(ras, layer=1)
## Then write the data from that layer to a .bil file ##
writeRaster(ras, filename='ras_example.bil', format="EHdr")
## The file will be saved in your current working directory ##
## The resulting file can be imported into ArcGIS ##
## by using the "Raster to Other Format" tool in the
## "To Raster" section of the ArcToolbox.
###########################################################
## NOTE: Weather data obtained from a Gaussian grid must
## first be resampled onto a regular grid !!! ##
## Here we use the interp.loess() function
## from the tgp package
## Using data from a T62 Gaussian grid queried above
## Interpolate the data from the first layer (i.e. datetime)
## onto a regular grid ##
library(tgp)
wx.reg <- interp.loess(x=rep(as.numeric(dimnames(wx.extent3)[[2]]),
each=length(dimnames(wx.extent3)[[1]])),
y=rep(as.numeric(dimnames(wx.extent3)[[1]]),
length(dimnames(wx.extent3)[[2]])),
z=as.vector(wx.extent3[,,1]), span=0.6,
gridlen=c(length(dimnames(wx.extent3)[[2]]),
length(dimnames(wx.extent3)[[1]])))
## Create a stacked raster object from the first layer
## (i.e. datetime) after interpolation ##
## Again, notice the offset (1/2 the resolution) ##
## Also note that the matrix (i.e. wx.reg$z) must be flipped
## along the y axis and transposed ##
## This is required b/c of the interpolation performed above ##
ras <- stack(raster(t(wx.reg$z[,length(wx.reg$y):1]),
crs="+proj=longlat +datum=WGS84",
xmn=min(as.numeric(wx.reg$x)) - abs(diff(wx.reg$x)[1]/2),
xmx=max(as.numeric(wx.reg$x)) + abs(diff(wx.reg$x)[1]/2),
ymn=min(as.numeric(wx.reg$y)) - abs(diff(wx.reg$y)[1]/2),
ymx=max(as.numeric(wx.reg$y)) + abs(diff(wx.reg$y)[1]/2)))
## Add each subsequent layer in the array to the raster stack ##
## after interpolating onto a regular grid ##
for(i in 2:length(dimnames(wx.extent3)[[3]])){
## Interpolate
t.wx.reg <- interp.loess(x=rep(as.numeric(dimnames(wx.extent3)[[2]]),
each=length(dimnames(wx.extent3)[[1]])),
y=rep(as.numeric(dimnames(wx.extent3)[[1]]),
length(dimnames(wx.extent3)[[2]])),
z=as.vector(wx.extent3[,,i]), span=0.6,
gridlen=c(length(dimnames(wx.extent3)[[2]]),
length(dimnames(wx.extent3)[[1]])))
## Note the offset ##
## Note flipping the matrix along the y axis and transposing ##
## Add layer to stack
ras <- addLayer(ras, raster(t(t.wx.reg$z[,length(t.wx.reg$y):1]),
crs="+proj=longlat +datum=WGS84",
xmn=min(as.numeric(t.wx.reg$x)) - abs(diff(t.wx.reg$x)[1]/2),
xmx=max(as.numeric(t.wx.reg$x)) + abs(diff(t.wx.reg$x)[1]/2),
ymn=min(as.numeric(t.wx.reg$y)) - abs(diff(t.wx.reg$y)[1]/2),
ymx=max(as.numeric(t.wx.reg$y)) + abs(diff(t.wx.reg$y)[1]/2)))
}
##################################################################
##################################################################
## Another option is to create a Spatial object
## using the sp package
## Again, data from a Gaussian grid may require special attention
## as the grid points are unevenly spaced
library(sp)
## Using the data from a query above
## Convert the array to a data.frame ##
wx.df <- NCEP.array2df(wx.extent2)
## 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')
#########################################################
## A Spatial object of a single datetime (i.e. layer) can
## be written to .asc, a format that may then be
## imported into ArcGIS.
## First, convert the first layer of the array to a data.frame ##
wx.df <- NCEP.array2df(wx.extent2[,,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. ##
##################################################################
## There are still more options for writing these data to files ##
## See e.g. writeMat() in the R.matlab package for writing Matlab files
## Also see the RSAGA package for GIS functionality in R
## One could even write the data array back to NetCDF
## (see packages RNetCDF and ncdf)
}
Run the code above in your browser using DataLab