Как бы вы ускорили этот код? Уменьшение разрешения данных из netcdf, а затем преобразование их в формат xyz для статистики

Question

Я беру netcdf урожая кукурузы и убранной площади, уменьшая разрешение с 2,5 угловых минут до.5 градусов, а затем преобразовываю все это в формат XYZ, чтобы мне было легче "разговаривать" с данными, которые я использую. мы получили в этом формате. (Я полагаю, что я мог бы превратить мои другие данные в матричную форму, но мне нравится xyz.)

Данные здесь.

Приведенный ниже код определяет несколько функций для расчета общего производства по убранной площади и средней урожайности, а затем он использует некоторые данные "фидера" при запросе netcdf, а затем использует plyr для циклического прохождения фидера, извлечения данных, применения функции, а затем сохранить в XYZ. Это работает, но для запуска только одного из этих файлов требуется около 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)

Заранее спасибо.

Answer 1

plyr функции, как известно, медленные, когда количество кусков увеличивается. Я действительно рекомендовал бы хранить данные в многомерном массиве. Это позволяет вам использовать apply например получить среднее для всех lat-lon комбинации и т. д. Многомерный массив требует меньше памяти RAM, поскольку метаданные хранятся не непосредственно в виде столбцов, а неявно в пределах измерений массива. К тому же, apply часто гораздо быстрее, чем при использовании plyr, ncdf Пакет изначально считывает данные в многомерные массивы, так что это также экономит вам шаг обработки (например, используя melt).

После сокращения набора данных я бы часто использовал melt перейти к тому, что вы называете XYZ формат для построения. Но к тому времени набор данных настолько мал, что это на самом деле не имеет значения.