{"id":833,"date":"2011-03-01T20:06:02","date_gmt":"2011-03-02T01:06:02","guid":{"rendered":"http:\/\/lukemiller.org\/?p=833"},"modified":"2015-11-30T16:00:20","modified_gmt":"2015-12-01T00:00:20","slug":"extracting-sea-surface-temperatures-from-noaas-oisstv2","status":"publish","type":"post","link":"https:\/\/lukemiller.org\/index.php\/2011\/03\/extracting-sea-surface-temperatures-from-noaas-oisstv2\/","title":{"rendered":"Extracting sea surface temperatures from NOAA’s OISSTv2"},"content":{"rendered":"

Update, 2015-11-30<\/strong> It appears that NOAA has gone through and upgraded all of the OISST files to the newer version of the NetCDF file format. As a result, the functions outlined in this post don’t work any longer. Instead, see the updated functions in my newer post, https:\/\/lukemiller.org\/index.php\/2014\/11\/extracting-noaa-sea-surface-temperatures-with-ncdf4\/<\/a>. The concepts are the same as described here, but the newer functions use the ncdf4 package to access the newer NetCDF file format. <\/em><\/p>\n

\"\"<\/a>
Sea surface temperatures<\/figcaption><\/figure>\n

NOAA’s Physical Sciences Division<\/a> produces a global map of weekly averaged sea surface temperatures, derived from satellite sources. The data can be found at http:\/\/www.esrl.noaa.gov\/psd\/data\/gridded\/data.noaa.oisst.v2.html<\/a>. I’ve put together a R script that will extract data from a region of interest and optionally plot it for you. The input data file is a 135+ MB netCDF file found at the link above.<\/p>\n

The OISSTv2 data are arranged in a 1\u00b0 x 1\u00b0 grid covering the entire world, and a weekly average has been calculated for each grid cell for each 7-day period starting Jan 1, 1990, in the file sst.wkmean.1990-present.nc<\/code>. You will also want to grab the land-sea mask file lsmask.nc<\/code> as well. The land-sea mask is necessary to remove the spurious temperature values that fall on the continents. In the script below I use the land-sea mask to simply convert all points over land into NA<\/code> values in the output matrix.<\/p>\n

This script makes use of two add-on packages available from the CRAN repositories. First, the ncdf<\/code><\/a> package, written by David Pierce and maintained by Brian Ripley, is used to open the netCDF file and extract the desired data. The fields<\/code><\/a> package, written by Reinhard Furrer, Douglas Nychka and Stephen Sain, provides the tools to plot the data with a color scale bar alongside the plot.<\/p>\n

The NOAA_OISST_extraction.R script can be downloaded here<\/a>.<\/del> Please see the newer functions described in https:\/\/lukemiller.org\/index.php\/2014\/11\/extracting-noaa-sea-surface-temperatures-with-ncdf4\/<\/a><\/p>\n

\n
\n
# Filename: NOAA_OISST_extraction.R<\/span>\r\n# <\/span>\r\n# Author: Luke Miller   Mar 1, 2011<\/span>\r\n###############################################################################<\/span>\r\nlibrary<\/span><\/a>(<\/span>ncdf<\/span><\/a>)<\/span>\r\n\r\n# For info on the NOAA Optimum Interpolated Sea Surface Temperature V2 (OISST):<\/span>\r\n# http:\/\/www.esrl.noaa.gov\/psd\/data\/gridded\/data.noaa.oisst.v2.html<\/span>\r\n\r\n# This script works on the OISST sea surface temperature netCDF file called <\/span>\r\n# sst.wkmean.1990-present.nc, which must be downloaded into the current <\/span>\r\n# working directory (file size is >135MB currently)<\/span>\r\n\r\n# The OISST grid layout is 1 degree latitudes from 89.5 to -89.5 degrees north.<\/span>\r\n# The longitude grid is 1 degree bins, 0.5 to 359.5 degrees east.<\/span>\r\n# The SST values are given in degrees Celsius.<\/span>\r\n# Time values are \"days since 1800-1-1 00:00\", starting from 69395<\/span>\r\n# Time delta_t is 0000-00-07 00:00, indices range from 0 to 1103 currently<\/span>\r\n# The time values in the sst field give the 1st day of the week for each <\/span>\r\n# weekly temperature value (though the data are meant to be centered on <\/span>\r\n# Wednesday of that week), starting at 69395 = 12\/31\/1989.<\/span>\r\n# Use as.Date(69395,origin = '1800-1-1') to convert the netCDF day value to a <\/span>\r\n# human readable form. <\/span>\r\n\r\nlats = seq<\/span><\/a>(<\/span>89.5<\/span>,<\/span>-<\/span>89.5<\/span>,<\/span>-<\/span>1<\/span>)<\/span> #generate set of grid cell latitudes (center of cell)<\/span>\r\nlons = seq<\/span><\/a>(<\/span>0.5<\/span>,<\/span>359.5<\/span>,<\/span>1<\/span>)<\/span> #generate set of grid cell longitudes (center of cell)<\/span>\r\n\r\n#Ask user to enter the boundaries of the search area<\/span>\r\ncat<\/span><\/a>(<\/span>\"Enter longitude of western edge of search area (degrees east 0 to 359)\\n<\/span>\"<\/span>)<\/span>\r\nlon1 = scan<\/span><\/a>(<\/span>what = numeric<\/span><\/a>(<\/span>0<\/span>)<\/span>,<\/span>n = 1<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>\"Enter longitude of eastern edge of search area (degrees east 0 to 359)\\n<\/span>\"<\/span>)<\/span>\r\nlon2 = scan<\/span><\/a>(<\/span>what = numeric<\/span><\/a>(<\/span>0<\/span>)<\/span>,<\/span>n = 1<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>\"Enter latitude of northern edge\\n<\/span>\"<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>\"of search area (degrees north, 89.5 to -89.5\\n<\/span>\"<\/span>)<\/span>\r\nlat1 = scan<\/span><\/a>(<\/span>what = numeric<\/span><\/a>(<\/span>0<\/span>)<\/span>,<\/span>n = 1<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>\"Enter latitude of southern edge\\n<\/span>\"<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>\"of search area (degrees north, 89.5 to -89.5\\n<\/span>\"<\/span>)<\/span>\r\nlat2 = scan<\/span><\/a>(<\/span>what = numeric<\/span><\/a>(<\/span>0<\/span>)<\/span>,<\/span>n = 1<\/span>)<\/span>\r\n\r\nlon1a = which.min<\/span><\/a>(<\/span>abs<\/span><\/a>(<\/span>lon1 -<\/span> lons)<\/span>)<\/span> #get index of nearest longitude value<\/span>\r\nlon2a = which.min<\/span><\/a>(<\/span>abs<\/span><\/a>(<\/span>lon2 -<\/span> lons)<\/span>)<\/span> #get index of nearest longitude value<\/span>\r\nlat1a = which.min<\/span><\/a>(<\/span>abs<\/span><\/a>(<\/span>lat1 -<\/span> lats)<\/span>)<\/span> #get index of nearest latitude value<\/span>\r\nlat2a = which.min<\/span><\/a>(<\/span>abs<\/span><\/a>(<\/span>lat2 -<\/span> lats)<\/span>)<\/span> #get index of nearest latitude value<\/span>\r\n#The lon\/lat 1a\/2a values should now correspond to indices in the netCDF file<\/span>\r\n#for the desired grid cell. <\/span>\r\nnlon = (<\/span>lon2a -<\/span> lon1a)<\/span> +<\/span> 1<\/span> #get number of longitudes to extract<\/span>\r\nnlat = (<\/span>lat2a -<\/span> lat1a)<\/span> +<\/span> 1<\/span> #get number of latitudes to extract<\/span>\r\n\r\n#Ask the user to enter the dates of interest<\/span>\r\ncat<\/span><\/a>(<\/span>\"Enter the starting date in the form: 1990-1-31\\n<\/span>\"<\/span>)<\/span>\r\ndate1 = scan<\/span><\/a>(<\/span>what = character<\/span><\/a>(<\/span>0<\/span>)<\/span>,<\/span>n = 1<\/span>)<\/span>\r\ncat<\/span><\/a>(<\/span>\"Enter the ending date in the form: 1990-1-31\\n<\/span>\"<\/span>)<\/span>\r\ndate2 = scan<\/span><\/a>(<\/span>what = character<\/span><\/a>(<\/span>0<\/span>)<\/span>,<\/span>n = 1<\/span>)<\/span>\r\n#date1 = '2007-08-01'<\/span>\r\n#date2 = '2007-09-01'<\/span>\r\ndate1 = as.Date<\/span><\/a>(<\/span>date1,<\/span> format<\/span><\/a> = \"%Y-%m-%d\"<\/span>)<\/span> #convert to Date object<\/span>\r\ndate2 = as.Date<\/span><\/a>(<\/span>date2,<\/span> format<\/span><\/a> = \"%Y-%m-%d\"<\/span>)<\/span> #convert to Date object<\/span>\r\n\r\n#Open the netCDF file for reading. <\/span>\r\nnc = open.ncdf(<\/span>\"sst.wkmean.1990-present.nc\"<\/span>)<\/span>\r\n#print.ncdf(nc) will show info about the structure of the netCDF file<\/span>\r\n\r\n#Extract available dates from netCDF file<\/span>\r\nncdates = nc$<\/span>dim$<\/span>time$<\/span>vals\r\nncdates = as.Date<\/span><\/a>(<\/span>ncdates,<\/span>origin = '1800-1-1'<\/span>)<\/span> #available time points in nc<\/span>\r\n\r\ndate1a = which.min<\/span><\/a>(<\/span>abs<\/span><\/a>(<\/span>date1 -<\/span> ncdates)<\/span>)<\/span> #get index of nearest time point<\/span>\r\ndate2a = which.min<\/span><\/a>(<\/span>abs<\/span><\/a>(<\/span>date2 -<\/span> ncdates)<\/span>)<\/span> #get index of nearest time point<\/span>\r\nndates = (<\/span>date2a -<\/span> date1a)<\/span> +<\/span> 1<\/span>#get number of time steps to extract<\/span>\r\n\r\n#Extract the data from the netCDF file to a matrix or array called 'sstout'. <\/span>\r\nsstout = get.var.ncdf(<\/span>nc,<\/span> varid = 'sst'<\/span>,<\/span> start<\/span><\/a> = c<\/span><\/a>(<\/span>lon1a,<\/span>lat1a,<\/span>date1a)<\/span>,<\/span>\r\n\t\tcount = c<\/span><\/a>(<\/span>nlon,<\/span>nlat,<\/span>ndates)<\/span>)<\/span>\r\n#If you only retrieve one time point, sstout will be a 2D matrix with <\/span>\r\n#longitudes running down the rows, and latitudes across the columns. <\/span>\r\n#For example, Column 1, sstout[1:nrow(sstout),1], will contain sst values for <\/span>\r\n#each of the longitude values at the northern-most latitude in your search <\/span>\r\n#area, with the first row, sstout[1,1], being the western-most longitude, and <\/span>\r\n#the last row being the eastern edge of your search area.<\/span>\r\n#If you retrieve multiple time points, sstout will be a 3D array, where time is<\/span>\r\n#the 3rd dimension. Lat\/lon will be arranged the same as the 2D case. <\/span>\r\n\r\n#The vector 'datesout' will hold the Date values associated with each set of <\/span>\r\n#sst values in the sstout array, should you need to access them.<\/span>\r\ndatesout = ncdates[<\/span>date1a:<\/span>date2a]<\/span>\r\n\r\n###############################################################################<\/span>\r\n###############################################################################<\/span>\r\n# The NOAA OISST files contain sea surface temperatures for the entire globe,<\/span>\r\n# including on the continents. This clearly isn't right, so they also supply a<\/span>\r\n# land-sea mask file in netCDF format at the website listed at the start of <\/span>\r\n# this script. <\/span>\r\n# We use the values (0 or 1) stored in the mask file to turn all of the <\/span>\r\n# continent areas into NA's. <\/span>\r\n\r\nnc2 = open.ncdf(<\/span>'lsmask.nc'<\/span>)<\/span> #open land-sea mask<\/span>\r\nmask = get.var.ncdf(<\/span>nc2,<\/span> varid = \"mask\"<\/span>,<\/span>start<\/span><\/a> = c<\/span><\/a>(<\/span>lon1a,<\/span>lat1a,<\/span>1<\/span>)<\/span>,<\/span>\r\n\t\tcount = c<\/span><\/a>(<\/span>nlon,<\/span>nlat,<\/span>1<\/span>)<\/span>)<\/span> #get land-sea mask values (0 or 1)<\/span>\r\n\r\nmask = ifelse<\/span><\/a>(<\/span>mask ==<\/span> 0<\/span>,<\/span>NA<\/span>,<\/span>1<\/span>)<\/span> #replace 0's with NA's<\/span>\r\n\r\nif<\/span> (<\/span>is.matrix<\/span><\/a>(<\/span>sstout)<\/span>)<\/span> {<\/span> #if there is only 1 time point (2D matrix)<\/span>\r\n\tsstout = sstout *<\/span> mask #all masked values become NA\t<\/span>\r\n}<\/span>\r\nif<\/span> (<\/span>!<\/span>is.matrix<\/span><\/a>(<\/span>sstout)<\/span>)<\/span> {<\/span> #if ssout is a 3D matrix<\/span>\r\n\tdims = dim<\/span><\/a>(<\/span>sstout)<\/span>\r\n\tfor<\/span> (<\/span>i in<\/span> 1<\/span>:<\/span>dims[<\/span>3<\/span>]<\/span>)<\/span>{<\/span>\r\n\t\tsstout[<\/span>,,<\/span>i]<\/span> = sstout[<\/span>,,<\/span>i]<\/span> *<\/span> mask #all masked values become NA<\/span>\r\n\t}<\/span>\r\n}<\/span>\r\n###############################################################################<\/span>\r\n###############################################################################<\/span>\r\n#Now we can produce a plot of the data. <\/span>\r\nlibrary<\/span><\/a>(<\/span>fields<\/span><\/a>)<\/span>\r\ngoodlons = lons[<\/span>lon1a:<\/span>lon2a]<\/span>\r\ngoodlats = lats[<\/span>lat1a:<\/span>lat2a]<\/span>\r\n\r\n#plot the 1st time point<\/span>\r\nif<\/span> (<\/span>is.matrix<\/span><\/a>(<\/span>sstout)<\/span>)<\/span>{<\/span> #if there is only one time point<\/span>\r\n\tpar<\/span><\/a>(<\/span>mar<\/span><\/a> = c<\/span><\/a>(<\/span>5<\/span>,<\/span>5<\/span>,<\/span>4<\/span>,<\/span>6<\/span>)<\/span>)<\/span> #widen margins<\/span>\r\n\timage<\/span><\/a>(<\/span>goodlons,<\/span>rev<\/span><\/a>(<\/span>goodlats)<\/span>,<\/span>sstout[<\/span>,<\/span>]<\/span>,<\/span>\r\n\t\t\tylim = rev<\/span><\/a>(<\/span>range<\/span><\/a>(<\/span>goodlats)<\/span>)<\/span>,<\/span> #reverse y-axis<\/span>\r\n\t\t\tcol<\/span><\/a> = tim.colors(<\/span>32<\/span>)<\/span>,<\/span>\r\n\t\t\tyaxt = \"n\"<\/span>,<\/span> xaxt = \"n\"<\/span>,<\/span>\r\n\t\t\tylab = \"Latitude (degrees north)\"<\/span>,<\/span>\r\n\t\t\txlab = \"Longitude (degrees east)\"<\/span>,<\/span>\r\n\t\t\tmain = paste<\/span><\/a>(<\/span>\"Week of \"<\/span>,<\/span>datesout[<\/span>1<\/span>]<\/span>,<\/span>sep = \"\"<\/span>)<\/span>)<\/span>\r\n\taxis<\/span><\/a>(<\/span>1<\/span>,<\/span>at = pretty<\/span><\/a>(<\/span>goodlons)<\/span>,<\/span>labels<\/span><\/a> = pretty<\/span><\/a>(<\/span>goodlons)<\/span>,<\/span> las = 1<\/span>)<\/span>\r\n\taxis<\/span><\/a>(<\/span>2<\/span>,<\/span>at = rev<\/span><\/a>(<\/span>pretty<\/span><\/a>(<\/span>goodlats)<\/span>)<\/span>,<\/span>labels<\/span><\/a> = pretty<\/span><\/a>(<\/span>goodlats)<\/span>,<\/span> las = 1<\/span>)<\/span>\r\n\r\n\t#image.plot from the fields package inserts a color bar<\/span>\r\n\timage.plot(<\/span>zlim = range<\/span><\/a>(<\/span>sstout[<\/span>,<\/span>]<\/span>,<\/span>na.rm = TRUE<\/span>)<\/span>,<\/span>nlevel = 32<\/span>,<\/span>\r\n\t\t\tlegend.only = TRUE<\/span>,<\/span> horizontal = FALSE<\/span>,<\/span> col<\/span><\/a> = tim.colors(<\/span>32<\/span>)<\/span>,<\/span>\r\n\t\t\tlegend.args = list<\/span><\/a>(<\/span>text<\/span><\/a> = \"Temperature, C\"<\/span>,<\/span> cex = 1.2<\/span>,<\/span>\r\n\t\t\t\t\tside = 4<\/span>,<\/span> line<\/span><\/a> = 2<\/span>)<\/span>)<\/span>\r\n}<\/span> else<\/span> {<\/span>\r\n\t#If there is more than one time point, just plot the first time point<\/span>\r\n\t<\/a>par<\/span><\/a>(<\/span>mar<\/span><\/a> = c<\/span><\/a>(<\/span>4.5<\/span>,<\/span>4.5<\/span>,<\/span>3<\/span>,<\/span>6<\/span>)<\/span>)<\/span> #widen margins<\/span>\r\n\timage<\/span><\/a>(<\/span>goodlons,<\/span>rev<\/span><\/a>(<\/span>goodlats)<\/span>,<\/span>sstout[<\/span>,,<\/span>1<\/span>]<\/span>,<\/span>\r\n\t\t\tylim = rev<\/span><\/a>(<\/span>range<\/span><\/a>(<\/span>goodlats)<\/span>)<\/span>,<\/span> #reverse y-axis<\/span>\r\n\t\t\tcol<\/span><\/a> = tim.colors(<\/span>32<\/span>)<\/span>,<\/span>\r\n\t\t\tyaxt = \"n\"<\/span>,<\/span> xaxt = \"n\"<\/span>,<\/span>\r\n\t\t\tylab = \"Latitude (degrees north)\"<\/span>,<\/span>\r\n\t\t\txlab = \"Longitude (degrees east)\"<\/span>,<\/span>\r\n\t\t\tmain = paste<\/span><\/a>(<\/span>\"Week of \"<\/span>,<\/span>datesout[<\/span>1<\/span>]<\/span>,<\/span>sep = \"\"<\/span>)<\/span>)<\/span>\r\n\taxis<\/span><\/a>(<\/span>1<\/span>,<\/span>at = pretty<\/span><\/a>(<\/span>goodlons)<\/span>,<\/span>labels<\/span><\/a> = pretty<\/span><\/a>(<\/span>goodlons)<\/span>,<\/span> las = 1<\/span>)<\/span>\r\n\taxis<\/span><\/a>(<\/span>2<\/span>,<\/span>at = rev<\/span><\/a>(<\/span>pretty<\/span><\/a>(<\/span>goodlats)<\/span>)<\/span>,<\/span>labels<\/span><\/a> = pretty<\/span><\/a>(<\/span>goodlats)<\/span>,<\/span> las = 1<\/span>)<\/span>\r\n\t#image.plot from the fields package inserts a color bar<\/span>\r\n\timage.plot(<\/span>zlim = range<\/span><\/a>(<\/span>sstout[<\/span>,,<\/span>1<\/span>]<\/span>,<\/span>na.rm = TRUE<\/span>)<\/span>,<\/span>nlevel = 32<\/span>,<\/span>\r\n\t\t\tlegend.only = TRUE<\/span>,<\/span> horizontal = FALSE<\/span>,<\/span> col<\/span><\/a> = tim.colors(<\/span>32<\/span>)<\/span>,<\/span>\r\n\t\t\tlegend.args = list<\/span><\/a>(<\/span>text<\/span><\/a> = \"Temperature, C\"<\/span>,<\/span> cex = 1.2<\/span>,<\/span>\r\n\t\t\t\t\tside = 4<\/span>,<\/span> line<\/span><\/a> = 2<\/span>)<\/span>)<\/span>\r\n}<\/span><\/pre>\n<\/div>\n<\/div>\n
Created by Pretty R at inside-R.org<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"
Update, 2015-11-30 It appears that NOAA has gone through and upgraded all of the OISST files to the newer version of the NetCDF file format. As a result, the functions outlined in this post don’t work any longer. Instead, see the updated functions in my newer post, https:\/\/lukemiller.org\/index.php\/2014\/11\/extracting-noaa-sea-surface-temperatures-with-ncdf4\/. The concepts are the same as described […]<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"closed","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[218],"tags":[76,74,75,64,77,58,73],"class_list":["post-833","post","type-post","status-publish","format-standard","hentry","category-r-project","tag-fields","tag-ncdf","tag-netcdf","tag-noaa","tag-oisst","tag-r-project","tag-sea-surface-temperature"],"_links":{"self":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts\/833","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/comments?post=833"}],"version-history":[{"count":13,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts\/833\/revisions"}],"predecessor-version":[{"id":2143,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/posts\/833\/revisions\/2143"}],"wp:attachment":[{"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/media?parent=833"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/categories?post=833"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/lukemiller.org\/index.php\/wp-json\/wp\/v2\/tags?post=833"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}