построение данных netcdf для выбранного региона

Question

Мне трудно работать с имеющимися у меня данными netcdf. Я просмотрел примеры здесь ( Создание файла netcdf с уровнями в R ), но я думаю, что я что-то упустил.

Я пытаюсь построить глубину смешанного слоя для региона к югу от 35 градусов. Данные можно найти здесь (последний файл внизу страницы): http://www.ifremer.fr/cerweb/deboyer/mld/Surface_Mixed_Layer_Depth.php

файл содержит 7 переменных, каждая из которых содержит широту, долготу, время (12 месяцев) и значение глубины смешанного слоя.

Пока у меня есть:

    MLD <- "mld_DReqDTm02_c1m_reg2.0.nc"
    MLD <- nc_open(MLD)
    print(MLD)

есть 7 переменных, и я хочу только 'mld'

lon <- ncvar_get(MLD, varid = "lon")
lat <- ncvar_get(MLD, varid = "lat")

summary(lon)
summary(lat)



MLD$dim$time$units

MLD_1.array <- ncvar_get(MLD, "mld")
dim(MLD_1.array)

length(lon)
length(lat)

ndvi.slice <- MLD_1.array[, , 12] 
dim(ndvi.slice)

mld.vec.long <- as.vector(MLD_1.array)
length(mld.vec.long)

nlon <- dim(lon)
nlat <- dim(lat)
lonlat <- expand.grid(lon, lat)

t <- ncvar_get(MLD, "time")
tunits <- ncatt_get(MLD, "time", "units")
nt <- dim(t)

dname <- "mld" 

tmp.mat <- matrix(mld.vec.long, nrow = nlon * nlat,    ncol = nt)
dim(tmp.mat)

head(na.omit(tmp.mat))

lonlat <- expand.grid(lon, lat)
tmp.df02 <- data.frame(cbind(lonlat, tmp.mat))
 names(tmp.df02) <- c("lon", "lat", "Jan", "Feb", "Mar", "Apr", "May", 
"Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

library(reshape)
tmp.df03 <- melt(tmp.df02, id=c("lat","lon"))

tmp.df04 <- subset(tmp.df03, lat >= "-35" & lat <="-80")

 tmp.df04[tmp.df04 ==1.000e+09] <- NA
 tmp.df04
 summary(tmp.df04)


 target = c(-180, 180, -90, -20)
 w <- crop(rgeos::gBuffer(spTransform(countriesLow, CRS(pprj)), width = 0))

Я могу построить это с помощью ggplot после извлечения данных в файл .csv, удаления значений маски, повторного импорта и повторного построения графика (это занимает некоторое время ...).

Есть ли способ, которым я могу извлечь одну переменную (mld), а затем найти максимальное значение для каждого месяца (времени) и построить график? Я знаю, что мой код очень грязный и нестабильный ...

любая помощь будет высоко ценится! благодарю вас!

Answer 1

Причина сбоя brick(f) состоит в том, что есть одна широта, которая не является регулярной (последняя, ​​которая должна была быть 90 или отсутствовала)

library(ncdf4)
library(raster)
MLD <- nc_open("c:/temp/mld_DReqDTm02_c1m_reg2.0.nc")
lat <- ncvar_get(MLD, varid = "lat")
tail(lat)
#[1] 80.0 82.0 84.0 86.0 88.0 89.5

Поскольку это только одна строка (самая северная), возможно, можно игнорировать это и делать, как показывает @rar. Здесь с немного другим подходом

mld <- ncvar_get(MLD, "mld")
a <- aperm(mld, c(2,1,3))
b <- flip(brick(a), 'y')
extent(b) <- c(0,360,-90,90)
b <- rotate(b)
b <- reclassify(b, cbind( 1e+09, NA))
names(b) = month.abb
plot(b,1)

Теперь вы можете использовать кроп

bb <- crop(b, extent(-180,180,-90,35))
plot(bb,1)

Или извлечь по координатам

extract(b, cbind(0,0))
#          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec
#[1,] 27.80343 23.17406 25.27729 19.04638 16.96231 18.96719 18.79073 21.04775 30.04144 38.48611 27.82375 20.18804
> 

Answer 2

Надеюсь, это поможет!

library(ncdf4)
library(raster)

#Reading netcdf file and extracting lat, lon and variable
MLD <- nc_open("D:/Personal/test/mld_DReqDTm02_c1m_reg2.0.nc")
lon <- ncvar_get(MLD, varid = "lon")
lat <- ncvar_get(MLD, varid = "lat")
mld <- ncvar_get(MLD, "mld")

# using 'raster' package to read the temporal values of variable into a raster
#r1<-flip(raster(t(matrix(mld[,,1], nrow = 180,    ncol = 90))),direction="y")
e<-extent(min(lon),max(lon),min(lat),max(lat))
#extent(r1)<-e
R<-stack()
# creating raster stack of the time series data
for(i in 1:12){
  r1<-flip(raster(t(matrix(mld[,,i], nrow = 180,    ncol = 90))),direction="y")
  extent(r1)<-e
  R<-stack(R,r1)
} # IMPORTANT: This raster stack could be additionally cropped to extract the user's area of interest.
plot(R)

# Extracting max/min values for each time (raster layer) into a dataframe.
df<-data.frame(Months=month.abb)
df$Months <- factor(df$Months, levels = df$Months)
df$months<-c(1:12)
df$MLD_max<-maxValue(R)
df$MLD_min<-minValue(R)
# using 'ggplot2' package
ggplot(df)+geom_point(aes(Months,MLD_min))+geom_line(aes(months,MLD_min))