How to properly plot projected gridded data in ggplot2?

Question

I've been using using ggplot2 to plot climatic gridded data for years. These are usually projected NetCDF files. Cells are square in model coordinates, but depending on which projection the model uses, might not be so in the real world.

My usual approach is to remap the data first on a suitable regular grid, and then plot. This introduces a small modification to the data, usually this is acceptable.

However, I've decided this is not good enough anymore: I want to plot the projected data directly, without remapping, as other programs (e.g. ncl) can, if I am not mistaken, do, without touching the model output values.

However, I am encountering some issues. I will detail the possible solutions step-by-step below, from the simplest to the most complex, and their problems. Can we overcome them?

EDIT: thanks to @lbusett's answer I got to this nice function that encompasses the solution. If you like it please upvote @lbusett's answer!

Initial setup

#Load packages
library(raster)
library(ggplot2)

#This gives you the starting data, 's'
load(url('https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata'))
#If you cannot download the data, maybe you can try to manually download it from http://s000.tinyupload.com/index.php?file_id=04134338934836605121

#Check the data projection, it's Lambert Conformal Conic
projection(s)
#The data (precipitation) has a 'model' grid (125x125, units are integers from 1 to 125)
#for each point a lat-lon value is also assigned
pr <- s[[1]]
lon <- s[[2]]
lat <- s[[3]]

#Lets get the data into data.frames
#Gridded in model units:
pr_df_basic <- as.data.frame(pr, xy=TRUE)
colnames(pr_df_basic) <- c('lon', 'lat', 'pr')
#Projected points:
pr_df <- data.frame(lat=lat[], lon=lon[], pr=pr[])

We created two dataframes, one with model coordinates, one with real lat-lon cross points (centres) for each model cell.

Optional: using a smaller domain

If you want to more clearly see the shapes of the cells, you can subset the data and extract only a small number of model cells. Just be careful that you might need to adjust point sizes, plot limits and other amenities. You can subset like this and then redo the above code part (minus the load()):

s <- crop(s, extent(c(100,120,30,50)))

If you want to fully understand the problem maybe you want to try both the big and the small domain. the code is identical, only the point sizes and map limits change. The values below are for the big complete domain. Ok, now let's plot!

Start with tiles

The most obvious solution is to use tiles. Let's try.

my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())
my_cols <- scale_color_distiller(palette='Spectral')
my_fill <- scale_fill_distiller(palette='Spectral')

#Really unprojected square plot:
ggplot(pr_df_basic, aes(y=lat, x=lon, fill=pr)) + geom_tile() + my_theme + my_fill

And this is the result:

Ok, now something more advanced: we use real LAT-LON, using square tiles

ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.2, height=1.2) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat)) #the result is weird boxes...

Ok, but those are not the real model squares, this is a hack. Also, model boxes diverge at the top of the domain, and are all oriented in the same way. Not nice. Let's project the squares themselves, even though we already know this is not the right thing to do... maybe it looks good.

#This takes a while, maybe you can trust me with the result
ggplot(pr_df, aes(y=lat, x=lon, fill=pr)) + geom_tile(width=1.5, height=1.5) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))

First of all, this takes a lot of time. Not acceptable. Also, again: these are not correct model cells.

Try with points, not tiles

Maybe we can use round or square points instead of tiles, and project them too!

#Basic 'unprojected' point plot
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) + geom_point(size=2) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_cols + my_theme +
    coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat))

We can use square points... and project! We get closer, even though we know that it is still not correct.

#In the following plot pointsize, xlim and ylim were manually set. Setting the wrong values leads to bad results.
#Also the lambert projection values were tired and guessed from the model CRS
ggplot(pr_df, aes(y=lat, x=lon, color=pr)) +
    geom_point(size=2, shape=15) +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_cols +
    coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75))

Decent results, but not totally automatic and plotting points is not good enough. I want real model cells, with their shape, mutated by the projection!

Polygons, maybe?

So as you can see I am after a way of correctly plotting the model boxes as projected in the correct shape and position. Of course the model boxes, which are squares in the model, once projected become shapes which are not regular anymore. So maybe I can use polygons, and project them? I tried to use rasterToPolygons and fortify and follow this post, but have failed to do so. I have tried this:

