I'm using projectRaster and rasterToPoints in R to get a dataframe with lat/lon coordinates for plotting in ggplot but I seem to be losing the right edge of the map (New Zealand). I'm not sure if this is an issue with my use of projectRaster or rasterToPoints or geom_raster

Here's my initial rasterLayer:

class       : RasterLayer 
band        : 1  (of  4  bands)
dimensions  : 578, 943, 545054  (nrow, ncol, ncell)
resolution  : 46261.72, 46261.72  (x, y)
extent      : -20134202, 23490599, -9583413, 17155861  (xmin, xmax, ymin, 
coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 
+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs 

Which produces: enter image description here

And here's my code to to project and convert to a dataframe:

wgs84 <- "+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs 
b.wgs <- projectRaster(dilmap, crs=CRS(wgs84))
b2.wgs <- rasterToPoints(b.wgs, spatial=TRUE)
b2.df <- as.data.frame(b2.wgs)


enter image description here

EDIT: Added plot of b.wgs

  • Does plotting b.wgs look correct? Also, if I try your code with a dilmap created from random data with the same size/extent etc as your but with random numbers I get "Error: Aesthetics must be either length 1 or the same as the data (362406): x, y, fill" which might actually be a clue. Is the nrows(b2.df) the same as the number of cells in dilmap? Because I think they have to match to make this work. In any case, I'd use tmap package for plotting spatial and avoid ggplot.
    – Spacedman
    Jan 9 '19 at 14:49
  • Aha, b.wgs is indeed clipped. nrows of 430 in b.wgs but 578 in dilmap. Any idea how I correct? Jan 9 '19 at 19:47
  • Even using a reference raster to align nrow, ncol and ncell still not solving the problem of the right edge not plotting. dimensions : 578, 943, 545054 (nrow, ncol, ncell) resolution : 0.3817603, 0.3114187 (x, y) extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 data source : in memory names : dilmap values : 0, 253.858 (min, max) Jan 9 '19 at 20:08
  • Can you edit your Q and show what b.wgs looks like? If it is clipped then all the stuff about making a data frame and using ggplot isn't relevant.
    – Spacedman
    Jan 9 '19 at 20:16
  • Yep, looks like that's the issue. Is there an easy way to correct for the stretching and clipping? Do I need to edit the CRS defined in wgs84 or add an argument to projectRaster()? Sorry for the ignorance, clearly I'm a GIS newbie. Jan 9 '19 at 21:44

Your source raster isn't from -180 to +180 in its coordinate system - it goes beyond +180 and this is being chopped when you transform it. To see:

Create a point at almost +180E on the equator:

> pt = SpatialPoints(coords=matrix(c(179.9,0),ncol=2),proj4string=CRS(wgs84))

Now transform to the coordinate system of your raster to see where it would be:

> coordinates(spTransform(pt, CRS(projection(dilmap))))
     coords.x1     coords.x2
[1,]  20026376 -7.081155e-10

But your raster extent has an x-coordinate max of:

> xmax(dilmap)
[1] 23490599

which is a lot bigger than the X coord of the transformed point. So the right hand portion of your raster is beyond +180E. When transforming the coordinates this actually gets wrapped round to the western hemisphere, and will overlap with the portion of raster on the left of your image, because the left edge of your raster is at 180W, or near enough:

> pt = SpatialPoints(coords=matrix(c(-179.9,0),ncol=2),proj4string=CRS(wgs84))
> coordinates(spTransform(pt, CRS(projection(dilmap))))
     coords.x1     coords.x2
[1,] -20026376 -7.081155e-10
> xmin(dilmap)
[1] -20134202

So there's something fundamentally wrong with your source raster, and you should check carefully you've got the CRS correct or that the data hasn't been manipulated or shifted or otherwise mashed up, or that its extent hasn't been set wrongly.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.