NetCDF в R: проблемы с созданием карты

Поэтому я использовал следующий учебник в R, чтобы привыкнуть к NetCDF (я использую пакет netcdf4): http://geog.uoregon.edu/bartlein/courses/geog607/Rmd/netCDF_01.htm

Однако я не могу создать хорошее изображение для какого-либо конкретного года в своем наборе данных, который можно найти здесь: https://www.ncdc.noaa.gov/paleo-search/study/19419

Вот мой код, в котором я пытаюсь получить карту на 100 год нашей эры:

library(lattice)
library(ncdf4)
library(chron)
library(RColorBrewer)
setwd('/Users/Nikki/Dropbox/Europe/Drought Maps')
ncname <- "owda"
ncfname <- paste(ncname,".nc",sep="")
dname <- "pdsi"
ncin <- nc_open(ncfname)
print(ncin)

lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
head(lat)

print(c(nlon, nlat))

t <- ncvar_get(ncin, "time")
nt <- dim(t)
head(t)

drought.array <- ncvar_get(ncin, dname)
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
fillvalue <- ncatt_get(ncin, dname, "_FillValue")
dim(drought.array)

creation_date <- ncatt_get(ncin, 0, "creation_date")
Description <- ncatt_get(ncin, 0, "Description")

nc_close(ncin)


m <- 100
drought.slice <- drought.array[m,,]
image(lon, lat, drought.slice, col = rev(brewer.pal(10, "RdBu")))

В частности, я получаю следующее сообщение об ошибке, когда пытаюсь запустить последнюю строку:

Ошибка в image.default(lon, lat, drought.slice, col = rev(brewer.pal(10, : размеры z) не равны длине(x)(-1) умноженной на длину(y)(-1)

Может кто-нибудь объяснить, что я делаю неправильно?

Кстати, чтобы дать вам некоторое представление о структуре моих данных, вот вам:

>print(ncin)
File owda.nc (NC_FORMAT_NETCDF4_CLASSIC):
 1 variables (excluding dimension variables):
    double pdsi[time,lat,lon]   
        longname: Palmer Drought Severity Index
        units: unitless

 3 dimensions:
    time  Size:2013
    lat  Size:88
        longname: latitude
        units: degrees
    lon  Size:114
        longname: longitude
        units: degrees

2 global attributes:
    creation_date: 10-Mar-2015 15:53:08
    Description: Reconstructed PDSI for Europe-Mediterranean Region (OWDA v1.0)

person GIS_newb    schedule 04.04.2018    source источник
comment
P.S. Кажется, что слои по какой-то странной причине не отсортированы по времени. Есть ли способ исправить это?   -  person GIS_newb    schedule 05.04.2018
comment
Я настоятельно рекомендую использовать пакет raster для простой обработки данных с координатной сеткой, включая netcdf. Пакет stars, который сейчас находится в зачаточном состоянии и доступен только на github, также является интересным вариантом.   -  person AF7    schedule 05.04.2018


Ответы (1)


Похоже, ваше изображение нужно повернуть:

rotate <- function(x) t(apply(x, 2, rev))
image(lon, lat, rotate(drought.slice), col = rev(brewer.pal(10, "RdBu")))

работает на меня.

person Amadou Kone    schedule 05.04.2018
comment
Благодарю вас! Это создает карту, но я не уверен, что это карта Европы, какой она должна быть. На севере должна быть Скандинавия, а на юге — Италия. - person GIS_newb; 05.04.2018
comment
Хорошо, это работает, когда я использую функцию поворота ‹- (x) t (apply (x, c (1,2), rev)) Большое спасибо! - person GIS_newb; 10.04.2018