1

I want to calculate area for each different class of landuse for raster image as shown below:-

enter image description here

The dimension of the raster image is as follows:

class      : RasterLayer 
dimensions : 16366, 29848, 488492368  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 96488.54, 991928.5, 2893376, 3384356  (xmin, xmax, ymin, 
ymax)
crs        : +proj=utm +zone=44 +datum=WGS84 +units=m +no_defs 
source     : LULC10.tif 
names      : LULC10 
values     : 0, 255  (min, max)

I want to calculate area for each class as shown in figure above. How can this be done? Please help me on it. How is this possible?

Data source: https://download.hermes.com.np/land-cover-map-of-nepal-2010/

2
  • Hey, focus in one question at time: 1- calculate are by class 2- reclass (you missed class 10) 3- change symbology Otherwise, the question will be closed
    – aldo_tapia
    Commented Jun 13, 2022 at 0:05
  • @aldo_tapia,thanks for the information. I now have tried to be more more specific to my question. Please look through it and help me on this. Thank you.
    – Walker
    Commented Jun 13, 2022 at 1:47

2 Answers 2

5

Read the file and set the "active category"

library(terra)
r <- rast("np_lc_2000_v2f.tif")
activeCat(r) <- "Class_Name"
r
#class       : SpatRaster 
#dimensions  : 15985, 29459, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : 96518.54, 980288.5, 2904806, 3384356  (xmin, xmax, ymin, ymax)
#coord. ref. : WGS_1984_Transverse_Mercator 
#source      : np_lc_2000_v2f.tif 
#color table : 1 
#categories  : Red, Green, Blue, Opacity, Class_Name 
#name        :    Class_Name 
#min value   :        Forest 
#max value   : Built-up area 

The naïve approach is to assume that the resolution is constant, in this case 900 m2 or 0.0009 km2.

f <- freq(r)
f$area <- f$count * 0.0009
f$percent = round(100 * f$area / sum(f$area), 1)
f
#  layer            value    count       area
#1     1           Forest 68313326 61481.9934
#2     1        Shrubland  3851371  3466.2339
#3     1        Grassland 18606350 16745.7150
#4     1 Agriculture area 45563649 41007.2841
#5     1      Barren area 15615609 14054.0481
#6     1       Water body   805826   725.2434
#7     1     Snow/glacier 10810771  9729.6939
#8     1    Built-up area   470013   423.0117

To account for varying cell sizes you can compute the actual size of each cell. You have 16366 x 29848 = 470902115 cells and that would take a long time to compute (in fact it will probably fail). But we can approximate the values by first aggregating and then resampling. And then use "zonal" to get the areas.

Compute cell size for many cells work-around

x <- aggregate(rast(r), 100)
a <- cellSize(x, unit="km") / 10000
b <- resample(a, r)
minmax(a)
#            area
#min 0.0008956130
#max 0.0009007204

That is, cell size varies between 0.896 and 0.901 m2. Not very much.

Now use zonal

z <- zonal(b, r, sum, na.rm=TRUE)
z
#        Class_Name       area
#1           Forest 61454.0155
#2        Shrubland  3465.0476
#3        Grassland 16740.2584
#4 Agriculture area 40990.4952
#5      Barren area 14055.4962
#6       Water body   724.8611
#7     Snow/glacier  9730.7667
#8    Built-up area   422.9286

The difference is very small, less than 1‰, so the naïve approach is good enough in this case

(1000 * (f$area - z$area) / f$area)  |> round(2)
#[1]  0.46  0.34  0.33  0.41 -0.10  0.53 -0.11  0.20
2
  • Using your code as:- 'library(terra) RupandehiLULC2010 <- rast('C:/Users/lenovo/Desktop/EM/RupandehiLULC.tif') activeCat(RupandehiLULC2010) <- "Class_Name" RupandehiLULC2010' shows an error :- Error: [activeCat] layer is not categorical. Why is this occurring??
    – Walker
    Commented Jun 23, 2022 at 0:25
  • You need to unzip all the files that are in the zip file. Not only the .tif file. Commented Jun 23, 2022 at 15:47
0

I recommend the use of dpylr package for data manipulation out of raster or terra environment. For this case, I also changed the layer name (in the case the name you have is different than my example).

library(dplyr)

names(LULC) <- 'class_2010'

Then, count the pixels by class and multiply it for pixel area. The output is the same than CRS unit (meter in this case), you can convert it later as you please:

result <- as.data.frame(LULC) %>%
  group_by(class_2010) %>% 
  summarise(pixels = n()) %>% 
  mutate(`area m^2` = pixels * 30 * 30, #multiply by pixel area
         `area km^2` = round(`area m^2`/1e6,2)) 

result
# # A tibble: 9 × 4
#   class_2010    pixels   `area m^2` `area km^2`
#        <int>     <int>        <dbl>       <dbl>
# 1          1  68931804  62038623600      62039.
# 2          2   3811656   3430490400       3430.
# 3          3  17125301  15412770900      15413.
# 4          4  44866591  40379931900      40380.
# 5          5  13995717  12596145300      12596.
# 6          6    872418    785176200        785.
# 7          7  13966450  12569805000      12570.
# 8          8    604000    543600000        544.
# 9         NA 324318431 291886587900     291887.
4
  • @aldo_tapia.I tried working same as this but its displaying : Must group by variables found in .data. ✖ Column RasterLayer is not found. Why is this happening?
    – Walker
    Commented Jun 13, 2022 at 3:31
  • For what function? did you use as.data.frame(LULC)?
    – aldo_tapia
    Commented Jun 13, 2022 at 3:42
  • how can we know which class name is representing different class name from your calculations ?
    – Walker
    Commented Jun 23, 2022 at 12:59
  • 1
    Check raster's metadata or the other answer.
    – aldo_tapia
    Commented Jun 23, 2022 at 18:24

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.