README.md

Land Use Change Calculus (lucCalculus)

This package implements the LUC Calculus for reasoning about land use change trajectories. Based on a set of classified time series, we build expressions to answer specific questions, such as Which events of "Forest" areas were replaced by "Pasture"?

With package "lucCalculus" is possible to build questions using Allen's interval temporal logic relationships and also others extended from their study. I suggest the reader read (Allen 1983) and (Allen 1984) for more details. Besides, is possible to generate graphics with event information and plot maps with results. Using these events the user can to perform analysis on time series data to discover important land use changes.

Prerequisites:

How to use the package:

Example 1

 
library(lucCalculus)
options(digits = 12)
#-----------------------
# 0. Open images and create a RasterBrick 
#-----------------------

# create a RasterBrick from individual raster GeoTIFF classified previously
lucC_create_RasterBrick(path_open_GeoTIFFs = c(system.file("extdata/raster/rasterItanhanga", package = "lucCalculus")),
                        path_save_RasterBrick = getwd())

# ------------- define variables to use in metadata -------------
# open files
file <- paste0(getwd(),"/rasterItanhanga.tif")
file

# create timeline with classified data from SVM method
timeline <- c("2001-09-01", "2002-09-01", "2003-09-01", "2004-09-01", "2005-09-01", "2006-09-01", 
              "2007-09-01", "2008-09-01", "2009-09-01", "2010-09-01", "2011-09-01", "2012-09-01", 
              "2013-09-01", "2014-09-01", "2015-09-01", "2016-09-01")
timeline

# new variable with raster object
rb_class <- raster::brick(file)

# ------------- define variables to plot raster -------------
# original label - see QML file, same order
label <- c("Degradation", "Double_cropping", "Single_cropping", "Forest", "Pasture", "Pasture", 
            "Pasture", "Double_cropping", "Double_cropping", "Double_cropping", "Double_cropping", 
            "Double_cropping", "Single_cropping", "Single_cropping", "Water", "Water")
label

# original colors set - see QML file, same order
colors_1 <- c("#BEEE53", "#cd6155", "#e6b0aa", "#228b22", "#7ecfa4", "#afe3c8",  "#64b376", 
              "#e1cdb6", "#b6a896", "#b69872", "#b68549", "#9c6f38", "#e5c6a0", "#e5a352", 
              "#0000ff", "#3a3aff")
colors_1

# plot rasterBrick
lucC_plot_raster(raster_obj = rb_class,
                 timeline = timeline, label = label,
                 custom_palette = TRUE, RGB_color = colors_1, plot_ncol = 4)

Fig. 1. Plot images classified from a RasterBrick

Fig. 1. Plot images classified from a RasterBrick

# Forest holds from 2001 to 2007
a <- lucC_pred_holds(raster_obj = rb_class, raster_class = "Degradation",
                     time_interval = c("2001-09-01","2007-09-01"),
                     relation_interval = "contains", label = label, timeline = timeline)
head(a)

# Cerrado holds from 2009 to 2013
b <- lucC_pred_holds(raster_obj = rb_class, raster_class = "Pasture",
                     time_interval = c("2009-09-01","2013-09-01"),
                     relation_interval = "equals", label = label, timeline = timeline)
head(b)

# only locations where b occurs after a - Before Allen's Relation
c <- lucC_relation_before(a, b)

# plot results
# individual locations over time
lucC_plot_sequence_events(c, custom_palette = FALSE, show_y_index = FALSE)

# barplot with result quantified
lucC_plot_bar_events(c, custom_palette = FALSE, pixel_resolution = 232, side_by_side = TRUE)

# plot RasterBrick with individual states over time
lucC_plot_raster_result(raster_obj = rb_class, data_mtx = c, timeline = timeline, 
                        label = label, custom_palette = TRUE, RGB_color = colors_1, 
                        relabel = FALSE, shape_point = ".", plot_ncol = 4) 
