0

我正在获取玉米产量和收获面积的 netcdf，将分辨率从 2.5 弧分缩小到 0.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)

提前致谢。

4

1 回答 1

2

plyr当块的数量变大时，函数的速度非常慢。我真的建议将数据保存在多维数组中。这允许您使用apply例如获取所有lat-lon组合的平均值等。多维数组占用更少的 RAM 存储空间，因为元数据不是直接存储为列，而是隐含在数组的维度中。此外，apply通常比使用plyr. 该ncdf包本机将数据读取到多维数组中，因此这也为您节省了一个处理步骤（例如使用melt）。

减少数据集后，我经常会使用melt您所说的 XYZ 格式进行绘图。但是到那时，数据集太小了，这并不重要。

于 2013-04-03T09:22:52.027 回答