0

Thanks to this answer I've been able to plot maps from NetCDF files with the correct projection and overlay map. However, when plotting filled contour maps using ggplot, I encounter an error. (see the linked answer for the input file)

#See https://gis.stackexchange.com/questions/120900/plotting-netcdf-file-using-lat-and-lon-contained-in-variables

# First, estimate the extent.
# We start with the lat and long, then project it to the Lambert Conformal projection
library(raster)
inputfile <- "file.nc"

# Grab the lat and lon from the data
lat <- raster(inputfile, varname="xlat")
lon <- raster(inputfile, varname="xlon")

# Convert to points and match the lat and lons
plat <- rasterToPoints(lat)
plon <- rasterToPoints(lon)
lonlat <- cbind(plon[,3], plat[,3])

# Specify the lonlat as spatial points with projection as long/lat
lonlat <- SpatialPoints(lonlat, proj4string = CRS("+proj=longlat +datum=WGS84"))

# Need the rgdal package to project it to the original coordinate system
library("rgdal")

# My best guess at the proj4 string from the information given
mycrs <- CRS("+proj=lcc +lat_1=35.00 +lat_2=51.00 +lat_0=39.00 +lon_0=14.00 +x_0=-6000. +y_0=-6000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs")
plonlat <- spTransform(lonlat, CRSobj = mycrs)
# Take a look
plonlat
extent(plonlat)

# Yay! Now we can properly set the coordinate information for the raster
pr <- raster(inputfile, varname="pr")
# Fix the projection and extent
projection(pr) <- mycrs
extent(pr) <- extent(plonlat)
# Take a look
pr
plot(pr)

# Project to long lat grid
r <- projectRaster(pr, crs=CRS("+proj=longlat +datum=WGS84"))

library(rasterVis)
#Plot with ggplot:
#Add an overlay map
library(maps)
world    <-data.frame(map(plot=FALSE)[c("x","y")])
gplot(r*86400) + 
  geom_tile(aes(fill=value)) + 
  scale_fill_discrete() + 
  geom_path(aes(x,y), data=world) + 
  coord_equal(xlim=c(xmin(r), xmax(r)), ylim=c(ymin(r), ymax(r)))

Error: Continuous value supplied to discrete scale
This is because I want to plot a discrete contour map, with discrete colours, not with a continuous colour scale. Of course it works if using scale_fill_continuous.
I cannot use fill=factor(value) beacuse it will take ages (and GBs of memory) and then return me with a color entry for each single different value in the dataset. I could manually create some "bins" and modify the dataset to fit into those, but I feel this would be a workaround for a much easier solution.

What is the simple elegant solution I'm missing? Shouldn't I use geom_tile (or geom_raster)?

EDIT: Thsi is a pdf example of what I'd like to get: https://copy.com/yMDzEt4i1ELMxpca That plot is created with GrADS. The presence or absence of a color legend is not important, sometimes I use it sometimes I don't.

EDIT2: cut(value) could do but as stated I'd like a ggplot solution and additionally...

...is an example of what it produces. I'd like the labels to be between the colors, not above, like this:

Community
  • 1
  • 1
AF7
  • 3,160
  • 28
  • 63
  • 1
    It is not completely clear how the map should finally look like. Apparently you don't want to have a color for each unique value in column `value` (otherwise you could use `fill = factor(value)`). For me the bin solution makes the most sense, and it should be reasonable fast (`?cut`). – thothal Nov 05 '14 at 14:20
  • I edited the question with a link to an example pdf which represents more or less what I want to obtain. `cut()` is indeed an option (thanks! I did not think about it tbh) but I'll be surprised if there is no ggplot way of dealing with this. – AF7 Nov 05 '14 at 14:46
  • 1
    Maybe this is what you are looking for: http://stackoverflow.com/a/25328524 – Oscar Perpiñán Nov 05 '14 at 16:34
  • @OscarPerpiñán, yes, using `reclassify()` seems a good option. At this point tho I thing `cut()` looks cleaner and simpler. @thothal Your comment led me in the correct direction, which seems to be to use `cut()`. Wuould you like to expand the answer below, and I'll accept it? If not, I'll write a complete answer myself or complete yours in comments. I'll still need a way to copy the colorbar of my last question edit but that maybe is material for another question. – AF7 Nov 06 '14 at 14:46
  • @AdrianoFantini `reclassify` and `cut` are similar functions. Beware that there is a `cut` method defined for `Raster*` objects. – Oscar Perpiñán Nov 07 '14 at 06:55
  • @OscarPerpiñán yes, I'm aware of the fact that there are two `cut()` definitions when using the `raster` package. Am I wrong in supposing that the `raster` function is automatically used if the package is loaded? – AF7 Nov 07 '14 at 09:00
  • 1
    @AdrianoFantini It depends on the object you use with the function. The `raster` method will only be used if you feed it with a `Raster*` object. Read the "Method Selection and Dispatch: Details" section [here](https://stat.ethz.ch/R-manual/R-devel/library/methods/html/Methods.html). – Oscar Perpiñán Nov 07 '14 at 09:28

1 Answers1

2

Based on your edits I would do the following:

r$value.bin <- cut(r$value, c(0:4, 2 * (3:5)))
gplot(r*86400) + 
    geom_tile(aes(fill = value.bin)) + 
    scale_fill_manual(
       values = colorRampPalette(brewer.pal(9, "Greens"))(nlevels(r$value.bin))) + 
    geom_path(aes(x,y), data=world) + 
    coord_equal(xlim=c(xmin(r), xmax(r)), ylim=c(ymin(r), ymax(r)))

Maybe you have to play with the colours in scale_fill_manual and the cutting a bit such that there are spots with an NA color.

thothal
  • 16,690
  • 3
  • 36
  • 71
  • Using `cut()`looks indeed to be the correct way. In your answer there's something wrong, since it's not "world" I want to plot (that's just the map overlay), but "r", but it helped me getting the colors right using `scale_fill_manual()` – AF7 Nov 06 '14 at 14:42
  • You are still missing a ) in the first line and using a "world$..." instead of "r$..." in the fifth ;) EDIT: just tred it and even if fixed it fails: `Error in eval(expr, envir, enclos) : object 'value.bin' not found` even though `r$value.bin` is correctly a RasterLayer – AF7 Nov 06 '14 at 15:08
  • 1
    Well, the problem is that I cannot redo your example as I do not have your data. Thus, I could not test the code. The problem occurs because `gplot` is a special function working with `raster` objects (which in turn are `S4` objects). Thus, `gplot` uses a `raster` object, extracts the relevant data and produces a `ggplot` object. In this transformation step a `data.frame` is produced where there is a column `value`, which needs to be transformed. `g <- gplot(r * 86400)` and `g$data$value.bin <- cut(g$data$value, c(0:4, 2 * (3:5)))` should do the trick. – thothal Nov 06 '14 at 16:27
  • I will try that tomorrow. The data is more or less included above, see the first line of the question. – AF7 Nov 06 '14 at 16:39