【问题标题】:NetCDF in R: having trouble creating a mapR中的NetCDF:创建地图时遇到问题
【发布时间】:2018-09-14 15:26:22
【问题描述】:

所以我一直在 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 中的错误(经度、纬度、干旱.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)

【问题讨论】:

  • P.S.由于某种奇怪的原因,这些层似乎不是按时间排列的。有没有办法解决这个问题?
  • 我强烈建议使用raster 包来轻松处理网格数据,包括netcdf。 stars 包现在还处于起步阶段,只能从 github 获得,也是一个有趣的选择。

标签: r gis netcdf netcdf4 ncdf4


【解决方案1】:

看起来您的图像需要旋转:

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

为我工作。

【讨论】:

  • 谢谢!这会生成一张地图,但我不确定这是否是欧洲地图,因为它应该是。北边应该是斯堪的纳维亚半岛,南边应该是意大利。
  • 好的,当我使用 "rotate
猜你喜欢
  • 1970-01-01
  • 2018-02-09
  • 1970-01-01
  • 1970-01-01
  • 2023-04-04
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多