UPDATE
THIS POST IS DEPRECATED. SEE /2015/07/climate-grids-and-thredds-server-experimenting-update
Abstract
- Climate grids have become readily available via netCDF and THREDDS
Methods
- Use the R package ncdf, read data by changing the URL passed to the ncdf reader
- I wanted to convert the matrix to a spatial raster object but struggled with the ‘flipped’ orientation of the input I got.
- Found this link [[http://stackoverflow.com/a/137]] which says “The reason is that the NetCDF interface you are using is very low-level, and all you have done is read out the variable without any of its dimension information. The orientation of the grid is really arbitrary, and the coordinate information needs to be understood in a particular context” AND SO:
- first transpose the matrix
- r<-raster(t(vals), …
- then in another test of a DIFFERENT grid product I had to use the flip tool
- d <- flip(r, direction = “y”)
- That flipped around “y”, keeping the georeferencing from the original context.
Results
eMAST Grids test
# ref http://www.emast.org.au/observations/climate/
#install.packages("ncdf", type = "source", configure.args="--with-netcdf-include=/usr/include")
require(ncdf)
## Loading required package: ncdf
#install.packages("raster")
require(raster)
## Loading required package: raster
## Loading required package: sp
# install.packages("rgdal")
require(rgdal)
## Loading required package: rgdal
## rgdal: version: 0.9-1, (SVN revision 518)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
## Path to GDAL shared files: /usr/share/gdal
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
# if extracting for points shapefile
#shp <- readOGR(dsn="test.shp", layer='test')
#plot(shp, add = T)
# a loop through days, see comment sections that print for debugging
strt <-'2012-01-01'
end <- '2012-01-04'
dates <- seq(as.Date(strt),as.Date(end),1)
dates
## [1] "2012-01-01" "2012-01-02" "2012-01-03" "2012-01-04"
# if extracting to shp then set up an output dataframe to collect
#dat_out <- as.data.frame(matrix(nrow = 0, ncol = 4))
# else just plots
par(mfrow = c(2,2))
for(i in 1:length(dates)){
# i=1
date_i <- dates[i]
infile <- sprintf("http://dapds00.nci.org.au/thredds/dodsC/rr9/Climate/eMAST/ANUClimate/0_01deg/v1m0_aus/day/land/tmin/e_01/2012/eMAST_ANUClimate_day_tmin_v1m0_%s.nc", gsub("-", "", date_i))
nc <- open.ncdf(infile)
vals <- get.var.ncdf(nc, varid="air_temperature")
nc.att <- nc$var$air_temperature
xmin <- min(nc.att$dim[[1]]$vals)
xmax <- max(nc.att$dim[[1]]$vals)
ymin <- min(nc.att$dim[[2]]$vals)
ymax <- max(nc.att$dim[[2]]$vals)
print(c(xmin,xmax))
print(c(ymin,ymax))
r <- raster(t(vals),
xmn=xmin, xmx=xmax,
ymn=ymin, ymx=ymax)
#str(r)
plot(r)
title(date_i)
#image(r)
#e <- extract(r, shp, df=T)
#str(e)
#e1 <- shp@data
#e1$date <- date_i
#e1$values <- e[,2]
#dat_out <- rbind(dat_out, as.data.frame(e1))
}
## [1] 112.905 153.995
## [1] -43.735 -9.005
## [1] 112.905 153.995
## [1] -43.735 -9.005
## [1] 112.905 153.995
## [1] -43.735 -9.005
## [1] 112.905 153.995
## [1] -43.735 -9.005
#dat_out
sessionInfo()
## R version 3.0.1 (2013-05-16)
## Platform: x86_64-redhat-linux-gnu (64-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=C LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rgdal_0.9-1 raster_2.3-12 sp_1.0-16 ncdf_1.6.8 knitr_1.8
##
## loaded via a namespace (and not attached):
## [1] evaluate_0.5.5 formatR_1.0 grid_3.0.1 lattice_0.20-15
## [5] stringr_0.6.2 tools_3.0.1
