r - このコードをどのように高速化しますか? netcdf からのデータの解像度を下げてから、統計用に xyz 形式に変換する

Question

トウモロコシの収量と収穫面積の netcdf を取得し、解像度を 2.5 分角から 0.5 度に縮小してから、全体を XYZ 形式に変換して、データとより簡単に「会話」できるようにします。その形式になっています。(他のデータを行列形式に変換できると思いますが、xyz が好きです。)

データはこちら。

以下のコードは、収穫された面積と平均収量から総生産量を計算するためのいくつかの関数を定義し、netcdf を照会するときにいくつかの「フィーダー」データを使用し、次に plyr を使用してフィーダーをループし、データを抽出し、関数を作成し、xyz に保存します。動作しますが、これらのファイルの 1 つだけを実行するのに約 30 分かかり、100 以上のファイルがあります。このコードを最適化する方法についての提案は素晴らしいものです。より大きなデータのチャンクを抽出してそれらに関数を適用する方が速いでしょうか? たぶん、地球のストライプ全体ですか？これがより速いかどうかをアプリオリに知るにはどうすればよいですか？

rm(list=ls()) 
library(ncdf)
library(reshape)
library(plyr)
library(sp)
library(maps)
library(rgeos)
library(maptools)
library(rworldmap)

getha = function(lat,size=lat[1]-lat[2]){
    lat1 = (lat-size/2)*pi/180
    lat2 = (lat+size/2)*pi/180
    lon1 = (0-size/2)*pi/180    #lon doesn't come in because all longitudes are great circles
    lon2 = (0+size/2)*pi/180
    6371^2 * abs(sin(lat1)-sin(lat2))*abs(lon1-lon2)*100    #6371 is the radius of the earth and 100 is the number of ha in a km^2
    }

gethamat = function(mat,latvec,blocksize=6){
    a = getha(latvec)
    areamat = matrix(rep(a,blocksize),blocksize)
    area = t(mat)*areamat   #The matrix is transposed because the dimensions of the Ramankutty's netcdf's are switched
    area
    }

getprod = function(yieldblock,areablock,latvec){
    b = gethamat(areablock,latvec)
    sum(t(yieldblock)*b,na.rm=TRUE)
    }

lat = as.matrix(seq(from=89.75,to=-89.75,by=-.5))
lon = as.matrix(seq(from=-179.75,to=179.75,by=.5))

lon = seq.int(from=1,to=4320,by=6)
lat = seq.int(from=1,to=2160,by=6)

lat = rep(lat,720)
lon = t(matrix(lon,720,360))
lon = as.data.frame(lon)
l = reshape(lon,direction="long",varying=list(colnames(lon)),v.names = "V")
coords = data.frame(cbind(l[,2],lat))
colnames(coords) = c("lng","lat")
feeder = coords
head(feeder)

maize.nc = open.ncdf('maize_5min.nc')

getcrops = function(feed,netcdf,var="cropdata"){
    yieldblock = get.var.ncdf(netcdf,varid=var,start = c(as.numeric(feed[1]),as.numeric(feed[2]),2,1),count = c(6,6,1,1))
    areablock = get.var.ncdf(netcdf,varid=var,start = c(as.numeric(feed[1]),as.numeric(feed[2]),1,1),count = c(6,6,1,1))
    lat = get.var.ncdf(netcdf,varid="latitude",start = feed[2],count = 6)
    prod = getprod(yieldblock,areablock,lat)
    lon = get.var.ncdf(netcdf,varid="longitude",start = feed[1],count = 6)
    #print(c(mean(lat),mean(lon)))
    data.frame(lat=mean(lat),lon = mean(lon),prod=prod)
    }

out = adply(as.matrix(feeder),1,getcrops,netcdf=maize.nc,.parallel=FALSE)

前もって感謝します。