NetCDF in R: having trouble creating a map

Question

So I've been using the following tutorial in R, to get accustomed to NetCDF (I use the netcdf4 package): http://geog.uoregon.edu/bartlein/courses/geog607/Rmd/netCDF_01.htm

However, I am not able to create a nice image for any particular year, in my dataset, which can be found here: https://www.ncdc.noaa.gov/paleo-search/study/19419

Here is my code, in which I attempt to get a map for the year 100 AD:

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")))

In particular, I get the following error message when I try to run the last line:

Error in image.default(lon, lat, drought.slice, col = rev(brewer.pal(10, : dimensions of z are not length(x)(-1) times length(y)(-1)

Could someone please explain what I'm doing wrong?

By the way, to give you some idea about my data's structure, here you go:

>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. It seems that the layers aren't arranged by time, for some strange reason. Is there a way to remedy this? — GIS_newb, Apr 04 '18 at 23:55
I would strongly suggest to use the `raster` package for easy handling of gridded data, including netcdf. The `stars` package, in its infancy now and available only from github, is also an interesting option. — AF7, Apr 05 '18 at 08:36

score 2 · Answer 1 · answered Apr 05 '18 at 00:04

2

It looks like your image needs to be rotated:

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

works for me.

answered Apr 05 '18 at 00:04

Amadou Kone

907
11
21

Thank you! This produces a map, but I'm not sure if this is a map of Europe, as it is supposed to be. It should have Scandinavia in the north, and Italy in the south. – GIS_newb Apr 05 '18 at 12:30
Okay, it works when I use "rotate <- function(x) t(apply(x, c(1,2), rev))" Thank you so much! – GIS_newb Apr 10 '18 at 18:34

NetCDF in R: having trouble creating a map

1 Answers1