Fig. 2.(a) Locations over time Fig. 2.(a) Locations over time Fig. 2.(b) Barplot with number of states in km² Fig. 2.(b) Barplot with number of states in km²

Fig. 3. Plot images classified from a RasterBrick and states from Before relation

Fig. 3. Plot images classified from a RasterBrick and states from Before relation

Example 2

#----------------------------
# 1. RECUR predicate indicates a class that appear again
#----------------------------

forest_recur <- lucC_pred_recur(raster_obj = rb_class, raster_class = "Forest",
                                time_interval1 = c("2001-09-01","2001-09-01"),
                                time_interval2 = c("2002-09-01","2016-09-01"),
                                label = label, timeline = timeline, 
                                remove_column = FALSE)
head(forest_recur)

#-------------------
# plot some results from RECUR
# lucC_plot_sequence_events(forest_recur, custom_palette = FALSE, show_y_index = FALSE)
# lucC_plot_bar_events(forest_recur, custom_palette = FALSE)

lucC_plot_raster_result(raster_obj = rb_class, data_mtx = forest_recur,
                        timeline = timeline, label = label, custom_palette = TRUE,
                        RGB_color = colors_1, relabel = FALSE, shape_point = ".", plot_ncol = 4)
#-------------------

- Plotted RasterBrick rb_class results

Fig. 4. Plot images classified from a RasterBrick and states from RECUR relation

Fig. 4. Plot images classified from a RasterBrick and states from RECUR relation

#----------------------------
# 2. EVOLVE to verify Forest class that occurs after a different class in 2001
#----------------------------
forest_evolve <- NULL

# classes without Forest based on original label
classes <- c("Degradation", "Double_cropping", "Single_cropping", "Pasture", "Pasture", 
              "Pasture", "Double_cropping", "Double_cropping", "Double_cropping", "Double_cropping", 
              "Double_cropping", "Single_cropping", "Single_cropping", "Water", "Water")

# percor all classes
for(i in seq_along(classes)){
    print(classes[i])
    temp <- lucC_pred_evolve(raster_obj = rb_class, raster_class1 = classes[i],
                             time_interval1 = c("2001-09-01","2001-09-01"), relation_interval1 = "equals",
                             raster_class2 = "Forest",
                             time_interval2 = c("2002-09-01","2016-09-01"), relation_interval2 = "contains",
                             label = label, timeline = timeline, remove_column = FALSE)
    forest_evolve <- lucC_merge(forest_evolve, temp)
  }


#-------------------
# plot some results from EVOLVE
# lucC_plot_sequence_events(forest_evolve, custom_palette = FALSE, show_y_index = FALSE)
# lucC_plot_bar_events(forest_evolve, custom_palette = FALSE, legend_text = "Legend:")

lucC_plot_raster_result(raster_obj = rb_class, data_mtx = forest_evolve,
                        timeline = timeline, label = label, custom_palette = TRUE,
                        RGB_color = colors_1, relabel = FALSE, shape_point = ".", plot_ncol = 4)
#-------------------
Fig. 5. Plot images classified from a RasterBrick and states from EVOLVE relation

Fig. 5. Plot images classified from a RasterBrick and states from EVOLVE relation

#----------------------------
# 3. Merge both forest_recur and forest_evolve data sets
#----------------------------

forest_secondary <- lucC_merge(forest_evolve, forest_recur)
head(forest_secondary)

# plot
lucC_plot_bar_events(forest_secondary, custom_palette = FALSE, pixel_resolution = 232, legend_text = "Legend:")

# 4. Remove column 2001 because it' is not used to replace pixels's only support column
forest_secondary <- lucC_remove_columns(data_mtx = forest_secondary, name_columns = c("2001-09-01"))
head(forest_secondary)

# plot
lucC_plot_bar_events(forest_secondary, custom_palette = FALSE, pixel_resolution = 232, legend_text = "Legend:")

