Superponer datos de elevación y datos batimétricos en R?

Question

Tengo dos capas que cubren aproximadamente, pero no exactamente, la misma área. Uno es la elevación de la tierra de algunas islas. El otro es la profundidad del agua alrededor de las islas. Puedo resolver los datos verticales entre las capas, así que eso no es un problema.

Me gustaría crear una nueva capa que consiste en la elevación de la tierra y la profundidad del agua circundante. Los problemas:

1. Las extensiones no son exactamente las mismas
2. Las resoluciones son muy diferentes
3. Las relaciones de aspecto son diferentes

¿Hay una manera fácil de crear una nueva capa mirando la capa de mayor resolución (tierra) y si el valor es menor que 75 suponer que es agua, luego mirar la otra capa ráster como la profundidad del agua de lo contrario aceptar la elevación de la tierra como el valor en la nueva capa?

Tierra

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)

Agua

> 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)

---

Pregunta Original

Deseo combinar datos de elevación en una capa ráster con datos batimétricos en otra capa ráster.

Creo que usar raster::overlay es el enfoque correcto aquí. Sin embargo, antes de hacer que funcione, necesito que ambas capas tengan la misma extensión.

Aquí está el código para aislar la parte de tierra en la que estoy interesado.

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

# Convertir a clase RasterLayer
img_raster <- raster(file)

# Transformar RasterLayer a sistema de coordenadas lat/long
elevation <- projectRaster(img_raster, crs="+proj=longlat +datum=WGS84")

# Crear caja delimitadora
# Longitud (xmin, xmax), Latitud (ymin, ymax)
island_rect <- as(raster::extent(-79.4, -79.33, 43.613, 43.635), "SpatialPolygons")
proj4string(island_rect) <- "+proj=longlat +datum=WGS84"

# Recortar solo las islas
to_island <- crop(elevation, island_rect)

Cuando veo la extensión, no es exactamente lo que ingresé:

> 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)

Esto se convierte en un problema porque cuando uso raster::overlay para combinar la elevación de la tierra y los datos batimétricos en una capa ráster overlay da un error porque las extensiones son diferentes: Error in compareRaster(x) : different extent

Suponga que en la capa de tierra cualquier valor menor a 75 es agua. No estoy 100% seguro de cómo usar el comando overlay.

to_island_combo <- raster::overlay(to_island, to_harbour, fun = function(x, y) {
  (x < 75) <- y
  return(x)
})

Answer 1

Primero configuraré algunos datos como los tuyos:

Extensiones y sistemas de coordenadas:

> 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"

Crea raster de tierra:

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

Ingresa algunos valores: esto crea dos bandas debajo del "nivel del mar" en la parte superior e inferior:

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

Ahora prepara un raster de "agua" como el tuyo:

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

Y coloca valores del 1 al N en él para que podamos ver dónde se reemplaza:

> water[]=1:ncell(water)

Los datos están listos, ahora podemos probar cosas.

Los índices de las celdas que queremos reemplazar son:

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

Lo cual parece correcto, son 100,000 ubicaciones en la parte superior y 100,000 en la parte inferior. A continuación, obtén las coordenadas de esos puntos en el sistema de coordenadas del raster de tierra:

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

Estas ubicaciones no están en el sistema de coordenadas del raster de agua, por lo que necesitamos transformarlas:

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

luego obtener los valores de raster de agua en esos puntos:

> ex = extract(water, xyt)

Y reemplaza las ubicaciones en el raster de tierra con los valores muestreados del raster de agua:

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

Esto muestra las dos bandas graduadas que han sido reemplazadas por los valores de agua.

Ten en cuenta que algunos píxeles de tierra por debajo de 75 en mi ejemplo están fuera de la extensión de mi raster de agua y por lo tanto devuelven NA - plot(is.na(land)) mostrará esto.

¿Por qué no envolver todo esto en una función y probarlo?