Как извлечь подмножество из файла netCDF, используя границы широты / долготы в R
У меня есть файл netCDF, из которого я хочу извлечь подмножество, определяемое границами широты / долготы (т. Е. Поле широты / долготы), используя пакет 'ncdf' в R.
Резюме моего файла netCDF ниже. Он имеет два измерения (широта и долгота) и 1 переменная (10U_GDS4_SFC). По сути, это сетка широта / долгота, содержащая значения ветра:
[1] "file example.nc has 2 dimensions:"
[1] "lat_0 Size: 1280"
[1] "lon_1 Size: 2560"
[1] "------------------------"
[1] "file example.nc has 1 variables:"
[1] "float 10U_GDS4_SFC[lon_1,lat_0] Longname:10 metre U wind component Missval:1e+30"
Переменная широты варьируется от +90 до -90, а переменная долготы от 0 до 360.
Я хочу извлечь подмножество общей сетки, используя следующие географические угловые границы:
нижний левый угол: широта: 34,5˚, длина: 355˚, левый верхний угол: широта: 44,5˚, длина: 355˚, правый верхний угол: широта: 44,5˚, длина: 12˚, правый нижний угол: широта: 34,5˚, Длинный: 12˚
Я знаю, что части переменной могут быть извлечены с помощью get.var.ncdf()
команда (пример ниже):
z1 = get.var.ncdf(example.nc, "10U_GDS4_SFC", start=c(11,26), count=c(5,5))
Тем не менее, я не могу понять, как можно включить широту / долготу, чтобы в итоге я получил заданную пространственную сетку, содержащую значения переменных. Я новичок в работе со значениями netCDF в R, и любые советы будут с благодарностью. Большое спасибо!
2 ответа
В принципе вы 2/3 пути туда. Конечно, вы можете создать начальные индексы, используя что-то вроде этого:
require(ncdf4)
ncFile <- nc_open( MyNetCDF )
LonStartIdx <- which( ncFile$dim$lon$vals == 355)
LatStartIdx <- which( ncFile$dim$lat$vals == 34.5)
Сделайте то же самое для подсчета. Затем прочитайте переменную, которую вы хотите
MyVariable <- ncvar_get( ncFile, varName, start=c( LonStartIdx, LatStartIdx), count=...)
Однако в твоем случае тебе не повезло, насколько я знаю. Процедуры чтения / записи netcdf делают свое дело последовательно. Ваша сетка оборачивается, так как у вас есть координаты, которые идут от 0 до 360 по долготе, и вы заинтересованы в поле, которое содержит нулевой меридиан.
Для вас (при условии, что у вас не слишком много данных) было бы более целесообразно прочитать всю таблицу в R, а затем использовать либо subset
или создать индексы, используя which
и вырезать свою "коробку" в R.
ncFile <- nc_open( MyNetCDF )
LonIdx <- which( ncFile$dim$lon$vals > 355 | ncFile$dim$lon$vals < 10)
LatIdx <- which( ncFile$dim$lat$vals > 34.5 & ncFile$dim$lat$vals < 44.5)
MyVariable <- ncvar_get( ncFile, varName)[ LonIdx, LatIdx]
nc_close(ncFile)
Замечание: я предпочитаю ncdf4
Я нахожу синтаксис немного легче запомнить (и было еще одно преимущество перед старым R-пакетом netcdf, которое я забыл...)
Хорошо. Комментарии не могут быть до тех пор, пока они мне понадобятся, поэтому я обновил ответ. Не беспокойтесь. Давайте рассмотрим вопросы шаг за шагом.
which
Функциональный способ будет работать. Я использую это сам.Данные будут в том же формате, что и в файле netcf, но я не слишком уверен, есть ли какая-то проблема с меридианом 0 (думаю, да). Возможно, вам придется поменять местами две половинки, выполнив что-то вроде этого (замените соответствующую строку во втором примере)
LonIdx <- c(which( ncFile$dim$lon$vals > 355) , which( ncFile$dim$lon$vals < 10) )
Это меняет порядок координатных индексов, так что сначала идет западная часть, а затем восточная.
Переформатирование всего в кадр данных 2x3 возможно. Возьмите данные, которые возвращает мой 2-й пример кода (будет матрица, [lon x lat]. Также получите значения координат из
lon <- ncFile$dim$lon$val[LonIdx]
(или как называется долгота в вашем примере, то же самое для
lat
). Затем собрать матрицу, используяcbind( rep(lat, each=length(lon)), rep(lon,length(lat)), c(myVariable) )
Координаты, конечно, будут такими же, как в файле netcdf...
Вам необходимо проверить последний cbind, так как я уверен только на 98%, что я не перепутал координаты. В сценариях R, которые я нашел на своем рабочем столе, я использую циклы, которые являются... злыми... Это должно быть (немного?) Быстрее, а также более разумно.
Вы также можете использовать CDO, чтобы сначала извлечь область из командной строки bash, а затем прочитать файл в R:
cdo sellonlatbox,-5,12,34.5,44.5 in.nc out.nc
Я отмечаю в вышеупомянутой дискуссии, что была проблема, касающаяся порядка широт. Вы также можете использовать команду CDO "invertlat", чтобы разобраться в этом за вас.
Если вы используете Linux, этого легко добиться с помощью nctoolkit ( https://nctoolkit.readthedocs.io/en/latest/):
import nctoolkit as nc
data = nc.open_data("example.nc")
data.clip(lon = [-12, -5], lat = [35.4, 44.5])