# 5. Plot secondary vegetation over raster without column 2001 because it' is not used to replace pixels's only support column
lucC_plot_raster_result(raster_obj = rb_class, data_mtx = forest_secondary, timeline = timeline, 
                        label = label, custom_palette = TRUE, RGB_color = colors_1, 
                        relabel = FALSE, shape_point = ".", plot_ncol = 4)


Fig. 6. Plot images classified from a RasterBrick and states from EVOLVE and RECUR relations, and also secondary vegetation locations

Fig. 6. Plot images classified from a RasterBrick and states from EVOLVE and RECUR relations, and also secondary vegetation locations

Example 3

#----------------------------
# 4. Update original raster to add new pixel value
#----------------------------
n_label <- length(label) + 1

# 1. update original RasterBrick with new class
rb_class_new <- lucC_raster_update(raster_obj = rb_class,
                                   data_mtx = forest_secondary,   # without 2001
                                   timeline = timeline,
                                   class_to_replace = "Forest",   # only class Forest
                                   new_pixel_value = n_label)     # new pixel value

head(rb_class_new)

lucC_plot_bar_events(data_mtx = rb_class_new, pixel_resolution = 232, custom_palette = FALSE)

# 2. save the update matrix as GeoTIFF images
lucC_save_GeoTIFF(raster_obj = rb_class,
                  data_mtx = rb_class_new,
                  path_raster_folder = paste0(getwd(),"/rasterItanhangaSecVeg"), 
                  as_RasterBrick = FALSE)

#------------
# create a RasterBrick from individual raster GeoTIFF, case saved as separate layers in lucC_save_GeoTIFF, as_RasterBrick = FALSE
lucC_create_RasterBrick(path_open_GeoTIFFs = paste0(getwd(),"/rasterItanhangaSecVeg"),
                        path_save_RasterBrick = getwd())

# open file RasterBrick
file <- c(paste0(getwd(),"/rasterItanhangaSecVeg.tif"))

# create timeline with classified data from SVM method
timeline <- c("2001-09-01", "2002-09-01", "2003-09-01", "2004-09-01", "2005-09-01", "2006-09-01", 
              "2007-09-01", "2008-09-01", "2009-09-01", "2010-09-01", "2011-09-01", "2012-09-01", 
              "2013-09-01", "2014-09-01", "2015-09-01", "2016-09-01")

# new variable with raster object
rb_class2 <- raster::brick(file)

# ------------- define variables to plot raster -------------
# original label - see QML file, same order and new class secondary vegetation
label2 <- c("Degradation", "Double_cropping", "Single_cropping", "Forest", "Pasture", "Pasture", 
            "Pasture", "Double_cropping", "Double_cropping", "Double_cropping", "Double_cropping", 
            "Double_cropping", "Single_cropping", "Single_cropping", "Water", "Water", 
            "Secondary_vegetation")

# original colors set - see QML file, same order
colors_2 <- c("#BEEE53" , "#cd6155", "#e6b0aa", "#228b22", "#7ecfa4", "#1e174d", "#afe3c8", 
              "#64b376", "#e1cdb6", "#b6a896", "#b69872", "#b68549", "#9c6f38", "#e5c6a0", 
              "#e5a352", "#0000ff", "#3a3aff")

# plot rasterBrick
lucC_plot_raster(raster_obj = rb_class2,
                 timeline = timeline, label = label2,
                 custom_palette = TRUE, RGB_color = colors_2, plot_ncol = 4)


- Plotted RasterBrick rb_class results with new class Secondary Vegetation

Fig. 7. Plot a RasterBrick and states with new secondary vegetation class

Fig. 7. Plot a RasterBrick and states with new secondary vegetation class

Example 4

#----------------------------
# 5. Discover Forest and Secondary vegetation - LUC Calculus
#----------------------------

secondary.mtx <- lucC_pred_holds(raster_obj = rb_class2, raster_class = "Secondary_vegetation",
                                 time_interval = c("2001-09-01","2016-09-01"),
                                 relation_interval = "contains", label = label2, timeline = timeline)
