Как изменить разрешение (или регрид) данных в R

У меня есть набор данных, состоящий из долготы, широты и среднемесячной переменной (например, температуры или осадков) за период с 1961 по 1970 год. Набор данных имеет разрешение 0,5 на 0,5 градуса долготы / широту и охватывает весь земной шар и был загружен как файл. Файл NC, данные из которого я извлек в R, используя:

library(ncdf)
f <- open.ncdf("D:/CRU/cru_ts3.21.1961.1970.tmp.dat.nc")
A <- get.var.ncdf(nc=f,varid="tmp")
B <- get.var.ncdf(nc=f,varid="lon")
C <- get.var.ncdf(nc=f,varid="lat")
D <- cbind(expand.grid(B, C))
E <- expand.grid(A)

Расширенная сетка (E) - это таблица данных, состоящая из 31 104 000 строк переменной, а расширенная сетка (D) - это таблица данных, состоящая из 259 200 строк долготы / широты. Если вы умножите 259 200 на 10 лет на 12 месяцев, вы получите 31 104 000. Следовательно, таблица E может быть разбита на месячные значения, используя:

Month <- 1
Start <- (Month-1)*(259200)+1
Finish <- (Month*259200)
G <- E[Start:Finish,]
H <- expand.grid(G)
I <- cbind(D,H) 

Таким образом, I теперь представляет собой таблицу данных за первый месяц (т.е. январь 1961 г.), состоящую из долготы, широты и переменной. Пример данных приведен ниже:

        lon    lat tmp
49184 -68.25 -55.75 7.5
49185 -67.75 -55.75 7.6
49186 -67.25 -55.75 7.6
49899 -70.75 -55.25 6.8
49900 -70.25 -55.25 7.0
49901 -69.75 -55.25 6.9
49902 -69.25 -55.25 7.1
49903 -68.75 -55.25 6.8
49904 -68.25 -55.25 7.6
49905 -67.75 -55.25 8.2

А теперь вопрос. Текущее разрешение сетки составляет 0,5 * 0,5 градуса, и я хотел бы «пересчитать» данные, чтобы разрешение было 0,25 * 0,25 градуса. Я не хочу делать ничего особенно умного с данными, поэтому я просто хочу, чтобы сетка 0,25 принимала значение сетки 0,5, в которой она находится, т.е. каждая сетка 0,5 * 0,5 содержит 4 сетки 0,25 * 0,25, и я просто хочу, чтобы 4 сетки 0,25 * 0,25, чтобы иметь то же значение, что и сетка 0,5 * 0,5.

Я посмотрел на растр, но, похоже, ничего не могу с ним поделать.


person AntonyDW    schedule 16.01.2014    source источник


Ответы (3)


Вот способ сделать это с помощью plyr::ddply() - вероятно, это будет немного медленным для вашего размера таблицы, в зависимости от того, как часто вы хотите менять сетку. Я подумаю, как это сделать с помощью data.table, который должен быть быстрее:

require(plyr)
# make your data frame
I<-data.frame(lat=seq(0.5,1000,0.5),lon=1,tmp=sample(1:100,2000,replace=T))

# make an adjustment grid
k<-expand.grid(c(0,0.25),c(0,0.25),0)

# use plyr:ddply() to expand out each entry into the correponding 4 entries
new_I<-ddply(I,.(lat,lon),function(x)as.list(x)+k)
colnames(new_I)<-c("lat","lon","newlat","newlon","tmp")

head(new_I)

  lat lon newlat newlon tmp
1 0.5   1   0.50   1.00  64
2 0.5   1   0.75   1.00  64
3 0.5   1   0.50   1.25  64
4 0.5   1   0.75   1.25  64
5 1.0   1   1.00   1.00  31
6 1.0   1   1.25   1.00  31

На самом деле, если подумать, вот лучший способ с точки зрения времени (хотя это немного похоже на взлом и дает вам меньше контроля для дополнительной обработки данных, которую вы, возможно, захотите выполнить в будущем), но это занимает 6,5 секунды для 2 м >> 8м рядов.

# make your data frame
I<-data.frame(lat=seq(0.5,1000000,0.5),lon=1,tmp=sample(1:100,2000000,replace=T))

# make an adjustment vector
v<-rep(0.25,times=2000000)

# make 3 new tables, apply the vector appropriately, and rbind
I_latshift<-I
I_lonshift<-I
I_bothshift<-I

