In order to assess the vegetation change, first the classified vegetation rasters of both years need to be available. These were previously created in the classifying vegetation tutorial. You can import these two rasters, and then proceed to calculate the changes for each pixel.
The first part of the function imports the classified vegetation rasters as tiff files from the working directory. It is important that you define year_1 as the first year, and year_2 as the second year of the same area, since this allows you to see the direction of the change. The last two lines define the projection as Rijksdriehoekstelsel, since R does not know in which projection the image was calculated. Adding the corresponding projection to the raster is important to prevent distortions, and allow it to be mapped in the correct location on the basemap.
# Import rasters of years to compare as .tiff files from the working directory
year_1 = raster("Tutorial_files/veg_klompenwaard_year1.tif", convert = TRUE)
year_2 = raster("Tutorial_files/veg_klompenwaard_year2.tif", convert = TRUE)
# Define projection of vegetation raster
year_1@crs = CRS("+init=epsg:28992")
year_2@crs = CRS("+init=epsg:28992")The rasters of both years should look something like the following two images, though they are generally not shown while running the script.
Once the vegetation classification is available for both years, they need to be compared. This is done by building a new raster, which has a unique value showing which vegetation class the area had in each pixel during year 1 and during year 2. This is build by multiplying the classes from year 1 with 10, and then adding this value to year 2. This means that for example, a pixel that was classified as shrub (which has a value of 2) in year 1, and grew into trees (which has a value of 3) in the time between the two years, will have a new value of 23. This should not be read as twenty-two, but as two and three, since the first number shows the class during year one, and the second number that of year two.
The result of this is an entirely new raster, with a lot more classes. The parts which do not have values for both years are filtered out of it, and then it is saved as a tiff file again, before we proceed to the next part.
# Multiply year 1 by 10
fun = function(x) {x*10}
sample_year_1 = calc(year_1, fun)
# Create and save the vegetation change map as a .tiff in working directory
change = mosaic(sample_year_1, year_2, fun=sum)
n = c(11, 11, 12, 12, 13, 13, 14, 14, 15, 15,
21, 21, 22, 22, 23, 23, 24, 24, 25, 25,
31, 31, 32, 32, 33, 33, 34, 34, 35, 35,
41, 41, 42, 42, 43, 43, 44, 44, 45, 45,
51, 51, 52, 52, 53, 53, 54, 54, 55, 55,
10, 0, 20, 0, 30, 0, 40, 0, 50, 0)
rclmat = matrix(n, ncol=2, byrow=TRUE)
change = reclassify(change, rclmat)
writeRaster(change, filename = "klompenwaard_change.tif", format="GTiff", overwrite = TRUE)From the previously imported vegetation rasters as shown above, we will create a new change map as shown below. Note the large amount of classes, which are combinations of the vegetation classes which we have imported earlier.
In order to quickly assess locations in which the vegetation has changed at all, no matter the direction, it is possible to map those parts of your study area which have and have not changed. These we call the stable and the unstable areas.
This is done by simply reclassifying the previously made vegetation map into two new classes. One of these classes shows the stable parts, made from those vegetation change classes which have two the same values, and the other one of all changing classes. This is then saved into a new tiff file.
# Build reclassification matrix
m = c(11, 1, 12, 2, 13, 2, 14, 2, 15, 2,
21, 2, 22, 1, 23, 2, 24, 2, 25, 2,
31, 2, 32, 2, 33, 1, 34, 2, 35, 2,
41, 2, 42, 2, 43, 2, 44, 1, 45, 2,
51, 2, 52, 2, 53, 2, 54, 2, 55, 1)
rclmat = matrix(m, ncol=2, byrow=TRUE)
stable = reclassify(change, rclmat)
writeRaster(stable, filename = "klompenwaard_stable.tif", format="GTiff", overwrite = TRUE)The resulting map is shown in the interface below. It shows the stable areas in green, and the parts which have changed in red.