pr2poly <- rasterToPolygons(pr)
#http://mazamascience.com/WorkingWithData/?p=1494
pr2poly@data$id <- rownames(pr2poly@data)
tmp <- fortify(pr2poly, region = "id")
tmp2 <- merge(tmp, pr2poly@data, by = "id")
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill

Ok, let's try to substitute lat-lons...

tmp2$long <- lon[]
tmp2$lat <- lat[]
#Mh, does not work! See below:
ggplot(tmp2, aes(x=long, y=lat, group = group, fill=Total.precipitation.flux)) + geom_polygon() + my_fill

(sorry that I changed the color scale in the plots)

Mmmmh, not even worth trying with a projection. Maybe I should try to compute the lat-lons of the model cell corners, and create polygons for that, and reproject that?

Conclusion

1. I want to plot projected model data on its native grid, but I was not able to do so. Using tiles is incorrect, using points is hackish, and using polygons seems not to work for unknown reasons.
2. When projecting via coord_map(), gridlines and axis labels are wrong. This makes the projected ggplots unusable for publications.

Winne

Have you tried ggmap? – Boarfish

@Boarfish I remember having use it in the past, but not recently. I forgot it existed. I'll look into it. Anything of particular interest to point out? – Winne

I couldn't load your data; are you sure the URL is correct? – Mayce

@Mayce I changed the donwload provider, please try again. This should work more reliably. – Winne

Your data is called from disk: "/home/adriano/EURO-CORDEX_59_tas_pr_prc_monmean_1998-2000.nc" So you save an object with no real data but an url into your machine. Instead RData can you save in a raster layer? – Zeldazelde

@Zeldazelde My bad, I did not do a readAll() before saving. It should work now, can you test it? – Winne

Hi. I think there's something strange with your dataset. Projection of s results as "+proj=lcc +lat_1=30. +lat_2=65. +lat_0=48. +lon_0=9.75 +x_0=-25000. +y_0=-25000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs", which is a "metric" projection. However, your coordinates are lat-long, and the resolution of the raster is 1. Where did you get the dataset ? Did you do any "spatial" elaborations on it ? – Lionellionello

Can't you just subsample the raster to a finer resolution prior to reprojecting? – Apogamy

@Winne Do yo need to plot using LCC projection? Can't you just change the projection on the raster data? – Evildoer

@LoBu Nice find. This is the data directly as it comes out of a model, and the projection string is given by the model itself. I have no idea why it says +units=m and I am no expert on projections, maybe there is a reason. Weird. Still, it seems to work if I reproject stuff. – Winne

@JacobSocolar The whole point of the question is that I do not want to touch the input data at all prior to plotting. If this was not a constraint, I could just remap with a proper weighted technique and be content with it, it does not affect the data so much in most cases. – Winne

@Evildoer The whole point of the question, as stated before, is that I do not want to introduce any modification to the data values. Hence I cannot remap, which is what I usually do when doing this kind of plots. – Winne

It seems to me that while the lat/long coordinates in s[[1]] and s[[2]] are correct, the "formatting" as a raster stack is completely off. In fact, if you try to visualize one of the layers using mapview, you end up with a completely wrong coordinates system. – Lionellionello

@LoBu I do not understand. Of course you end up with a wrong coordinate system, as the x-y coordinates are just model coordinates and go from 1 to 125. This is the point of having separate variables for lat-lon mapping – Winne

@Winne my point is just that then you shouldn't save the dataset as a RasterStack object. This may lead to confusion, since for example standard raster processing routines such as projectRaster are not going to work properly. Saving it as a SpatialPixelsDataFrame would be better, in my opinion. – Lionellionello

@LoBu the original dataset is netCDF and I usually read that with raster. If a projection is necessary inside R, I follow these excellent instructions: gis.stackexchange.com/questions/120900/… – Winne

But the plot will never show exactly the original data. If nothing else, its limited by the resolution of the graphics device. Is the question just academic, or is there an underlying efficiency issue, or... – Apogamy

@JacobSocolar Well if you plot polygons to PDF then the resolution of the graphics device is non-important. The main reason why I want to do this is that I want to reach the same level of precision other softwares have. I love plotting with R and ggplot and knowing that the pretty terrible NCL can do more precise plots is quite annoying :) Other NetCDF plotting softwares (such as GrADS and I believe Ferret), on the other hand, just automatically remap on their internal high-res LAT-LON grid without letting the user know it. These plots are considered acceptable, but I want to do better. – Winne

