Расчет взвешенных пространственных глобальных среднегодовых значений по ячейкам сетки с использованием набора данных netcdf в R - PullRequest
1 голос
/ 19 марта 2019

В настоящее время я работаю над проектом, который включает данные климатической модели, хранящиеся в файле NetCDF.В настоящее время я пытаюсь рассчитать «взвешенные» пространственные годовые «глобальные» средние значения осадков.Я должен сделать это для каждого из 95-летних глобальных данных об осадках, которые у меня есть.Идея состояла бы в том, чтобы как-то применить веса к каждой ячейке сетки, используя косинус ее широты (что означает, что ячейки сетки широты на экваторе будут иметь вес 1 (т. Е. Косинус 0 градусов равен 1), а полюса будут иметьзначение 1 (косинус 90 равен 1)).Затем я смог бы рассчитать средневзвешенные годовые значения на основе усреднения каждой ячейки сетки.

У меня есть идея, как сделать это концептуально, но я не уверен, с чего начать писать скрипт на R, чтобы применить веса ко всем ячейкам сетки, а затем усреднить их для каждого из 95 лет.Я был бы очень признателен за любую помощь в этом, или любые ресурсы, которые могут быть полезны !!!

По крайней мере, я открыл файл .nc и прочитал переменные NetCDF, как показано ниже:

ncfname<-"MaxPrecCCCMACanESM2rcp45.nc"
Prec<-raster(ncfname)
print(Prec)
Model<-nc_open(ncfname)
get<-ncvar_get(Model,"onedaymax")
longitude<-ncvar_get(Model, "lon")
latitude<-ncvar_get(Model, "lat")
Year<-ncvar_get(Model, "Year")

Кроме того, скажем, что я хотел создать временной ряд этих вновь полученных весовых усреднений для определенного местоположения или региона, следующий код, который я ранее использовал, чтобы показать тенденции за 95 лет для одногодневное максимальное количество осадков работает, но мне просто нужно немного изменить его, чтобы использовать среднегодовое взвешенное значение?:

r_brick<-brick(get, xmn=min(latitude), xmx=max(latitude), ymn=min(longitude),                               
ymx=max(longitude), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84   
+no_defs+ towgs84=0,0,0"))
r_brick<-flip(t(r_brick), direction='y')
randompointlon<-13.178
randompointlat<--59.548
Random<-extract(r_brick, 
SpatialPoints(cbind(randompointlon,randompointlat)),method='simple')
df<-data.frame(year=seq(from=1, to=95, by=1), Precipitation=t(Hope))
ggplot(data=df, aes(x=Year, y=Precipitation,   
group=1))+geom_line()+ggtitle("One-day maximum precipitation (mm/day) trend  
for Barbados for CanESM2 RCP4.5")

Кроме того, если это поможет, вот что содержит файл .nc:

 3 variables (excluding dimension variables):
    double onedaymax[lon,lat,time]   (Contiguous storage)  
        units: mm/day
    double fivedaymax[lon,lat,time]   (Contiguous storage)  
        units: mm/day
    short Year[time]   (Contiguous storage)  

 3 dimensions:
    time  Size:95
    lat  Size:64
        units: degree North
    lon  Size:128
        units: degree East

Опять же, любая помощь будет чрезвычайно полезна в этом!Я с нетерпением жду вашего ответа!

Ответы [ 2 ]

1 голос
/ 19 марта 2019

Пожалуйста, задавайте один четкий вопрос за один раз и предоставьте пример данных (через код).

Я не думаю, что вы правильно читаете данные ncdf.Я думаю, что вы должны сделать

library(raster)
ncfname <- "MaxPrecCCCMACanESM2rcp45.nc"
Prec <- brick(ncfname, var="onedaymax")

(делать не использовать nc_open и т. Д.)

Чтобы получить глобальное средневзвешенное значение

Пример данных

library(raster)
r <- abs(init(raster(), 'y'))
s <- stack(r, r, r)

s - это растровый стек со значением 90 на полюсах и 0 на экваторе

Невзвешенное глобальное среднее.Сначала усредните слои, затем ячейки (в этом случае также будет работать обратный порядок)

sm <- mean(s, na.rm=TRUE)
cellStats(sm, mean, na.rm=TRUE)
[1] 45

Теперь используйте взвешивание (чтобы получить меньшее число, высокие высоты получают меньший вес)

# raster with latitude cell values 
w <- init(s, 'y')
# cosine after transforming to radians
w <- cos(w  * (pi/180))
# multiply weights with values
x <- sm * w
# compute weighted average
cellStats(x, sum) / cellStats(w, sum)
#[1] 32.70567

Альтернативное и, возможно, более простое решение заключается в использовании площади каждой ячейки (которая пропорциональна cos(lat)).Результат, возможно, немного более точен (поскольку площадь учитывает не только центр ячейки).

a <- area(s) / 10000
y <- sm * a
cellStats(y, sum) / cellStats(a, sum)
#[1] 32.72697

Позже:

Для временного ряда простоиспользуйте s.

невзвешенный

cellStats(s, mean) 
#layer.1 layer.2 layer.3 
# 45      45      45 

взвешенный

a <- area(s) / 10000
y <- s * a
cellStats(y, sum) / cellStats(a, sum)
# layer.1  layer.2  layer.3 
#32.72697 32.72697 32.72697 
0 голосов
/ 19 марта 2019

Не то, чтобы я хотел отвлечь вас от R, но такого рода вычисления - абсолютный хлеб cdo (операторов климатических данных) прямо из командной строки!

Рассчитать средневзвешенное по пространству (оно учитывает грех широты и может также обрабатывать уменьшенные гауссовы сетки и т. Д.):

cdo fldmean input.nc fldmean.nc 

Рассчитать среднегодовое значение

cdo yearmean input.nc yearmean.nc

Рассчитать временные ряды глобальных среднегодовых значений, комбинируя их (т.е. используйте fldmean.nc в качестве ввода для второй команды), или вы можете сделать это в одной строке, отправив по трубопроводу:

cdo yearmean -fldmean input.nc yearglobal.nc

Что это? Вы хотите рассчитать его для региона латышского бокса, который вы говорите, а не для глобальных средних? Нет проблем, сначала используйте sellonlatbox, чтобы вырезать область

cdo sellonlatbox,lon1,lon2,lat1,lat2 in.nc out.nc 

так что, об этом:

cdo  yearmean -fldmean -sellonlatbox,lon1,lon2,lat1,lat2 in.nc yearregion.nc

Но подождите! Вы хотите конкретное местоположение, а не средний регион? Ну, вы можете выбрать ближайший gridbox к месту с

cdo remapnn,lon=mylon/lat=mylat in.nc out.nc 

, чтобы вы могли получить свою серию годовых усредненных значений с помощью:

cdo yearmean -remapnn,lon=mylon/lat=mylat in.nc yearmylocation.nc

возможностей много ... установите его с

sudo apt install cdo 

и посмотрите документацию здесь: https://code.mpimet.mpg.de/projects/cdo/wiki/Cdo#Documentation

...