I have two layers that cover approximately, but not exactly, the same area. One is land elevation of some islands. The other is the depth of the water around the islands. I can resolve the vertical datums between the layers so that is not an issue.

I would like to create a new layer that consists of the land elevation and the depth of the surrounding water. The problems:

  1. The extents aren't exactly the same
  2. The resolutions are very different
  3. The aspect ratios are different

Is there an easy way to create a new layer by looking at the higher resolution layer (land) and if the value is less than 75 say assume it's water then look to the other raster layer as to the depth of water otherwise accept the land elevation as the value in the new layer?


class       : RasterLayer 
dimensions  : 488, 1129, 550952  (nrow, ncol, ncell)
resolution  : 6.2e-05, 4.5e-05  (x, y)
extent      : -79.39998, -79.32998, 43.61302, 43.63498  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : dem_620483 
values      : 74.89935, 87.08463  (min, max)


> to_harbour
class       : RasterLayer 
dimensions  : 26, 84, 2184  (nrow, ncol, ncell)
resolution  : 0.0008333333, 0.0008333333  (x, y)
extent      : -79.39958, -79.32958, 43.61292, 43.63458  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : in memory
names       : ontario_lld 
values      : -15.9335, 5.764997  (min, max)

Original Question

I wish to combine elevation data in one raster layer with bathymetric data in another raster layer.

I think using raster::overlay is the correct approach here. However before I get that to work I require both layers to have the same extent.

Here is the code to isolate the part of land I am interested in.

file <- "data/GTA_elevation _data/GTA_DEM/6a7a7e71-f502-4336-bba2-364c7eefd950-SW/dem_620483/dem_620483.flt"

# Convert to class RasterLayer
img_raster <- raster(file)

# Transform RasterLayer to lat/long coordinate system
elevation <- projectRaster(img_raster, crs="+proj=longlat +datum=WGS84")

# Create bounding box
# Longitude (xmin, xmax), Latitude (ymin, ymax)
island_rect <- as(raster::extent(-79.4, -79.33, 43.613, 43.635), "SpatialPolygons")
proj4string(island_rect) <- "+proj=longlat +datum=WGS84"

# Crop just the islands
to_island <- crop(elevation, island_rect)

When I look at the extent it is not exactly what I entered:

> to_island
class       : RasterLayer 
dimensions  : 488, 1129, 550952  (nrow, ncol, ncell)
resolution  : 6.2e-05, 4.5e-05  (x, y)
extent      : -79.39998, -79.32998, 43.61302, 43.63498  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : dem_620483 
values      : 74.89935, 87.08463  (min, max)

This becomes an issue because when I use raster::overlay to combine land elevation and bathymetric data into one raster layer overlay gives an error because the extents are different: Error in compareRaster(x) : different extent

Assume that in the land layer anything lower than 75 is water. I'm not 100% sure how to use the overlay command.

to_island_combo <- raster::overlay(to_island, to_harbour, fun = function(x, y) {
  (x < 75) <- y
  • crop will only return whole pixels, so if you apply crop with an extent that is not on a pixel boundary you'll get something truncated (or possibly extended, I'm not sure if it rounds to the nearest pixel boundary). I think you need to show us summary(elevation) and summary(to_harbour) because to use overlay not only do the projections need to be the same but the extents and resolutions do too. – Spacedman Apr 10 at 7:01
  • Thanks Spacedman. I've added info to the question - different extents and resolutions. – ixodid Apr 11 at 0:13
  • One approach that springs to mind is to get the coordinates of all the pixels in the land raster that you want to replace with water values, extract the water value from the water raster, and then put those values into the land raster. Your result is at the resolution of the land raster. – Spacedman Apr 11 at 7:16

First I'll set up some data like yours:

Extents and coordinate systems:

> eland = extent(-79.39998, -79.32998, 43.61302, 43.63498)
> cland = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
> ewater=extent(-79.39958, -79.32958, 43.61292, 43.63458)
> cwater = "+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

Make land raster:

> mland = matrix(0,ncol=1129,nrow=488)
> land = raster(mland, crs=cland)
> extent(land)=eland

Put some values in - this creates two bands below "sea level" at the top and bottom:

> land[] = c(rep(74, 100000),rep(89,350952),rep(72,100000))

Now prep a "water" raster like yours:

> mwater = matrix(0,ncol=84,nrow=26)
> water = raster(mwater, crs=cwater)
> extent(water) = ewater

And put values of 1 to N in it so we can see where it gets replaced:

> water[]=1:ncell(water)

Data is all ready, now we can try things.

The indexes of the cells we want to replace are:

> w = which(land[]<75)
> length(w)
[1] 200000

Which looks correct - that's 100,000 locations at the top and 100,000 at the bottom. Next get the coordinates of those points in the land raster coordinate system:

> xy = SpatialPoints(land)[w]
> proj4string(xy) = projection(land)

These locations aren't in the coordinate system of the water raster so we need to transform them:

> xyt = spTransform(xy, projection(water))

then get the water raster values at those points:

> ex = extract(water, xyt)

And replace the locations in the land raster with the values sample from the water raster:

> land[w] = ex
> plot(land)

water and land

This shows the two graduated bands that have been replaced by the water values.

Note that some of the land pixels below 75 in my example are outside the extent of my water raster and so return NA - plot(is.na(land)) will show this.

Why not wrap this all up in a function and try it?

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.