создать растровый кирпич из ncdf нерегулярно расположенных x и y - PullRequest
Новая
       75

создать растровый кирпич из ncdf нерегулярно расположенных x и y

0 голосов
/ 09 января 2019

Я довольно новичок в использовании R для целей ГИС. У меня есть файл netcdf, содержащий несколько переменных с несколькими измерениями (x, y, z, значение и время). Я пытаюсь превратить это в растровый кирпич. Данные довольно большие, поэтому мне нужно извлечь данные из указанного временного окна и z (глубины). Это не было проблемой, и извлеките массив с соответствующими размерами, используя приведенный ниже код. 

library(ncdf4)
library(raster)
t <- ncvar_get(nc, "model_time")
t1<-ncvar_get(nc,"model_time_step")
tIdx<-t[t> 20120512 & t < 20120728]
tIdx2<-which(t> 20120512 & t < 20120728)
# Depths profiles < 6 meters
dIdx<-which(nc$dim$depthu$vals <6)
# ncdf dimension lengths
T3      <- nc$var[[7]]
varsize <- T3$varsize
# Define the data (depths,time,etc.) you wish to extract from the ncdf
start <- c(x = 1, y= 1,depthu=1, time_counter = min(tIdx2))
count <- c(x = max(varsize[1]), y = max(varsize[2]),depthu=1, time_counter = 
max(tIdx2)-min(tIdx2)+1)
# order of the dimensions
dim.order <- sapply(nc$var$votemper$dim, function(x) x$name)
temp<-ncvar_get(nc,"votemper",start=start[dim.order],count=count[dim.order])
nc$var$votemper

Пример моих данных (отбрасывая глубину / z и временные измерения) 

   temp<-structure(c(0,0,0,0,0,0,0,15.7088003158569,15.3642873764038,14.9720048904419,,15.9209365844727,14.9940872192383,15.0184164047241,15.0260219573975, 0,15.7754755020142, 15.424690246582, 15.6697931289673,15.6437339782715, 0,15.6151847839355, 15.5979156494141, 15.6487197875977,15.432520866394), .Dim = c(x = 5L, y = 5L))

Широты и долготы, извлеченные из ncdf, расположены нерегулярно и имеют два измерения каждое (т.е. нерегулярно разнесенные значения широты и долготы для каждой ячейки) 

lon<-structure(c(-71.2870483398438,-71.2038040161133,-71.1205596923828,-71.0373153686523, -70.9540710449219, -71.2887954711914, -71.2055587768555,-71.122314453125, -71.0390701293945,-70.9558258056641,-71.2905654907227,-71.2073211669922,-71.1240844726562,-71.0408401489258,-70.9576034545898,-71.292350769043,-71.209114074707, -71.1258773803711, -71.0426330566406,-70.9593963623047, -71.2941513061523, -71.2109222412109, -71.127685546875,-71.0444488525391, -70.9612045288086), .Dim = c(5L, 5L))

lat<-structure(c(38.5276718139648, 38.529125213623, 38.5305824279785,38.532039642334, 38.5334968566895, 38.5886116027832, 38.5900802612305,38.591552734375, 38.5930252075195, 38.5944976806641, 38.6494789123535,38.6509628295898, 38.6524467468262, 38.6539344787598, 38.6554222106934,38.7102699279785, 38.7117652893066, 38.713264465332, 38.7147674560547,38.7162704467773, 38.7709808349609, 38.7724952697754, 38.7740097045898,38.7755241394043, 38.777042388916), .Dim = c(5L, 5L))

Обычно я генерирую растровый кирпич из этих данных, используя 

Temp_brick <- brick(temp, xmn=min(lat), xmx=max(lat), ymn=min(lon), ymx=max(lon),transpose=T)
Temp_brick<-t(flip(Temp_brick,1))

Это, однако, не учитывает нерегулярный интервал, и значения растровых ячеек расположены в неправильном положении (lon, lat). Я искал переполнение стека и другие источники справки ГИС, и я не могу найти аналогичную проблему с решением, или я не задаю правильный вопрос. Я не совсем уверен, как это сделать. Не уверен, следует ли это учитывать при извлечении данных из netcdf или же это следует делать после создания растрового кирпича без определенного экстента. Я попытался найти способ определить lon lats для растра без какой-либо удачи. Попытался преобразовать lon, lat и value в три столбца данных, а затем использовать функцию raster :: rasterFromXYZ. Это не будет работать достаточно быстро для размера данных, с которыми я имею дело, что в действительности составляет 197 (x) * 234 (y) * 2 (z) * 900 (время) * 5 (переменные) * 12 ( лет (отдельные файлы netcdf).

Любая помощь очень ценится

1 Ответ

0 голосов
/ 10 января 2019

опция с акимой, чтобы сначала преобразовать данные в регулярную сетку, а затем превратить их в растр: 

    # define the regular lon lat or just pass the nx, ny param to interp functions
    lonlat_reg <- expand.grid(lon = seq(min(lon), max(lon), length.out = 5), 
        lat = seq(min(lat), max(lat), length.out = 5))

    # interp irregular data to a regular grid
    # both solution return the same results because 
    # i've define the regular grid as akima default
    test <- interp(x = as.vector(lon), y = as.vector(lat), z = as.vector(temp), 
        xo = unique(lonlat_reg[,"lon"]), yo = unique(lonlat_reg[,"lat"]), 
        duplicate = "error", linear = FALSE, extrap = FALSE)

    test <- interp(x = as.vector(lon), y = as.vector(lat), z = as.vector(temp), 
        nx = 5, ny = 5, linear = FALSE, extrap = FALSE)

    # turn into a raster
    test_ras <- raster(test)

Проверьте аргументы функции , чтобы выбрать выполненную интерполяцию и т. Д., И будьте осторожны, если вы используете экстраполяцию!

Я также видел, что метод

Приветствия

...