head(secondary.mtx)

forest.mtx <- lucC_pred_holds(raster_obj = rb_class2, raster_class = "Forest",
                              time_interval = c("2001-09-01","2016-09-01"),
                              relation_interval = "contains", label = label2, timeline = timeline)
head(forest.mtx)

Forest_secondary.mtx <- lucC_merge(secondary.mtx, forest.mtx)
head(Forest_secondary.mtx)

# plot results
lucC_plot_bar_events(data_mtx = Forest_secondary.mtx, custom_palette = TRUE, 
                     RGB_color = c("black", "gray60"), #c("#228b22", "#7ecfa4"),
                     pixel_resolution = 231.656, side_by_side = TRUE,
                     relabel = TRUE, original_labels = c("Forest", "Secondary_vegetation"),
                     new_labels = c("Forest", "Secondary vegetation"))

# Compute values
measuresFor_Sec <- lucC_result_measures(data_mtx = Forest_secondary.mtx, pixel_resolution = 231.656)
measuresFor_Sec

#------------
   Years              Classes Pixel_number        Area_km2   Cumulative_Sum Relative_Frequency Cumulative_Relative_Frequency
1   2001               Forest        46647 2503.2880404674  2503.2880404674    10.463261831648               10.463261831648
2   2002               Forest        43152 2315.7306048031  4819.0186452705     9.679307877447               20.142569709096
3   2003               Forest        38585 2070.6448226346  6889.6634679050     8.654896515835               28.797466224931
4   2004               Forest        33751 1811.2306183423  8700.8940862474     7.570595109653               36.368061334583
5   2005               Forest        29849 1601.8317302273 10302.7258164746     6.695348091257               43.063409425841
6   2006               Forest        27909 1497.7225956954 11800.4484121700     6.260191962173               49.323601388013
7   2007               Forest        27226 1461.0697405999 13261.5181527700     6.106990087861               55.430591475875
8   2008               Forest        25062 1344.9397575448 14606.4579103148     5.621589127377               61.052180603252
9   2009               Forest        24618 1321.1127185076 15927.5706288225     5.521996693711               66.574177296963
10  2010               Forest        24134 1295.1390993770 17222.7097281995     5.413431968723               71.987609265685
11  2011               Forest        22913 1229.6147420248 18452.3244702243     5.139552776139               77.127162041824
12  2012               Forest        21345 1145.4688023619 19597.7932725862     4.787838956343               81.915000998167
13  2013               Forest        20781 1115.2020230444 20712.9952956306     4.661329648712               86.576330646880
14  2014               Forest        20236 1085.9548692713 21798.9501649019     4.539082179459               91.115412826339
15  2015               Forest        19879 1066.7966419373 22865.7468068392     4.459004479416               95.574417305756
16  2016               Forest        19730 1058.8006310893 23924.5474379285     4.425582694244              100.000000000000
17  2001 Secondary_vegetation            0    0.0000000000     0.0000000000     0.000000000000                0.000000000000
18  2002 Secondary_vegetation          624   33.4866494577    33.4866494577     0.870256474624                0.870256474624
19  2003 Secondary_vegetation         1244   66.7586409060   100.2452903636     1.734934382104                2.605190856728
20  2004 Secondary_vegetation         2781  149.2409809964   249.4862713601     3.878498807581                6.483689664310
21  2005 Secondary_vegetation         2858  153.3731476763   402.8594190364     3.985886225123               10.469575889433
22  2006 Secondary_vegetation         4259  228.5571154490   631.4165344854     5.939779367669               16.409355257102
23  2007 Secondary_vegetation         6334  339.9109577962   971.3274922816     8.833661074153               25.243016331255
24  2008 Secondary_vegetation         3902  209.3988881151  1180.7263803967     5.441892249976               30.684908581231
25  2009 Secondary_vegetation         7631  409.5138173260  1590.2401977227    10.642511470929               41.327420052160
26  2010 Secondary_vegetation         5969  320.3234144436  1910.5636121663     8.324616822169               49.652036874329
27  2011 Secondary_vegetation         4785  256.7846436778  2167.3482558440     6.673360947241               56.325397821570
28  2012 Secondary_vegetation         5551  297.8916524671  2465.2399083112     7.741656555514               64.067054377083
29  2013 Secondary_vegetation         6753  362.3963842750  2827.6362925862     9.418015982595               73.485070359678
30  2014 Secondary_vegetation         6068  325.6362001748  3153.2724927610     8.462686359009               81.947756718687
31  2015 Secondary_vegetation         6417  344.3651114901  3497.6376042511     8.949416342412               90.897173061099
32  2016 Secondary_vegetation         6527  350.2682067471  3847.9058109982     9.102826938901              100.000000000000

