braz

Monday, December 28, 2009

R for Global Wind Wave / Wave Swell Forecasting

As a follow up to using freely available existing climate datasets, here's one that might be of interest to those surfers and sailors out there. It is similar approach to that used for the temperature forecasting so go read that post first, it will help.

The free information from the U.S. NOAA/National Weather Service's National Centers for Environmental Prediction, in particular its Environmental Modeling Center and its WaveWatch dataset. It is stored as a GRIB2 files and this needs to be pre-processed using wgrib2 (available on MacPorts) to the netcdf format. This stage is similar to the previous post.

There were a few little glitches with this dataset that are immediately visible and the fact there is a big white patch in the Atlantic is due to the maximum height being stored as 9.999e+20. This normally won't be a problem but unfortunately the data also stores dry land with this value. A workaround would be to know or compute longitudes and latitude such that we could tell if a given position was over land or over sea prior to processing it. Unfortunately as this is a bit of a quick and dirty hack, I didn't have the time so here is the basic approach warts and all.

# See http://www.nco.ncep.noaa.gov/pmb/products/wave/nww3.t00z.grib.grib2.shtml
# Data was from ftp ftpprd.ncep.noaa.gov dir: /pub/data/nccf/com/wave/prod/wave.20091227
# wgrib2 was installed using macports
# /opt/local/bin/wgrib2 -s nww3.t00z.grib.grib2 | grep "HTSGW:surface:an" | /opt/local/bin/wgrib2 -i nww3.t00z.grib.grib2 -netcdf WAVE.nc

library(ncdf)
waveFrac <- open.ncdf("WAVE.nc")
wave <- get.var.ncdf(waveFrac, "HTSGW_surface")
# Dirty hack to fix input model, look for another better solution
wave[wave>9.99999]<- -1
x <- get.var.ncdf(waveFrac, "longitude")
y <- get.var.ncdf(waveFrac, "latitude")
library(fields)
rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")
image.plot(x,y,wave,col=rgb.palette(255),axes=F,main=as.expression(paste("Significant Height of Combined Wind Waves and Swell in Meters 2009-12-27", sep="")),legend.lab="Meters")
# Add a rough outline for islands, countries, and continents
contour(x,y,wave,add=TRUE,lwd=0.25,levels=0.2,drawlabels=FALSE,col="grey30")
# Add the source of the file and ftp location
text(130,-75,"Data source ftpprd.ncep.noaa.gov pub/data/nccf/com/wave/prod/wave.20091227")


NWW3WindWaves2009122700UTC

Labels: , , , ,

0 Comments:

Post a Comment

<< Home