Answer 1

After digging a bit more, it seems that your model is based on a 50Km regular grid in "lambert conical" projection. However, the coordinates you have in the netcdf are lat-lon WGS84 coordinates of the center of the "cells".

Given this, a simpler approach is to reconstruct the cells in the original projection and then plot the polygons after converting to an sf object, eventually after reprojection. Something like this should work (notice that you'll need to install the devel version of ggplot2 from github for it to work):

load(url('https://files.fm/down.php?i=kew5pxw7&n=loadme.Rdata'))
library(raster)
library(sf)
library(tidyverse)
library(maps)
devtools::install_github("hadley/ggplot2")

#   ____________________________________________________________________________
#   Transform original data to a SpatialPointsDataFrame in 4326 proj        ####

coords = data.frame(lat = values(s[[2]]), lon = values(s[[3]]))
spPoints <- SpatialPointsDataFrame(coords, 
                                   data = data.frame(data = values(s[[1]])), 
                                   proj4string = CRS("+init=epsg:4326"))

#   ____________________________________________________________________________
#   Convert back the lat-lon coordinates of the points to the original      ###
#   projection of the model (lcc), then convert the points to polygons in lcc
#   projection and convert to an `sf` object to facilitate plotting

orig_grid = spTransform(spPoints, projection(s))
polys = as(SpatialPixelsDataFrame(orig_grid, orig_grid@data, tolerance = 0.149842),"SpatialPolygonsDataFrame")
polys_sf = as(polys, "sf")
points_sf = as(orig_grid, "sf")

#   ____________________________________________________________________________
#   Plot using ggplot - note that now you can reproject on the fly to any    ###
#   projection using `coord_sf`

# Plot in original  projection (note that in this case the cells are squared): 
my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())

ggplot(polys_sf) +
  geom_sf(aes(fill = data)) + 
  scale_fill_distiller(palette='Spectral') +
  ggtitle("Precipitations") +
  coord_sf() + 
  my_theme 

# Now Plot in WGS84 latlon projection and add borders: 

ggplot(polys_sf) +
  geom_sf(aes(fill = data)) + 
  scale_fill_distiller(palette='Spectral') +
  ggtitle("Precipitations")  +
  borders('world', colour='black')+
  coord_sf(crs = st_crs(4326), xlim = c(-60, 80), ylim = c(15, 75))+
  my_theme 

To add borders also wen plotting in the original projection, however, you'll have to provide the loygon boundaries as an sf object. Borrowing from here:

Converting a "map" object to a "SpatialPolygon" object

Something like this will work:

library(maptools)
borders  <- map("world", fill = T, plot = F)
IDs      <- seq(1,1627,1)
borders  <- map2SpatialPolygons(borders, IDs=borders$names, 
                               proj4string=CRS("+proj=longlat +datum=WGS84")) %>% 
            as("sf")

ggplot(polys_sf) +
  geom_sf(aes(fill = data), color = "transparent") + 
  geom_sf(data = borders, fill = "transparent", color = "black") +
  scale_fill_distiller(palette='Spectral') +
  ggtitle("Precipitations") +
  coord_sf(crs = st_crs(projection(s)), 
           xlim = st_bbox(polys_sf)[c(1,3)],
           ylim = st_bbox(polys_sf)[c(2,4)]) +
  my_theme

As a sidenote, now that we "recovered" the correct spatial reference, it is also possible to build a correct raster dataset. For example:

r <- s[[1]]
extent(r) <- extent(orig_grid) + 50000

will give you a proper raster in r:

r
class       : RasterLayer 
band        : 1  (of  36  bands)
dimensions  : 125, 125, 15625  (nrow, ncol, ncell)
resolution  : 50000, 50000  (x, y)
extent      : -3150000, 3100000, -3150000, 3100000  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=30. +lat_2=65. +lat_0=48. +lon_0=9.75 +x_0=-25000. +y_0=-25000. +ellps=sphere +a=6371229. +b=6371229. +units=m +no_defs 
data source : in memory
names       : Total.precipitation.flux 
values      : 0, 0.0002373317  (min, max)
z-value     : 1998-01-16 10:30:00 
zvar        : pr 

