-1

I'm trying to plot precipitation data which has a 2.5 x 2.5 grid with the country contour on top, the data is available in this link: https://www.esrl.noaa.gov/psd/data/gridded/data.cmap.html "Mean (Enhanced Monthly)"

I was using the answer from: R - Plotting netcdf climate data. However I get an error.

This is what I have done:

library(ncdf4)
ncpath <- "C:/Users/"
ncname <- "precip.mon.mean"
ncfname <- paste(ncpath,ncname,".nc",sep="")
ncin <- nc_open(ncfname)

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

lat <- ncvar_get(ncin, "lat")
nlat <- dim(lat)

dname <-"precip"
ppt_array <- ncvar_get(ncin,dname)
dim(ppt_array)

pres <- ppt_array[ , ,25:444] 
precip <- array(pres, , dim=c(nlon, nlat, 12, ano)) 
prec <- precip[97:115,21:34, ,1:ano] #I just want a piece of the map

Here is where I have the problem:

latlat <- rev(lat)
precipit <- prec[ , ,1,1] %Just to see if it works
lonlon <- lon-180
image(lonlon,latlat,precipit) 
library(maptools)
data(wrld_simpl) 

#however I don't know if this will work to plot just a portion of the map  
plot(wrld_simpl,add=TRUE)

I get several errors, could someone please help?

EDIT: The errors I got were these:

> image(lonlon,latlat,precipit)
Error in image.default(lonlon, latlat, precipit) : 
  increasing 'x' and 'y' values expected
> library(maptools)
> data(wrld_simpl)
> plot(wrld_simpl,add=TRUE)
Error in polypath(x = mcrds[, 1], y = mcrds[, 2], border = border, col = col,  : 
  plot.new has not been called yet
Ann M
  • 239
  • 2
  • 14
  • Can you post all the error messages that you saw? – Tung Jan 11 '18 at 07:17
  • I would strongly suggest to use package `raster` to handle simple netCDFs like this. It is so handy. The upcoming `stars` package will be even better, but it is still in heavy development. – AF7 Jan 11 '18 at 07:54
  • I'll give it a try with the raster, thank you! – Ann M Jan 12 '18 at 05:20

1 Answers1

1

There's several things that need to be fixed:

1) ano does not seem to be defined anywhere. Perhaps it was defined interactively?

precip <- array(pres, , dim=c(nlon, nlat, 12, ano))

2) It appears you intended to add a comment but used an infix operator instead - replace this with a #, like so:

precipit <- prec[ , ,1,1] # Just to see if it works

3) If you want to only have part of the map, you can either ensure that both the lat and lon arrays match the area that you want to show (essentially cropping the world map) or define NAs outside the region you want to highlight (which will appear similar to the map here)

Ming Chia
  • 51
  • 1
  • 4
  • You are right I forgot to copy it here, it's 35. Thank you for the link I'm gonna chek it out. – Ann M Jan 12 '18 at 05:19