I_latshift$lat<-I_latshift$lat+v
I_lonshift$lon<-I_lonshift$lon+v
I_bothshift$lat<-I_bothshift$lat+v
I_bothshift$lon<-I_bothshift$lon+v

I<-rbind(I,I_bothshift,I_latshift,I_lonshift)

# sort it for neatness
I<-I[with(I, order(lat, lon)), ]


head(I)

         lat  lon tmp
1       0.50 1.00   3
6000001 0.50 1.25   3
4000001 0.75 1.00   3
2000001 0.75 1.25   3
2       1.00 1.00  88
6000002 1.00 1.25  88
person Troy    schedule 16.01.2014

В пакете R raster есть решение. Это выглядит следующим образом

library("ncdf4")
library("raster")
nc <- nc_open("my_file.nc")
lon <- ncvar_get(nc, "lon")
lat <- ncvar_get(nc, "lat")
time <- ncvar_get(nc, "time")
dname <- "pre"        ## pre for the short name of precpitation 
nlon <- dim(lon)
nlat <- dim(lat)
nt <- dim(time)
lonlat <- expand.grid(lon, lat)    # make grid of given longitude and latitude 
pr.array <- ncvar_get(nc, dname)
dlname <- ncatt_get(nc, dname, "long_name")
dunits <- ncatt_get(nc, dname, "units")
fillvalue <- ncatt_get(nc, dname, "_FillValue")

pr.vec.long <- as.vector(pr.array)
pr.mat <- matrix(pr.vec.long, nrow = nlon * nlat, ncol = nt)
pr.df <- data.frame(cbind(lonlat, pr.mat))

pr_c <- pr.df[ ,-c(1:2)]
 ### Specific region have been clipped out from global datafile by 
## selecting lon and lat range and extract regridded data at 1lon 1lat
 ## resolution.  

x0 <- seq(67.5, 98.5, by = 1) ## choose different resolution, eg. by = 0.5 
y0 <- seq(6.5, 37.5, by = 1)


m <- cbind(x0, y0)
m <- as.data.frame(m)
s <- rasterFromXYZ(m)
pts <- expand.grid(x0, y0)
pos <- pr.df[ ,c(1:2)]
l_pr <- apply(pr_c, 2, function(x) cbind(pos, x))
colnm = c("x","y","z")
for (j in seq_along(l_pr)){
  colnames(l_pr[[j]]) <- colnm
}

pr_rstr <- lapply(l_pr, function(x) rasterFromXYZ(x))
## Use resample command to regrid the data, here nearest neighbor method can also be chosen by setting method = "ngb"
pr_bn <- lapply(pr_rstr, function(x) resample(x, s, method = "bilinear"))
pr_extr <- lapply(pr_bn, function(x) extract(x, pts))
df_pr <- do.call("cbind", lapply(pr_extr, data.frame))
## write dataframe in csv format
write.csv(df_pr, "my_data_regridded_1.csv")

Надеюсь, это послужит цели.

person Pankaj    schedule 16.08.2017
comment
Как сохранить файл в netcdf? - person user2543; 30.07.2018
comment
@ user2543, вы можете посетить ссылку, чтобы преобразовать файл .csv в формат NetCDF. - person Pankaj; 30.07.2018

Это не решение R, но просто чтобы указать, что вы можете использовать CDO для очень простой регистрации файлов netcdf из командной строки в среде Linux / MAC OS. Из вашего описания кажется, что вы хотите использовать интерполяцию ближайшего соседа, которая для регулярной сетки 0,25 градуса будет

cdo remapnn,r1440x720 in.nc out.nc

Однако вы также можете использовать консервативное переназначение первого или второго порядка. Например, для первого заказа:

cdo remapcon,r1440x720 in.nc out.nc

Затем вы можете прочитать в R зарегистрированное поле так же, как вы это делаете сейчас.

person Adrian Tompkins    schedule 17.01.2018
comment
как получить эквивалент разрешения (r) в разных степенях (0,25, 0,5, 1, 1,5, 2, 2,5 градуса) - person Bappa Das; 30.10.2019
comment
легко, просто определите, как точки, что эквивалентно. Например, обычная сетка широты и долготы с разрешением 1 градус имеет 360 точек долготы и 180 точек широты, поэтому вам просто нужно использовать cdo remapcon, r360x180 in.nc out.nc - person Adrian Tompkins; 30.10.2019