I'm working with raster images. I need to add the european coastline over this plot(fdiff_north) in order to get a more usable image`: plot(fdiff)

I got this plot by using these codes:

ext <- c(4e+06, 5200000, 2250000, 2800000)
fdiff_north <- crop(fdiff, ext)
cl <- colorRampPalette(c("red","white","blue"))(100)
plot(fdiff_north, col=cl)


class      : RasterLayer 
dimensions : 5500, 12000, 6.6e+07  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 4e+06, 5200000, 2250000, 2800000  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs 
source     : memory
names      : layer 
values     : -1, 1  (min, max)

thats the summary(fdiff_north):

Min.       -1
1st Qu.     0
Median      0
3rd Qu.     0
Max.        1
NA's        0

I found the ne_coastline function but I can't figure how to use it. These are the packages I'm using to do that:

    > install.packages("rnaturalearthdata") 
    > library("rnaturalearthdata")
    > install.packages("devtools") 
    > devtools::install_github("ropensci/rnaturalearthhires")
    > library("rnaturalearthhires")

On internet I found that example code but I really have non idea how to use it:

if (requireNamespace("rnaturalearthdata")) {
   sldf_coast <- ne_coastline()

   if (require(sp)) {
  • 1
    Short answer: looks like your data isn't in lat-long so you need to transform your outlines to its coordinate system (not vice-versa or you'll warp your neat grid raster)
    – Spacedman
    Commented Feb 4, 2022 at 16:36
  • Can you edit your Q and show summary(fdiff) so we can get a handle on the projection?
    – Spacedman
    Commented Feb 4, 2022 at 20:58
  • I did it. I'm trying to understand how to transform my outlines to its coordinate system.. thanks Commented Feb 5, 2022 at 11:52
  • Ah, maybe not summary, maybe print, some way that we can see the coordinate system (not that it really matters because we can get it with functions but its useful anyway...)
    – Spacedman
    Commented Feb 5, 2022 at 12:01
  • 1
    Saying "this is not working" doesn't help - you should always explain how something isn't working by showing the error message or explaining why the output isn't right. You can always edit your Q and add extra info and error messages.
    – Spacedman
    Commented Feb 7, 2022 at 19:12

1 Answer 1


First make a reproducible sample data set. You should do something like this in every Q here so we can quickly get to the problem without having to make sample data ourselves. I've reduced the size by a factor of 10x10 but kept the extent and coordinate system:

fdiff_north = raster(matrix(1:(550*1200), 550,1200))
extent(fdiff_north) = c(4e+06, 5200000, 2250000, 2800000)
crs(fdiff_north) = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"

Now get the coastline, and transform it to the coordinate projection of the raster above:

sldf_coast = rnaturalearth::ne_coastline()
sldf_prj = spTransform(sldf_coast, projection(fdiff_north))

You'll probably get warnings. These are generally okay to ignore at this point, and on the scale of your data probably safe to ignore. Don't confuse warnings with errors and think things "dont work".

Now plot the raster and overlay the coastline:

plot(sldf_prj, add=TRUE)

enter image description here

Note its better to use sf packages and classes, you can do that with:

coast = rnaturalearth::ne_coastline(returnclass="sf")
coast_t = st_transform(coast, crs(fdiff_north))
plot(st_geometry(coast_t), add=TRUE)

And you get the same map but probably bit quicker and with less memory.

I'll also show the other operation, transforming the raster to the coast coordinate system:

fdiff_trans = projectRaster(fdiff_north, crs=projection(sldf_coast))
plot(sldf_coast, add=TRUE)

enter image description here

Notice how the output grid is warped to a non-rectangular shape. This means the values are interpolated and its a non-perfectly-reversible transformation - some information is lost. That's why its preferable to warp the lines (which are all defined by coordinate points and so transformations are perfectly reversible) and keep the raster in its own coordinate system. Only if you want to relate to rasters in different coordinate systems might you need to transform one to another.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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