See that now the resolution is 50Km, and the extent is in metric coordinates. You could thus plot/work with r using functions for raster data, such as:

library(rasterVis)
gplot(r) + geom_tile(aes(fill = value)) + 
  scale_fill_distiller(palette="Spectral", na.value = "transparent") +
  my_theme  

library(mapview)
mapview(r, legend = TRUE)  

Lionellionello

object 'transf' not found. Also, how did you get the 'tolerance'? And why is it necessary to perform a remap to epsg:4326? Does this introduce an (albeit small) error? – Winne

ops, sorry. I edited the post. Concerning the tolerance, it is required beacuase to create a SpatialPixelsDataFrame you need a regular grid, and the "reprojection" of the cells leads probably to some rounding error. If you call SpatialPixelsDataFrame(orig_grid, orig_grid@data) without the tolerance argument, a "suggested" tolerance is given. Concerning the projection: No, converting to 4326 is not needed. It was just to show you that you can easily reproject to any CRS. – Lionellionello

oh, sorry. I probably misinterpreted your question about the "remap" to epsg:4326. When creating the SpatialPixelsDataFrame I am not remapping. I just assign to the coordinates their correct reference system, which is 4326. The fact that that is the "proper" coordinate system is testified by the fact that when you transform to "lcc" projection you obtain a regular grid (you can see that by issuing head(coordinates(orig_grid)) – Lionellionello

Yeah, thanks for clarifying that. It was important for me that if cell N had a value X, this did not get changed in the plot, even slightly. Your answer looks awesome. Let me explore it for a bit before accepting and rewarding it. – Winne

I am having problems plotting with geom_sf() only, on the model grid, and also plotting the borders with borders() from what I gather all layers (including the borders) should be automatically transposed to the input data CRS (just like when using coord_map()), but it does not seem to work. I should get a square plot like yours, with borders on top. I only get a plot like yours, with a small black dot in the middle. What is wrong? – Winne

that's because borders("world",..... is not an sf object, so that evidently the "reprojection" is not applied. See the updated answer for a workaround. – Lionellionello

Thanks. I filed a ggplot github issue about this, as this seems to be a bug. github.com/tidyverse/ggplot2/issues/2118 – Winne

By the way to get borders you can just do borders <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE)) – Winne

I also have another question: why extent(orig_grid) + 50000? Why not 25k, half the grid box size? Why adding anything at all? – Winne

You need to extend the extent by half a pixel wrt the points because the extent of a raster object is computed on corners of the cells. I use 50000 because when "summing" to a extent object, the extent is expanded by half the value in all directions (I found that out by trial and error). – Lionellionello

Excellent. I suggest this way of calculating the 50000: extog <- extent(orig_grid) ; target.res <- (extog[2] - extog[1]) / (nrow(r) -1) – Winne

Good idea. I will incorporate this in the answer. Also: in the function you added to your question, consider you can remove the line points_sf = as(orig_grid, "sf"). It was a leftover from testing if the cells' polygons where properly aligned with points. – Lionellionello

Answer 2

"Zooming in" to see the points that are the cell centres. You can see they are in rectangular grid.

I calculated the vertices of the polygons as follows.

• Convert 125x125 latitudes & longitudes to a matrix

• Initialize 126x126 matrix for cell vertices (corners).

• Calculate cell vertices as being mean position of each 2x2 group of points.

• Add cell vertices for edges & corners (assume cell width & height is equal to width & height of adjacent cells).

• Generate data.frame with each cell having four vertices, so we end up with 4x125x125 rows.

Code becomes

 pr <- s[[1]]
lon <- s[[2]]
lat <- s[[3]]

#Lets get the data into data.frames
#Gridded in model units:
#Projected points:
lat_m <- as.matrix(lat)
lon_m <- as.matrix(lon)
pr_m <- as.matrix(pr)

#Initialize emptry matrix for vertices
lat_mv <- matrix(,nrow = 126,ncol = 126)
lon_mv <- matrix(,nrow = 126,ncol = 126)


#Calculate centre of each set of (2x2) points to use as vertices
lat_mv[2:125,2:125] <- (lat_m[1:124,1:124] + lat_m[2:125,1:124] + lat_m[2:125,2:125] + lat_m[1:124,2:125])/4
lon_mv[2:125,2:125] <- (lon_m[1:124,1:124] + lon_m[2:125,1:124] + lon_m[2:125,2:125] + lon_m[1:124,2:125])/4

