writing a loop for upscaling precipitation for USA

Question

I am writing a code to calculate the mean amount of precipitation for different regions of conterminous USA. My total data has 300 times 120 (lon*lat) grids in Netcdf format. I want to write a loop in R to take the average of each 10 by 10 number of grids and assign that value (average) to all of the grids inside the region and repeat this for the next region. At the end instead of a 120 by 300 grids I will have 12 by 30 grids. So this is kind a upscaling method I want to apply to my data. I can use a for-loop for each region separately but It makes my code very huge and I don’t want to do that. Any idea would be appreciated. Thanks. P.S: Here is the function I have written for one region (10by10) lat*lon.

upscaling <- function(file, variable, start.time=1, count.time=1)
{
  library(ncdf)   # load ncdf library to manipulate ncdf data
  ncdata <- open.ncdf(file);      # open ncdf file
  lon  <- get.var.ncdf(ncdata, "lon");
  lat  <- get.var.ncdf(ncdata, "lat");
  time <- get.var.ncdf(ncdata, "time");
  start.lon <- 1
  end.lon   <- length(lon)
  start.lat <- 1
  end.lat   <- length(lat)
  count.lon <- end.lon - start.lon + 1;   # count number of longitude
  count.lat <- end.lat - start.lat + 1;   # count number of latitude
  dat <- get.var.ncdf(ncdata, variable, start=c(start.lon, start.lat, 1),               
                      count=c(count.lon, count.lat, 1))
  temp.data<- array(0,dim=c(10,10))
  for (i in 1:10)
  {
    for (j in 1:10)
    {
      temp.data <- mean(dat[i,j,])
    }
  }
}

Answer 1

There is no need to make a messy loop to spatially aggregate your data. Just use the aggregate function in the raster package:

library(raster)
a=matrix(data=c(1:100),nrow=10,ncol=10)
a=raster(a)
ra  <-  aggregate(a,  fact=5,  fun=mean) #fact=5 will aggregate using a 5x5 window
ra=as.matrix(ra)
ra

Now for your netcdf data, use raster's rasterFromXYZ to create the raster that can then be aggregated with the above method. Bonus includes the option to define your projection as an argument in the function so you end up with a georeferenced object at the end. This is important because if you aggregate your data without it you will then have to figure out by hand how to georeference the resulting matrix.

EDIT: If you want a resulting raster with the same dimensions as the original one, disaggregate the data right after aggregating it. While this seems redundant, these raster methods are very fast.

library(raster)
a=matrix(data=c(1:100),nrow=10,ncol=10)
a=raster(a)
ra  <-  aggregate(a,  fact=5,  fun=mean) #fact=5 will aggregate using a 5x5 window
ra  <-  disaggregate(ra, fact=5)
ra=as.matrix(ra)
ra

Answer 2

If you grid definitions follow standard netcdf conventions, then you might be able to remap using the CDO remapping functions. For first order conservative remapping you can try

cdo remapcon,grid_specification_here in.nc out.nc 

Note that the answer given above is approximate, and not quite correct as the grid cell size is not the same as a function of latitude. The size of the error is likely small for this particular task as the cell sizes are fine, but nevertheless the answer will be slightly off.