Agustin.Lobo@ictja.csic.es
20140428
file: /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/vectorplotv2:log.R
Variables that describe local displacement, such as wind fields, can be very effectively represented with function rasterVis:vectorplot() by
Oscar Perpinan
This material does not substitute the standard documentation provided with the package and the documents provided by Oscar Perpinan, i.e
require(utils)
require(colorRamps)
require(RNetCDF)
## Loading required package: RNetCDF
require(rasterVis)
We will use a climatological wind field dowloaded from the NOAA Blended Sea Winds
Details of the example file:
# system('wget -q
# ftp://eclipse.ncdc.noaa.gov/pub/seawinds/SI/uv/clm/uvclm95to05.nc')
download.file("ftp://eclipse.ncdc.noaa.gov/pub/seawinds/SI/uv/clm/uvclm95to05.nc",
"uvclm95to05.nc", method = "curl")
Most climate data are in NetCDF format. There are several packages in R for dealing with this format, we will be using RNetCDF see RNetCDF – A Package for Reading and Writing NetCDF Datasets
mincwind <- open.nc("uvclm95to05.nc")
print.nc(mincwind)
## dimensions:
## zlev = 1 ;
## time = 12 ;
## lat = 719 ;
## lon = 1440 ;
## variables:
## int time(time) ;
## time:long_name = "Center Time of the Data" ;
## time:units = "days since 0000-01-01 00:00:00" ;
## float zlev(zlev) ;
## zlev:long_name = "height above sea level" ;
## zlev:units = "meters" ;
## float lat(lat) ;
## lat:long_name = "latitude" ;
## lat:units = "degrees_north" ;
## lat:grids = "uniform grids from -89.75 to 89.75 by 0.25" ;
## float lon(lon) ;
## lon:long_name = "longitude" ;
## lon:units = "degrees_east" ;
## lon:grids = "uniform grids from 0.00 to 359.75 by 0.25" ;
## float u(lon, lat, zlev, time) ;
## u:long_name = "Sea Surface Wind: x-component" ;
## u:units = "m/s" ;
## u:_FillValue = -9999 ;
## float v(lon, lat, zlev, time) ;
## v:long_name = "Sea Surface Wind: y-component" ;
## v:units = "m/s" ;
## v:_FillValue = -9999 ;
## float w(lon, lat, zlev, time) ;
## w:long_name = "Sea Surface Wind: magnitude as scalar mean" ;
## w:units = "m/s" ;
## w:_FillValue = -9999 ;
##
## // global attributes:
## :Conventions = "COARDS, CF-1.0, Unidata Dataset Discovery v1.0" ;
## :title = "NOAA Multiple-Satellite Blended 0.25-degree Sea Winds - climatological monthly means" ;
## :source = "Multiple satellite observations: DMSP SSMI F08, F10, F11, F13,F14 F15; TMI; QuikSCAT; AMSR-E; Direction from NCEP Reanalysis-2" ;
## :summary = "Gridded and blended sea surface vector winds from multiple satellites with direction from NCEP Reanalysis-2; Global ocean coverage with a 0.25-degree resolution; The whole datasets covers from July 1987 to present. Climatological (1995-2005) monthlies in this dataset; 6-hourly, daily & monthly resolutions are also available in other directories; See http://www.ncdc.noaa.gov/oa/rsad/blendedseawinds.html and links within for details. Include (u,v) means and scalar mean speed w for comparison" ;
## :Keywords = "sea winds, ocean winds, sea surface winds, air-sea interaction, air-sea flux, wind-driven circulation, Ekman pumping, Ekman transport, ocean upwelling, wind stress, windstress" ;
## :references = "links at http://www.ncdc.noaa.gov/oa/rsad/blendedseawinds.html" ;
## :History = "Simple spatiotemporally weighted Interpolation (SI), V.1.2. Version 1.2 uses updated satellite retrievals by Remote Sensing System, released in September 2006: SSMI V06, TMI V04, QSCAT V03a. AMSRE V05 was also updated using the new SSMI rain rate" ;
## :institution = "NOAA NESDIS National Climatic Data Center" ;
## :Contact = "Huai-Min.Zhang AT noaa.gov or ncdc.satorder AT noaa.gov; ph:1+828-271-4090 or 271-4800" ;
## :Acknowledgment = "The gridded data were generated from the multiple satellite observations of DOD, NOAA and NASA (and future others) and wind retrievals of the Remote Sensing Systems, Inc. (http://www.remss.com), using scientific methods such as objective analysis (OA). The OA is only truly objective when the needed statistics are completely known, which may not be always the case." ;
## :Comment = "The time in the netCDF should be the center of a calendar month, and chosen as Day 15 for all months for simplicity. The climatolog period is 1995 - 2005" ;
Wind vectors come as two fields with the u and v component for each cell
u <- var.get.nc(mincwind, "u")
class(u)
## [1] "array"
dim(u)
## [1] 1440 719 12
u[u == -9999] <- NA
v <- var.get.nc(mincwind, "v")
dim(v)
## [1] 1440 719 12
v[v == -9999] <- NA
Each u and v component has been recorded at 12 times (months), we get just the one of September (position 9) for our example.
Also, note that the arrays have 1440 rows and 719 columns, thus the planisphere is rotated 90 deg.
plot(raster(u[, , 9]))
plot(raster(t(u[, , 9])[ncol(u):1, ]))
u9 <- raster(t(u[, , 9])[ncol(u):1, ])
v9 <- raster(t(v[, , 9])[ncol(v):1, ])
w <- brick(u9, v9)
Now we need the coordinates:
wlon <- var.get.nc(mincwind, "lon")
wlat <- var.get.nc(mincwind, "lat")
range(wlon)
## [1] 0.0 359.8
This seems odd (Longitude growing from Greenwich to 359.75, instead of the common -180 - 180 range) but the map is centered on the Pacific.
projection(w) <- CRS("+init=epsg:4326")
extent(w) <- c(min(wlon), max(wlon), min(wlat), max(wlat))
w
## class : RasterBrick
## dimensions : 719, 1440, 1035360, 2 (nrow, ncol, ncell, nlayers)
## resolution : 0.2498, 0.2497 (x, y)
## extent : 0, 359.8, -89.75, 89.75 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer.1, layer.2
## min values : -22.01, -24.65
## max values : 22.70, 27.06
plot(w[[1]])
plot(w[[2]])
vectorplot(w * 10, isField = "dXY", region = FALSE, margin = FALSE, narrows = 10000)
slope <- sqrt(w[[1]]^2 + w[[2]]^2)
aspect <- atan2(w[[1]], w[[2]])
vectorplot(w * 10, isField = "dXY", region = slope, margin = FALSE, par.settings = rasterTheme(region = matlab.like(n = 10)),
narrows = 10000, at = 0:10)
vectorplot(stack(slope * 10, aspect), isField = TRUE, region = FALSE, margin = FALSE)
download.file("http://biogeo.ucdavis.edu/data/world/countries_shp.zip", "countries_shp.zip",
method = "curl")
unzip("countries_shp.zip", exdir = "countries_shp")
Read the shapefile as SpatialPolygonsDataFrame using rgdal
cntry <- readOGR(dsn = "./countries_shp", layer = "countries", stringsAsFactors = TRUE)
## OGR data source with driver: ESRI Shapefile
## Source: "./countries_shp", layer: "countries"
## with 265 features and 18 fields
## Feature type: wkbPolygon with 2 dimensions
projection(cntry)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
getClass(cntry)
## class : SpatialPolygonsDataFrame
## features : 265
## extent : -180, 180, -90, 83.63 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## variables : 18
## names : OBJECTID, NAME, ISO3, ISO2, FIPS, COUNTRY, ENGLISH, FRENCH, SPANISH, LOCAL, FAO, WAS_ISO, SOVEREIGN, CONTINENT, UNREG1, ...
## min values : 1, Afghanistan, ABW, AD, AA, Afghanistan, Afghanistan, Afghanistan, Afganistán, Afghanestan, Afghanistan, ANT, Afghanistan, Africa, Antartica, ...
## max values : 99, Zimbabwe, ZWE, ZW, ZI, Zimbabwe, Zimbabwe, Zimbabwe, Zimbabwe, Zimbabwe, Zimbabwe, YUG, Zimbabwe, South America, Western Europe, ...
vectorplot(w * 10, isField = "dXY", region = slope, margin = FALSE, par.settings = rasterTheme(region = matlab.like(n = 10)),
narrows = 10000, at = 0:10) + layer(sp.polygons(cntry))
oops! This is what happens when people do not conform to geographic standards.
Let's go to the standard representation with Greenwich Meridian at the center:
We first make the E and W parts:
we <- crop(w, extent(c(0, (720 * res(w)[1]), min(wlat), max(wlat))))
ww <- crop(w, extent(c((720 * res(w)[1]), 1440 * res(w)[1], min(wlat), max(wlat))))
extent(ww) <- c(-720 * res(w)[1], 0, min(wlat), max(wlat))
plot(ww[[1]])
plot(we[[1]])
extent(ww)
## class : Extent
## xmin : -179.9
## xmax : 0
## ymin : -89.75
## ymax : 89.75
extent(we)
## class : Extent
## xmin : 0
## xmax : 179.9
## ymin : -89.75
## ymax : 89.75
Now we merge
w2 <- merge(ww, we)
plot(w2[[1]])
And the overlay of the boundaries is now correct:
slope2 <- sqrt(w2[[1]]^2 + w2[[2]]^2)
vectorplot(w2 * 10, isField = "dXY", region = slope2, margin = FALSE, par.settings = rasterTheme(region = matlab.like(n = 10)),
narrows = 10000, at = 0:10) + layer(sp.polygons(cntry))