#Top edge
lat_mv[1,2:125] <- lat_mv[2,2:125] - (lat_mv[3,2:125] - lat_mv[2,2:125])
lon_mv[1,2:125] <- lon_mv[2,2:125] - (lon_mv[3,2:125] - lon_mv[2,2:125])

#Bottom Edge
lat_mv[126,2:125] <- lat_mv[125,2:125] + (lat_mv[125,2:125] - lat_mv[124,2:125])
lon_mv[126,2:125] <- lon_mv[125,2:125] + (lon_mv[125,2:125] - lon_mv[124,2:125])

#Left Edge
lat_mv[2:125,1] <- lat_mv[2:125,2] + (lat_mv[2:125,2] - lat_mv[2:125,3])
lon_mv[2:125,1] <- lon_mv[2:125,2] + (lon_mv[2:125,2] - lon_mv[2:125,3])

#Right Edge
lat_mv[2:125,126] <- lat_mv[2:125,125] + (lat_mv[2:125,125] - lat_mv[2:125,124])
lon_mv[2:125,126] <- lon_mv[2:125,125] + (lon_mv[2:125,125] - lon_mv[2:125,124])

#Corners
lat_mv[c(1,126),1] <- lat_mv[c(1,126),2] + (lat_mv[c(1,126),2] - lat_mv[c(1,126),3])
lon_mv[c(1,126),1] <- lon_mv[c(1,126),2] + (lon_mv[c(1,126),2] - lon_mv[c(1,126),3])

lat_mv[c(1,126),126] <- lat_mv[c(1,126),125] + (lat_mv[c(1,126),125] - lat_mv[c(1,126),124])
lon_mv[c(1,126),126] <- lon_mv[c(1,126),125] + (lon_mv[c(1,126),125] - lon_mv[c(1,126),124])


pr_df_orig <- data.frame(lat=lat[], lon=lon[], pr=pr[])

pr_df <- data.frame(lat=as.vector(lat_mv[1:125,1:125]), lon=as.vector(lon_mv[1:125,1:125]), pr=as.vector(pr_m))
pr_df$id <- row.names(pr_df)

pr_df <- rbind(pr_df,
               data.frame(lat=as.vector(lat_mv[1:125,2:126]), lon=as.vector(lon_mv[1:125,2:126]), pr = pr_df$pr, id = pr_df$id),
               data.frame(lat=as.vector(lat_mv[2:126,2:126]), lon=as.vector(lon_mv[2:126,2:126]), pr = pr_df$pr, id = pr_df$id),
               data.frame(lat=as.vector(lat_mv[2:126,1:125]), lon=as.vector(lon_mv[2:126,1:125]), pr = pr_df$pr, id= pr_df$id))

Same zoomed image with polygon cells

Labels Fix

ewbrks <- seq(-180,180,20)
nsbrks <- seq(-90,90,10)
ewlbls <- unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste(abs(x), "°W"), ifelse(x > 0, paste(abs(x), "°E"),x))))
nslbls <- unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste(abs(x), "°S"), ifelse(x > 0, paste(abs(x), "°N"),x))))

Replacing geom_tile & geom_point with geom_polygon

ggplot(pr_df, aes(y=lat, x=lon, fill=pr, group = id)) + geom_polygon() +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_quickmap(xlim=range(pr_df$lon), ylim=range(pr_df$lat)) +
scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
scale_y_continuous(breaks = nsbrks, labels = nslbls, expand = c(0, 0)) + 
labs(x = "Longitude", y = "Latitude")

ggplot(pr_df, aes(y=lat, x=lon, fill=pr, group = id)) + geom_polygon() +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 39), ylim=c(19, 75)) +
scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
scale_y_continuous(breaks = nsbrks, labels = nslbls, expand = c(0, 0)) + 
labs(x = "Longitude", y = "Latitude")

Edit - work around for axis ticks

I've been unable to find any quick solution for the grid lines and labels for the latitude. There is probably an R package out there somewhere which will solve your problem with far less code!

Manually setting nsbreaks required and creating data.frame

