Извлечение всех уровней из файла NetCDF в R - PullRequest
Новая
       36

Извлечение всех уровней из файла NetCDF в R

0 голосов
/ 22 мая 2019

Я пытаюсь извлечь все уровни из определенного файла NetCDF в R. Я могу сделать это вручную, извлекая каждый уровень в виде одной строки кода, а затем объединяя их в виде фрейма данных. Но это очень долго, когда у меня много файлов. Можно ли извлечь все 43 слоя в один файл?

Я использовал это Как извлечь все уровни из файла netcdf с помощью растрового пакета? и Создание файла netcdf с уровнями в R в качестве руководства

По существу данные о нитратах
https://www.nodc.noaa.gov/cgi-bin/OC5/woa18/woa18oxnu.pl имеет концентрацию на 43 различных глубинах. Можно ли извлечь все глубины для определенного места?

Я могу сделать это для одного уровня. Но каждый уровень представляет глубину. Можно ли пройти все уровни?

Я также не понимаю 3-го предупреждающего сообщения: В .getCRSfromGridMap4 (atts): не могу обработать эти части CRS: epsg_code = EPSG: 4326

Я получаю другой результат (0,5 за январь 1-го уровня), но мой коллега получает 1,4 за январь 1-го уровня). Моя ошибка связана с приведенным выше предупреждением? 

#this works
Nit_Jan <- brick("~woa18_all_n01_01.nc", stopIfNotEqualSpaced =    
FALSE, varname = "n_an", level = 1)

#this doesn't
Nit_Jan <- brick("~woa18_all_n01_01.nc", stopIfNotEqualSpaced = 
FALSE, varname = "n_an", level = 1:43)

Warning messages:
1: In if (level <= 0) { :
the condition has length > 1 and only the first element will be used 
2: In if (oldlevel != level) { :
the condition has length > 1 and only the first element will be used
3: In .getCRSfromGridMap4(atts) : cannot process these parts of the   
CRS:epsg_code=EPSG:4326

Я бы хотел построить нитрат по глубине

1 Ответ

0 голосов
/ 22 мая 2019

Здесь есть некоторая путаница, потому что файл имеет «уровни» (четвертое измерение), но количество уровней равно одному (поэтому четвертого измерения нет).Код, вероятно, должен это обнаружить, но сейчас вам нужно добавить lvar=4, чтобы получить нужный объект. 

library(raster)
f <- "woa18_all_n01_01.nc"
b <- brick(f, var="n_oa", lvar=4)
b
#class      : RasterBrick 
#dimensions : 180, 360, 64800, 43  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : woa18_all_n01_01.nc 
#names      : X0, X5, X10, X15, X20, X25, X30, X35, X40, X45, X50, X55, X60, X65, X70, ... 
#meters     : 0, 800 (min, max)
#varname    : n_oa 
#level      : 1

И теперь вы можете сделать 

pt <- cbind(121.5, 27.5)
e <- extract(b, pt)
e[1:5]
#[1] 10.43725 10.37617 10.23662 13.76292 13.65862

Warning # 3

3: в .getCRSfromGridMap4 (atts): невозможно обработать эти части CRS: epsg_code = EPSG: 4326

можно игнорировать;но я исправлю это в следующей версии.Я думаю, что лучшая вещь здесь - 

crs(b) <- "+init=EPSG:4326"

PS: версия для разработки растра теперь ведет себя лучше: 

f <- "woa18_all_n01_01.nc"
brick(f, var="n_oa")
#class      : RasterBrick 
#dimensions : 180, 360, 64800, 43  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +init=EPSG:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#source     : woa18_all_n01_01.nc 
#names      : X0, X5, X10, X15, X20, X25, X30, X35, X40, X45, X50, X55, X60, X65, X70, ... 
#depth (meters): 0, 800 (min, max)
#varname    : n_oa 
#level      : 1
...