- Quantity of Forest and Secondary vegetation after LUC Calculus application

Fig. 8. Number of Forest and Secondary Vegetation

Fig. 8. Number of Forest and Secondary Vegetation

Example 5

#----------------------------
# 5. Discover direct land use transitions
#----------------------------

class1 <- c("Forest")
classes <- c("Degradation", "Double_cropping", "Pasture", "Single_cropping")

direct_transi.df <- NULL

# along of all classes
for(x in 2:length(timeline)){
  t_1 <- timeline[x-1]
  t_2 <- timeline[x]
  cat(paste0(t_1, ", ", t_2, sep = ""), "\n")

  # moves across all classes
  for(i in seq_along(classes)){
    cat(classes[i], collapse = " ", "\n")
    temp <- lucC_pred_convert(raster_obj = rb_class2, raster_class1 = class1,
                              time_interval1 = c(t_1,t_1), relation_interval1 = "equals",
                              raster_class2 = classes[i],
                              time_interval2 = c(t_2,t_2), relation_interval2 = "equals",
                              label = label2, timeline = timeline)
    direct_transi.df <- lucC_merge(direct_transi.df, temp)
  }
  cat("\n")
}

Forest_others <- direct_transi.df
head(Forest_others)
str(Forest_others)

# plot results
lucC_plot_frequency_events(data_mtx = Forest_others,
                           pixel_resolution = 232, custom_palette = FALSE)

lucC_plot_bar_events(Forest_others, custom_palette = FALSE, pixel_resolution = 232, 
                      legend_text = "Legend:", side_by_side = FALSE)

# replace classes by new legend
Forest_others[c(3:ncol(Forest_others))] <- 
            ifelse(Forest_others[c(3:ncol(Forest_others))] == "Degradation", "Forest_Degradation",
            ifelse(Forest_others[c(3:ncol(Forest_others))] == "Double_cropping", "Forest_Double_cropping",
            ifelse(Forest_others[c(3:ncol(Forest_others))] == "Pasture", "Forest_Pasture",
            ifelse(Forest_others[c(3:ncol(Forest_others))] == "Single_cropping", "Forest_Single_cropping"
            , ""))))

lucC_plot_bar_events(Forest_others, pixel_resolution = 231.6465, custom_palette = TRUE, 
                      RGB_color = c("#6e8b3d", "#a2cd5a", "#006400", "#b4eeb4"), 
                      legend_text = "Land use transitions:", relabel = TRUE, 
                      original_labels = c("Forest_Degradation", "Forest_Double_cropping", "Forest_Pasture", "Forest_Single_cropping"), 
                      new_labels = c("Degradation", "Forest to Double Cropping", "Forest to Pasture", "Forest to Single Cropping"), 
                      side_by_side = FALSE)

- Direct land use transition from Forest class to others in Itanhangá municipality

Fig. 9 Quantity of direct land use transitions

Fig. 9. Quantity of direct land use transitions



ammaciel/lucCalculus documentation built on June 13, 2020, 4:57 a.m.