Блог

Я имею дело с файлом netcdf4, загруженным из службы мониторинга морской среды Copernicus, и мне кажется, я что-то здесь упускаю…

Рассмотрите мой подход в следующем коде для извлечения значений из файла nc. Моя задача — повторно проанализировать данные, хранящиеся во временных рядах, с помощью линейной регрессии, чтобы получить среднее изменение SST в день на пиксель (не совсем уверен, получаю ли я статистику, которую я хочу здесь! Я открыт для предложений! в любом случае, это не тема моей темы).

 library(ncdf4)
library(raster)
nc_data <- nc_open('L4_analyzed_sst_005res_25081981_31122018.nc')
print(nc_data) #it shows metadata

lon <- ncvar_get(nc_data, "lon")
lat <- ncvar_get(nc_data, "lat", verbose = F)
t <- ncvar_get(nc_data, "time")
t1 <- as.POSIXct(t, origin="1981-01-01", format="%Y-%m-%d", tz = "UTC")

sst.array <- ncvar_get(nc_data, "analysed_sst") # store the data in a 3-dimensional array (long, lat, time)
dim(sst.array)

fillvalue <- ncatt_get(nc_data, "analysed_sst", "_FillValue")
fillvalue
nc_close(nc_data)

sst.array[sst.array == fillvalue$value] <- NA #fill the no value with NA

str(sst.array)

##################################################################################################
#1. sst change for the entire period
#######################################################################################
slope_matrix <- matrix(nrow = length(lon), ncol = length(lat))
for (i in 1:dim(sst.array)[1]){
  for (j in 1:dim(sst.array)[2]){
    value <- sst.array[i,j,]
    if (anyNA(value)==FALSE){ #this "if" speed up the process because in my case the only NA are noValue
    extracted_data <- data.frame(value=value, t=t1)
    extracted_data$quarter <- "null"
    extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="01"|substr(extracted_data$t, 6,7)=="02"|substr(extracted_data$t, 6,7)=="03", "Q1", extracted_data$quarter)
    extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="04"|substr(extracted_data$t, 6,7)=="05"|substr(extracted_data$t, 6,7)=="06", "Q2", extracted_data$quarter)
    extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="07"|substr(extracted_data$t, 6,7)=="08"|substr(extracted_data$t, 6,7)=="09", "Q3", extracted_data$quarter)
    extracted_data$quarter <- ifelse(substr(extracted_data$t, 6,7)=="10"|substr(extracted_data$t, 6,7)=="11"|substr(extracted_data$t, 6,7)=="12", "Q4", extracted_data$quarter)
    model <- lm(value ~ time(t)   quarter, data = extracted_data)
    slope_matrix[i,j] <- mean(predict(model) - extracted_data$value)
    }
  }
}

slope_raster <- raster(t(slope_matrix), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS(" proj=longlat  ellps=WGS84  datum=WGS84  no_defs  towgs84=0,0,0"))
plot(slope_raster)
slope_raster <- flip(slope_raster, direction='y')
writeRaster(slope_raster, "sst_change_1981-2018.tif", "GTiff", overwrite=TRUE)
  

Я использовал файл nc в качестве массива и построил вложенный цикл for для моделирования временных рядов в каждом пикселе. Затем я преобразовал матрицу значений в растр, перевернул его, чтобы переориентировать, и записал его как GeoTIFF для дальнейшего использования.
Теперь, когда я открываю его в ArcGIS, он не полностью перекрывается с растром, созданным с помощью инструмента «создать растровый слой NetCDF». Растр, созданный в R, имеет немного меньшее разрешение (0,049 против 0,050).

Вы знаете, в чем здесь проблема?

Спасибо за помощь!

РЕДАКТИРОВАТЬ: здесь вы можете загрузить данные.

Существуют гораздо более простые способы обработки файлов ncdf с пространственными данными. Ниже я открываю файл с raster пакетом (который имеет код на основе R для этого) и с terra пакетом (который использует библиотеку GDAL). Они дают одинаковое разрешение; поэтому можно с уверенностью предположить, что это правильно.

 f <- "L4_analyzed_sst_005res_25081981_31122018.nc"
library(raster)
b <- brick(f)
b
#class      : RasterBrick 
#dimensions : 34, 48, 1632, 13643  (nrow, ncol, ncell, nlayers)
#resolution : 0.05004599, 0.05015772  (x, y)
#extent     : 13.22879, 15.631, 39.52957, 41.23494  (xmin, xmax, ymin, ymax)
#crs        :  proj=longlat  datum=WGS84  no_defs 
#source     : L4_analyzed_sst_005res_25081981_31122018.nc 
#names      : X1981.08.25, X1981.08.26, X1981.08.27, X1981.08.28, X1981.08.29, X1981.08.30, X1981.08.31, X1981.09.01, X1981.09.02, X1981.09.03, X1981.09.04, X1981.09.05, X1981.09.06, X1981.09.07, X1981.09.08, ... 
#Date/time  : 1981-08-25, 2018-12-31 (min, max)
#varname    : analysed_sst 
  

Или с terra

 x <- rast(f)
x
#class       : SpatRaster 
#dimensions  : 34, 48, 13643  (nrow, ncol, nlyr)
#resolution  : 0.05004599, 0.05015772  (x, y)
#extent      : 13.22879, 15.631, 39.52957, 41.23494  (xmin, xmax, ymin, ymax)
#coord. ref. :  proj=longlat  datum=WGS84  no_defs 
#data source : L4_analyzed_sst_005res_25081981_31122018.nc 
#names       : L4__1, L4__2, L4__3, L4__4, L4__5, L4__6, ... 
 
  

Используемое вами разрешение неверно, потому что вы:

 raster(nrow=34, ncol=48, xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat))
#class      : RasterLayer 
#dimensions : 34, 48, 1632  (nrow, ncol, ncell)
#resolution : 0.04900336, 0.04868249  (x, y)
#extent     : 13.25381, 15.60598, 39.55465, 41.20986  (xmin, xmax, ymin, ymax)
#crs        :  proj=longlat  datum=WGS84  no_defs 
  

Вы используете значения lon / lat для центров угловых ячеек. Но что вам нужно указать, так это координаты внешнего края угловых ячеек.

Не имеет отношения к вашему вопросу, но теперь вы можете продолжить и вычислить кварталы (вне любого цикла — нет необходимости делать это снова и снова для каждой ячейки?)

 z <- getZ(b)
month <- as.numeric(substr(z, 6 ,7 ))
quarter <- ((month-1) %% 4 )   1
# or quarter <- trunc((month-1) / 3)   1
table(quarter)
#quarter
#   1    2    3    4 
#3434 3333 3435 3441 
  

И смотрите ?calc Примеры выполнения регрессии на основе ячеек