Для цикла открытия и доступа к данным спутникового netcdf
Я использую спутниковые данные, чтобы посмотреть на уровень освещенности на поверхности океана. (любые данные из здесь работают, так как я просто пытаюсь понять, как открывать и работать с данными)
В приведенном ниже коде я сохранил все файлы .nc в папке на рабочем столе, затем импортировал их все в files и запустил цикл открытия файла.
Есть ли способ 1) обрезать данные extent(-180, -176, -40, -36) , а затем 2) извлечь значения широты, долготы и номинала? Я предполагаю, что доступно несколько глубин, но я не знаю, как это проверить.
library(raster)
setwd("/Users/AquaModis_PAR")
nc.data<-list()
files<-list.files(pattern = ".nc",full.names=TRUE)
df=NULL
for(i in seq_along(files)) {
par[i] = brick(files[i], stopIfNotEqualSpaced = FALSE, varname = "par")
}
Ответы 2
Чтобы загрузить и обрезать данные, сделайте что-то вроде этого:
library(raster)
files <- list.files(pattern = ".nc", full.names=TRUE)
region <- extent(-180, -176, -40, -36)
data <- lapply(files, function(path) {
# Load and crop.
d <- brick(path, stopIfNotEqualSpaced = FALSE, varname = "chlor_a") |>
crop(region)
# Get raster values.
values <- getValues(d)
# Get coordinates.
coords <- coordinates(d)
cbind(coords, values) |>
data.frame() |>
setNames(c("lon", "lat", "value")) |>
na.omit()
})
Посмотрите на данные из первого файла:
> head(data[[1]])
lon lat value
28 -178.8542 -36.02084 0.06653178
29 -178.8125 -36.02084 0.06850998
30 -178.7708 -36.02084 0.07330192
31 -178.7292 -36.02084 0.07878597
34 -178.6042 -36.02084 0.07993947
35 -178.5625 -36.02084 0.08084643
Вероятно, существуют более элегантные способы извлечения данных в фрейм данных, но этот справляется со своей задачей.
na.omit() в конце предназначен для удаления записей, в которых значение отсутствует (в данных, которые я просматривал, их было довольно много).
Для полноты вот файлы, на которых я тестировал:
AQUA_MODIS.20030101.L3m.DAY.CHL.chlor_a.4km.nc
AQUA_MODIS.20030102.L3m.DAY.CHL.chlor_a.4km.nc
AQUA_MODIS.20030103.L3m.DAY.CHL.chlor_a.4km.nc
Это очень неэффективно, поскольку вы выбрасываете 99,9% только что прочитанных данных, немедленно обрезая их. Другого варианта, когда вы используете raster, я полагаю... нет.
Продукты MODIS AQUA L3m имеют глобальное покрытие. Вас интересует небольшая область мира, (4 / 360) * (4 / 180) = 0.0247% и если вы используете raster::brick(), вы считываете 100% данных, а затем отбрасываете 99,9753%. Не очень эффективно.
С пакетом ncdfCF вы сможете сделать это гораздо эффективнее. Это также поможет вам ответить на другие ваши вопросы. Но сначала некоторые пояснения.
Данные MODIS распространяются в файлах NetCDF. Эти файлы имеют атрибуты, которые вам следует проверить для интерпретации данных. К сожалению, пакеты швейцарского армейского ножа, такие как raster (а также его преемники terra и stars), не обеспечивают (простого) доступа к этим атрибутам, поэтому вам понадобится что-то более специально адаптированное для файлов NetCDF. RNetCDF и ncdf4 предоставляют доступ к атрибутам, но не очень удобным для пользователя способом. ncdfCF, с другой стороны, дает вам легкий доступ:
> ds <- open_ncdf("AQUA_MODIS.20100101.L3m.DAY.PAR.par.4km.nc")
> ds[["par"]]
Variable: [0] par | Photosynthetically Available Radiation, R. Frouin
Dimensions:
id axis name long_name length values
1 X lon Longitude 8640 [-179.97917 ... 179.97917]
0 Y lat Latitude 4320 [89.97916 ... -89.97917]
Attributes:
id name type length value
0 long_name NC_CHAR 49 Photosynthetically Available Radiation, R. Frouin
1 scale_factor NC_FLOAT 1 0.0020000000949949
2 add_offset NC_FLOAT 1 65.5
3 units NC_CHAR 20 einstein m^-2 day^-1
4 standard_name NC_CHAR 53 surface_downwelling_photosynthetic_photon_flux_...
5 _FillValue NC_SHORT 1 -32767
6 valid_min NC_SHORT 1 -32750
7 valid_max NC_SHORT 1 32250
8 reference NC_CHAR 208 Frouin, R., Ligner, D.W., and Gautier, C., 1989...
9 display_scale NC_CHAR 6 linear
10 display_min NC_FLOAT 1 0
11 display_max NC_FLOAT 1 76.1999969482422
Это ответ на ваш вопрос о нескольких глубинах. Переменная «par» имеет размеры lon и lat, поэтому имеется только 1 уровень. (Это имеет смысл: PAR — это падающее излучение (см. атрибут «standard_name»), на уровне моря для этого океанического продукта ослабление в результате фотосинтеза или поглощение в воде не учитывается.)
Поскольку вы хотите перебрать несколько файлов, вам следует различать даты отдельных файлов. Предполагая, что вы используете «дневные» данные, вы можете получить даты, используя sapply(strsplit(files, '.', fixed = T), function(d) d[2]). Обратите внимание, что «дата» — это гибкая конструкция в данных MODIS. Атрибуты набора данных покажут вам, что он охватывает период, превышающий 24 часа (атрибуты «time_coverage_start» и «time_coverage_end»). Взять дату из имени файла, вероятно, будет разумным выбором.
Объединив это в функцию, вы получите:
read_par_array <- function(fn, extent) {
out <- lapply(fn, function(f) {
ds <- open_ncdf(f)
subset(ds[["par"]], extent)
})
out <- abind(out, along = 3)
dimnames(out)[[3]] <- sapply(strsplit(fn, '.', fixed = T), function(d) d[2])
out
}
read_par_df <- function(fn, extent) {
out <- lapply(fn, function(f) {
ds <- open_ncdf(f)
data <- subset(ds[["par"]], extent)
as.data.frame(data)
})
names(out) <- sapply(strsplit(fn, '.', fixed = T), function(d) d[2])
out
}
А затем вы вызываете его со своими файлами и степенью интереса:
> library(ncdfCF)
> library(abind)
> files <- list.files(pattern = ".nc$", full.names = TRUE)
> extent <- list(X = c(-180, -176), Y = c(-40, -36))
> par <- read_par_array(files, extent)
read_par_array() создает трехмерный массив с размерами долготы, широты и даты. Для этого вам понадобится пакет abind.
> par[1:3, 1:3, ]
, , 20100101
-36.02084 -36.0625 -36.10417
-179.97917 NA NA NA
-179.9375 NA NA NA
-179.89583 NA NA NA
, , 20100102
-36.02084 -36.0625 -36.10417
-179.97917 66.924 67.634 68.238
-179.9375 67.118 68.226 68.530
-179.89583 67.770 68.330 68.716
, , 20100103
-36.02084 -36.0625 -36.10417
-179.97917 39.712 45.216 47.190
-179.9375 38.406 39.846 42.994
-179.89583 37.066 33.886 39.686
, , 20100104
-36.02084 -36.0625 -36.10417
-179.97917 NA NA NA
-179.9375 NA NA NA
-179.89583 NA NA NA
, , 20100105
-36.02084 -36.0625 -36.10417
-179.97917 68.226 68.160 68.162
-179.9375 68.254 68.184 68.230
-179.89583 68.302 68.312 68.318
read_par_df() дает список с 1 элементом на файл. Каждый элемент назван в честь даты в имени файла MODIS. В каждом списке есть 1 data.frame с именами широты в столбцах и долготы в строках:
> names(par)
[1] "20100101" "20100102" "20100103" "20100104" "20100105"
> par[["20100102"]][1:3, 1:3]
-36.02084 -36.0625 -36.10417
-179.97917 66.924 67.634 68.238
-179.9375 67.118 68.226 68.530
-179.89583 67.770 68.330 68.716
В данных довольно много NA, и я предполагаю, что они вызваны облачностью в течение дня: алгоритм предназначен только для условий ясного неба. В пасмурную погоду количество ФАР на поверхности моря значительно снижается, поэтому, если вы можете смириться с этим в своих анализах, пожалуйста, продолжайте. В противном случае обратите внимание на продукты MCD18C2 (5-минутные гранулы, поэтому требуется немного больше усилий...)
Отказ от ответственности: я являюсь разработчиком пакета ncdfCF. У меня также есть докторская степень по вычислению PAR с использованием данных MODIS, полученная в 2004 году. Я начал свои исследования за несколько недель до того, как MODIS Terra появилась на свет.
спасибо, это сработало! У меня 2 вопроса: 1) для
head(data[[1]])означает ли это, что он принимает только 1-й слой (т. е. поверхностное значение)? как мне узнать, сколько там глубин/слоев? и 2) есть ли способ добавить время? Я не уверен, что это есть в файлах, потому что, когда я использовал nc_open, я получал 4 измерения: широту, долготу, ребро и восьмибитный цвет. Если вы заметили, файлы сохраняются по дате. Это что-то, что мне нужно добавить вручную?