How to merge land cover map of raster data using R, terra

ghz 8months ago ⋅ 112 views

I would like to merge land cover maps into one, and I found the code below:

tif_files <- list.files(path = "LandCover",
                          pattern = ".tif$",
                          full.names = TRUE,
                          ignore.case = TRUE
  )
  tif_rasts <- list()
  for (i in 1:length(tif_files)){
    # print(tif_file_name[i])
    tif_rasts <- c(tif_rasts,  rast(tif_files[i]))
  }
  
  
  tif_mosaic <- do.call(mosaic, c(tif_rasts, fun="max"))
  plot(tif_mosaic)
  
  writeRaster(tif_mosaic,
              filename = "land_cover.tif"),
              overwrite = TRUE
  )

However, I found the final land cover data are into 3 bands (layers), and these layers have the same and different values. I would like to know how to combine them into a single layer in order to identify the land use type represented by each value.

unique(tif_mosaic)
   CHI_1_LULC_1 CHI_1_LULC_2 CHI_1_LULC_3
2             0            0            0
3             0            0          255
4            34          139           34
5           128            0            0
6           128          236          104
7           133          133          133
8           255            0            0
9           255            0          192
10          255          193           37
11          255          255          255

Answers

It seems like your land cover data has been merged into a multi-layer raster, where each layer represents a different aspect of land cover. To combine these layers into a single layer that represents the dominant land use type for each pixel, you can use the overlay function from the raster package in R.

Here's how you can do it:

library(raster)

# Read the multi-layer raster
tif_mosaic <- brick("land_cover.tif")

# Combine the layers by selecting the maximum value for each pixel
combined_layer <- overlay(tif_mosaic, fun=max)

# Plot the combined layer
plot(combined_layer)

# Write the combined layer to a new raster file
writeRaster(combined_layer, filename = "combined_land_cover.tif", overwrite = TRUE)

In this code:

  • We use brick to read the multi-layer raster created by merging the land cover maps.
  • We use overlay with the max function to combine the layers. The max function selects the maximum value for each pixel across all layers.
  • We plot the combined layer to visualize the result.
  • Finally, we use writeRaster to save the combined layer as a new raster file.

This will create a new raster file (combined_land_cover.tif) where each pixel value represents the dominant land use type at that location.

#r