ewbrks <- seq(-180,180,20)
nsbrks <- c(20,30,40,50,60,70)
nsbrks_posn <- c(-16,-17,-16,-15,-14.5,-13)
ewlbls <- unlist(lapply(ewbrks, function(x) ifelse(x < 0, paste0(abs(x), "° W"), ifelse(x > 0, paste0(abs(x), "° E"),x))))
nslbls <- unlist(lapply(nsbrks, function(x) ifelse(x < 0, paste0(abs(x), "° S"), ifelse(x > 0, paste0(abs(x), "° N"),x))))
latsdf <- data.frame(lon = rep(c(-100,100),length(nsbrks)), lat = rep(nsbrks, each =2), label = rep(nslbls, each =2), posn = rep(nsbrks_posn, each =2))

Remove the y axis tick labels and corresponding gridlines and then add back in "manually" using geom_line and geom_text

ggplot(pr_df, aes(y=lat, x=lon, fill=pr, group = id)) + geom_polygon() +
    borders('world', xlim=range(pr_df$lon), ylim=range(pr_df$lat), colour='black') + my_theme + my_fill +
    coord_map('lambert', lat0=30, lat1=65, xlim=c(-20, 40), ylim=c(19, 75)) +
    scale_x_continuous(breaks = ewbrks, labels = ewlbls, expand = c(0, 0)) +
    scale_y_continuous(expand = c(0, 0), breaks = NULL) + 
    geom_line(data = latsdf, aes(x=lon, y=lat, group = lat), colour = "white", size = 0.5, inherit.aes = FALSE) +
    geom_text(data = latsdf, aes(x = posn, y = (lat-1), label = label), angle = -13, size = 4, inherit.aes = FALSE) +
    labs(x = "Longitude", y = "Latitude") +
    theme( axis.text.y=element_blank(),axis.ticks.y=element_blank())

Stouthearted

I like the approach, but I think there is something wrong with the result, since the spatial patterns in your output map are very different from the original. – Lionellionello

I noticed that. I will verify, but believe it could be down to different data being downloaded with load(url.. or maybe original plots having overlapping tiles or points that lost some of the information. – Stouthearted

I see what you mean. I can't fix it just now, sorry, as I'm at work. But I'm sure it will be fairly straightforward. Probably just a case of doing a matrix transposition. – Stouthearted

I'll suggest a quick work around, that may work, but I'm unable to test at the moment. pr_m <- as.matrix(pr). Then in data.frame use pr = as.vector(pr_m) I'm thinking that treating pr exactly the same as lat and lon will ensure they don't get swapped. – Stouthearted

Will do, either this afternoon or evening (Cyprus time). – Stouthearted

Let us continue this discussion in chat. – Stouthearted

I've fixed the bug. – Stouthearted

@JeremyVoisey awesome. Now the only thing missing would be fixing the ggplot problem with the axis labels and grid lines. I think the way to go would be converting the polygons to sf type and using geom_sf() instead of geom_polygon(). geom_sf() uses st_graticule() and should work much better. I'll look into it but not soon unfortunately. – Winne

I've amended the labels, thanks to an answer found elsewhere #33302924 – Stouthearted

@JeremyVoisey That's an answer to a question I posed, I know how top fix the label that way :). It was not what I meant, I should have been more precise, sorry. Look at the last plot. The labels are misplaced compared to the grid lines. 40N, for example, does not match the 40N line. Same for -20E. Another thing is that the 60N line for some reason does not extend to the borders of the plot, too. – Winne

I see what you mean. The labels are actually horizontally aligned with the intercepts with the -20E Longitude line. Which makes sense in one way, as this is the specified limit on that axis. But obviously doesn't work for a Lambert projection. – Stouthearted

Well, it could be sufficient to "manually" remove the first and last break label if plotting in Lambert, then. By the way: there's an error in the reassignemnt of labels: labels east of Greenwhich should be 10 E, 20 E, etcetera. – Lionellionello

@LoBu Mixing up E & W, that's embarrassing. I got the code from stackoverflow, I assumed it was reliable, my bad. – Stouthearted

I've added a fix for the grid lines (OK) and labels (needs a bit more work). I think @LoBu has come up with a better answer by some way! – Stouthearted

@Jeremy, still very nice effort, I learned something. If I could award both, I would. Thanks. – Winne

Thanks, @LoBu is the clear winner and had my upvote. I've also learnt alot. – Stouthearted