Developing a method to characterize land use history using Landsat [PDF]

Laboratory of Geo-Information Science and Remote Sensing Thesis Report GIRS-2015-19

Developing a method to characterize land use history using Landsat time series as an idle

land

early

detection

method

in

Indonesia

22nd May 2015

Marthalina Indhawati

i

Developing a method to characterize land use history using Landsat time series as an idle land early detection method in Indonesia

Marthalina Indhawati

Supervisors: dr. ir. Lammert Kooistra dr. ir. Jan Verbesselt

A thesis submitted in partial fulfillment of the degree of Master of Science at Wageningen University and Research Centre, The Netherlands.

22nd May 2015 Wageningen, The Netherlands

Thesis code number: GRS-80436 Thesis report: GIRS-2015-19 Wageningen University and Research Centre Laboratory of Geo-Information Science and Remote Sensing ii

Acknowledgements “He has made everything beautiful in its time.” (Ecclesiastes 3:11)

First of all, I would like to thank to my Lord, Jesus Christ, for giving me strength and grace for every single second in my life. My deepest thanks to dr. ir. Lammert Kooistra and dr. ir. Jan Verbesselt who always gave me an extensive supervision, discussion, and advices during my thesis. It helped me a lot to become critical in scientific direction. I always believe that a teacher is not only the one who transfers the knowledge, but also who teaches the value of life and keep supporting their students. Both are great supervisors and teachers for me. I am also grateful to Prof.dr. Martin Herold and dr.ir. Jan Clevers for giving me valuable feedback during my thesis. Special thanks to dr. Harm Bartholomeus and Benjamin DeVries who support me by answering all my questions related to technical and substantial matters. For Nikos Tziolas, Tom Bewernick, and Latifa, thanks very much for all the technical discussions we had. Next, I also want to say thanks for my colleagues from BPN, LAPAN, and BIG for all the data and the cooperation when I was in Indonesia and in Wageningen. Many thanks to all my friends in Indonesia and Wageningen, especially to Marnix, Sukmo, and Ratih, who always support me during this period, we laugh and struggle in togetherness. Lastly, I would like to say thank you to my beloved parents, brother, and sister who always keep praying for my health and success in my study. I hope this thesis would give a contribution in Indonesia land management in which land should give prosperity to all Indonesia citizens.

“We abuse land because we see it as a commodity belonging to us. When we see land as a community to which we belong, we may begin to use it with love and respect” (Aldo Leopold)

iii

Abstract Nowadays, idle land (land abandonment) is a serious problem in Indonesian land management. About 7.3 million hectare is estimated as idle land spread over the whole country. It causes many disadvantages and increases land conflicts in Indonesia. Hence, the national government gave instructions to the National Land Agency of the Republic of Indonesia to control idle land in Indonesia. However, idle land controlling has some major obstacles during indicated idle land inventory; one of them is that Hak Guna Usaha land parcels (HGU - right to cultivate) are located far away from the capital city, resulting in difficulties to enforce controlling mechanisms. Remote sensing has been used successfully for many years in detection of changes and together with time series analysis it will lead to a better understanding of land use changes. In this study we provide a possible method for performing indicated idle land inventory based on remote sensing methods in which, we aimed: 1) to find out whether there is any land clearance/deforestation on HGU land parcels or not and 2) to find what is the HGU land cover, four years after the publication of idle land controlling regulation and compare, whether it is the same as the purpose of the entitlement. Firstly, we used Landsat satellite images for the years 2008 until 2014 and applied the Break For Additive Season and Trend Monitor (BFAST Monitor) with a monitoring period from the 1st of January 2013 to the 1st of July 2014 to detect spatial and temporal changes for five HGU land parcels in the Riau Province – Indonesia. We combined the breakpoints and set a threshold of 0, -0.025, -0.05, -0.075, -0.1 for the magnitude value to evaluate which threshold would obtain the best accuracy in detecting land clearance for the five study areas. We validated our BFAST monitor results, by using SPOT 6 images of the years 2013 and 2014. The purpose of this validation was to validate, the accuracy of the BFAST monitor in detecting land clearance. We chose the cross point between user’s and producer’s accuracy as an unbiased estimator in order to get the best threshold for detecting land clearance for each study area. In general, a magnitude threshold of -0.075 gave the best results in land clearance detection for this particular study area. This threshold value corresponded to the highest overall accuracy of 88.24%. We concluded that the BFAST Monitor could be used for detecting land clearance in the five HGU land parcels by combining breakpoints and magnitude values. Secondly, we performed a maximum likelihood supervised classification for a Landsat 8 image using a SPOT 6 image acquired on 17 June 2014 as the ground truth for observing the existing land cover. We validated our classification results by using a SPOT 6 image for the year 2014 in combination with the Kampar land cover map of the year 2014. We performed random sampling and combined the results into one confusion matrix for the five study areas. The results of this accuracy assessment showed an overall accuracy of 86.13% and a Kappa coefficient of 0.78, which corresponds to a substantial agreement. We concluded that detecting bare soil is more accurate than vegetation i.e. oil palm, forest, and grassland. In general, it can be concluded that Landsat time series are suitable to be used for idle land early detection. However, high cloud cover was a major aspect causing a low data availability, which limited the usability of our methodology. By combining all the results in this study, it gives many opportunities for the National Land Agency of the Republic of Indonesia to have an idle land early detection method to decide whether a land parcel will be targeted in the idle land controlling or not. A continuously development in the idle land controlling stages is necessary for the National Land Agency of the Republic of Indonesia to gain a proper land management in Indonesia. Keywords: idle land, Landsat, time series, BFAST Monitor, land cover classification, Indonesia. iv

Table of Contents Acknowledgements ............................................................................................................................... iii Abstract ................................................................................................................................................. iv List of tables......................................................................................................................................... viii List of figures ......................................................................................................................................... ix List of acronyms, special terms and abbreviations ................................................................................. x 1. Introduction......................................................................................................................................... 1 1.1

Background information........................................................................................................... 1

1.2

Problem statement .................................................................................................................. 2

1.3

Research objectives and research questions ............................................................................. 4

1.4

Overview of report ................................................................................................................... 5

2. Theoretical background ....................................................................................................................... 6 2.1

The National Land Agency of the Republic of Indonesia ............................................................ 6

2.2

Idle land controlling in Indonesia.............................................................................................. 8

2.3

Oil palm plantation in Indonesia ............................................................................................... 9

2.4

Disturbance detection in Landsat time series using BFAST Monitor ........................................ 11

2.5

Land cover classification ......................................................................................................... 13

2.6

Accuracy assessment.............................................................................................................. 14

3. Materials and Methods ..................................................................................................................... 16 3.1

Study area .............................................................................................................................. 16

3.2

Materials................................................................................................................................ 17

3.2.1

HGU land parcels and Kampar regency land cover map .................................................. 17

3.2.2

Satellite images datasets ................................................................................................ 17

3.2.3

Software......................................................................................................................... 21

3.3

Data analysis methods ........................................................................................................... 21

3.3.1

Preprocessing ................................................................................................................. 23

3.3.2

The determination of forest mask................................................................................... 23 v

3.3.3

Land clearance detection in Landsat time series using BFAST Monitor ............................ 23

3.3.4

Land clearance validation and accuracy assessment ....................................................... 25

3.3.5

Land cover classification and analysis ............................................................................. 25

3.3.6

Land cover classification accuracy assessment ................................................................ 28

4. Results ............................................................................................................................................... 30 4.1

Observing land clearance in HGU land parcels ........................................................................ 30

4.1.1

Landsat multi-temporal NDVI profile .............................................................................. 30

4.1.2

Land clearance detection using BFAST monitor ............................................................... 32

4.1.3

Land clearance validation and accuracy assessment ....................................................... 35

4.2

Observing land cover classification ......................................................................................... 38

4.2.1

Land cover based on Landsat maximum likelihood supervised classification ................... 38

4.2.2

Land cover classification validation and accuracy assessment ......................................... 41

5. Discussions ........................................................................................................................................ 44 5.1

Landsat time series for idle land early detection..................................................................... 44

5.1.1

Observing NDVI time series for land use history.............................................................. 44

5.1.2

BFAST monitor in detecting land clearance ..................................................................... 45

5.1.3

Potential errors and limitations from BFAST monitor analysis ......................................... 45

5.2

Landsat images classification to support for idle land controlling............................................ 48

5.3

The combination of NDVI time series profile, BFAST monitor, and Landsat classification as idle land early detection ........................................................................................................ 50

5.4

Future research needs............................................................................................................ 52

6. Conclusions........................................................................................................................................ 54 References ............................................................................................................................................ 55 Appendixes ........................................................................................................................................... 62 Appendix I. The annual median NDVI of HGU A, B, C, and E for period 2008 until 2014. ...................... 62 Appendix II. The forest mask of HGU A, B, C, D, and E. ........................................................................ 66 Appendix III. The BFAST monitor output of HGU A, B, C, D, and E. ...................................................... 67 Appendix IV. The accuracy assessment of BFAST monitor in detecting land clearance. ....................... 69 vi

Appendix V. The BFAST monitor plots for single pixel time series. ...................................................... 77

vii

List of tables Table 1. Favorable soil conditions for oil palm plantations ........................................................ 10 Table 2. General overview of HGU land parcels........................................................................ 17 Table 3. The overview of Landsat TM, ETM+, Landsat 8 spectral bands ................................... 18 Table 4. Monthly frequency of Landsat data from 2008 to 2014 for path 127 row 60 ................ 19 Table 5. The overview of SPOT 6 spectral bands ...................................................................... 20 Table 6. The descriptions of land cover classes for Landsat classification.................................. 26 Table 7. The basic characteristics of image interpretation .......................................................... 26 Table 8. The determination of reference classes based on SPOT 6 images visual assessment ..... 27 Table 9. The mean and standard deviation clear sky pixel observation of the forest NDVI time series stack for history and monitoring period for five HGU land parcels ............. 32 Table 10. The number of pixels for land clearance class for five different magnitude thresholds... ............................................................................................................... 33 Table 11. The BFAST monitor accuracy assessment in detecting land clearance for HGU B ..... 36 Table 12. The number of pixels in each class type as reference data for training and validating the Landsat classification for the five study areas ....................................................... 39 Table 13. The number of pixels and the area (ha) after performing a land cover classification for the five selected study areas ................................................................................. 39 Table 14. The accuracy assessment of classification using SPOT 6 images for five HGU land parcels ................................................................................................................ 42 Table 15. The BFAST monitor accuracy assessment in detecting land clearance for HGU A..... 70 Table 16. The BFAST monitor accuracy assessment in detecting land clearance for HGU C ..... 72 Table 17. The BFAST monitor accuracy assessment in detecting land clearance for HGU D..... 74 Table 18. The BFAST monitor accuracy assessment in detecting land clearance for HGU E ..... 76

viii

List of figures Figure 1. The organizational structure of NLA. ............................................................................... 6 Figure 2. The management of high resolution images in Indonesia and the National Land Agency of the Republic of Indonesia ................................................................................ 7 Figure 3. The idle land controlling main steps. ................................................................................ 8 Figure 4. Oil palm plantation in Riau Province .............................................................................. 10 Figure 5. BFAST monitor plot output for single pixel time series .................................................. 12 Figure 6. The principle of maximum likelihood classifier .............................................................. 13 Figure 7. Study area of five HGU land parcels in Riau Province – Indonesia ................................. 16 Figure 8. The number of Landsat TM/ETM+ images from 2008 to 2014 for path 127 row 60........ 19 Figure 9. Oil palm plantation from SPOT 6 satellite images .......................................................... 20 Figure 10. The steps to detect land use history as part of idle land early detection by using the BFAST Monitor ..................................................................................................... 22 Figure 11. Demonstration of the sequential monitoring approach used by DeVries et al., 2015 ...... 24 Figure 12. Annual median NDVI of HGU D for period 2008 until 2014. ....................................... 31 Figure 13. Land clearance detected by BFAST monitor for HGU B .............................................. 34 Figure 14. Land clearance detected by BFAST monitor for single pixel time series in HGU B ...... 35 Figure 15. Overall accuracy (OA), producer’s accuracy (PA), and user’s accuracy (UA) of land clearance detection as a function of magnitude for HGU A, B, C, D, and E .......... 37 Figure 16. The comparison of Landsat before (left) and after (right) applying CFmask.................. 38 Figure 17. Land covers classification of HGU A, B, C, D, and E from Landsat 8 image of 12 June 2014. ................................................................................................................... 40 Figure 18. Overlay between HGU A, B, C, D, E and Kampar land cover map of 2014 with a scale 1:50,000 ........................................................................................................... 43 Figure 19. Example of commission errors arising from low data availability in the temporal history period ................................................................................................ 47 Figure 20. Spectral signature plot of land cover classes in Landsat 8 classification ........................ 49 Figure 21. The combination of BFAST monitor and Landsat classification as idle land early detection ...................................................................................................................... 51 Figure 22. The annual median NDVI of HGU A for period 2008 until 2014 .................................. 62 Figure 23. The annual median NDVI of HGU B for period 2008 until 2014 .................................. 63 Figure 24. The annual median NDVI of HGU C for period 2008 until 2014 .................................. 64 Figure 25. The annual median NDVI of HGU E for period 2008 until 2014 .................................. 65 Figure 26. The determination of HGU A’s forest mask .................................................................. 66 Figure 27. The forest masks of HGU B, C, D, and E ..................................................................... 66 Figure 28. The output of BFAST monitor for HGU A ................................................................... 67 Figure 29. The output of BFAST monitor for HGU B ................................................................... 67 Figure 30. The output of BFAST monitor for HGU C ................................................................... 67 Figure 31. The output of BFAST monitor for HGU D ................................................................... 68 Figure 32. The output of BFAST monitor for HGU E.................................................................... 68 Figure 33. Land clearance detected by BFAST monitor for HGU A .............................................. 69 Figure 34. Land clearance detected by BFAST monitor for HGU C .............................................. 71 Figure 35. Land clearance detected by BFAST monitor for HGU D .............................................. 73 Figure 36. Land clearance detected by BFAST monitor for HGU E............................................... 75 Figure 37. Some examples of BFAST monitor plot for single pixel time series in detecting disturbances/land clearances .......................................................................................... 77

ix

List of acronyms, special terms and abbreviations AOI BAL BFAST BFAST Monitor BIG BPN RI CCDC EROS ESPA ETM+ GPS HGU IR LAPAN LCCS-UNFAO LEDAPS MODIS MODIS VCF NDVI NIR NLA NPP OA OLI PA RMSE SLC SPOT SWIR TIRS TM UA USGS UTM WGS

: Area of Interest : Basic Agrarian Law : Breaks for Additive Seasonal and Trend : Breaks for Additive Seasonal and Trend Monitor : Badan Informasi Geospasial : Badan Pertanahan Nasional Republik Indonesia : Continuous Change Detection and Classification : Earth Resources Observation and Science : EROS Science Processing Architecture : Enhanced Thematic Mapper plus : Global Positioning System : Hak Guna Usaha : Infra Red : Lembaga Penerbangan dan Antariksa Nasional : Land Cover Classification System United Nation – Food and Agriculture Organization : Landsat Ecosystem Disturbance Adaptive Processing System : Moderate Resolution Imaging Spectroradiometer : Moderate Resolution Imaging Spectroradiometer Vegetation Continuous Fields : Normalized Difference Vegetation Index : Near Infra Red : National Land Agency : Net Primary Productivity : Overall Accuracy : Operational Land Imager : Producer’s Accuracy : Root Mean Square Error : Scan Line Corrector : Satellite Pour l’Observation de la Terre : Short Wave Infra Red : Thermal Infrared Sensor : Thematic Mapper : User’s Accuracy : The United States Geological Survey : Universal Transverse Mercator : World Geodetic System

x

1. Introduction 1.1

Background information

Indonesia, as an agrarian country, considers land as a natural resource which is very valuable for the whole society. Land is a basis for the cultural and social development, and a source of prosperity. Therefore, the government puts high attention to land development and management in Indonesia (Winoto, 2009). Indonesia is an archipelago country with approximately 17,504 islands. It has a total area of 9.8 million km2, from which about 1.9 million km2 is land area (Badan Pusat Statistik - National Statistic Agency of the Republic of Indonesia, 2014). Administratively, since July 2013, Indonesia has 34 provinces, 412 districts, and 93 cities (Kementerian Dalam Negeri Republik Indonesia -Ministry of Home Affairs of the Republic of Indonesia, 2014a). State policies related to land have been defined in two basic regulations: Undang-undang Dasar Tahun 1945 (the Constitution 1945) and Undang-undang Nomor 5 Tahun 1960 Tentang Peraturan Dasar Pokok-pokok Agraria (the Basic Agrarian Law No. 5/1960). The Constitution 1945 article 33 sub-article 3 mentions that earth, water and natural richness inside are controlled by state and must be utilized for welfare of the people. Furthermore, as its implementation, the Basic Agrarian Law No. 5/1960 (BAL No. 5/1960) regulates the basic rules about land management in Indonesia. From the total Indonesia land area, 70% is zoned as forest land and 30% is non-forest land which is administrated by the National Land Agency (Bell and Srinivas, 2013). Badan Pertanahan Nasional Republik Indonesia (BPN) or the National Land Agency of the Republic of Indonesia (NLA) is a governmental institution, non-ministry, who is responsible to the President of the Republic of Indonesia in accordance with Peraturan Presiden Nomor 6 Tahun 2006 Tentang Badan Pertanahan Nasional (Presidential Decree Number 10/2006). In September 2013, this regulation was replaced with Peraturan Presiden Nomor 63 Tahun 2013 Tentang Badan Pertanahan Nasional Republik Indonesia (Presidential Decree Number 63/2013), and again in Januari 2015 this regulation was replaced with Peraturan Presiden Nomor 20 Tahun 2015 Tentang Badan Pertanahan Nasional Republik Indonesia (Presidential Decree Number 20/2015). Generally, the main task of NLA is supporting the President in land administration and land management at national, regional and sectoral levels, which is about land surveying and mapping, land rights and registration, land planning, land controlling, land acquisitioning for development of public interest, and land disputes assisting. Based on the Basic Agrarian Law No. 5/1960, NLA as the representation of the state is able to give land rights to Indonesian citizens, one of the land rights is Hak Guna Usaha (HGU) or right to cultivate ( BAL No.5/1960 chapter II article 16 sub article 1). HGU is the right to cultivate the land which is directly controlled by the state for 25 until 35 years for agriculture, fisheries or farms (BAL No.5/1960 chapter IV article 28 sub article 1). For HGU holders, they have to carry out agriculture, plantation, fishery and / or farms as intended and act according to the requirements as defined in the decision of rights determination (Government Regulation No. 40 Year 1996 article 12 sub article 1); otherwise, HGU is able to be taken back by NLA (BAL 1

No.5/1960 chapter IV article 34). Furthermore, land rights have social function (BAL No.5/1960 chapter I article 6) which means land should be useful not only for the owner but also for other people, society, nation and state. HGU lands come from state land or land from forest conversion which means that the mechanism should be carried out in accordance with the provisions of all related regulations. All land rights have to be utilized in three years after the publication of the land certificate; otherwise, NLA has the possibility to identify that particular land, as indicated idle land (Regulation of the Head Office of National Land Agency Number 4/2010). Indonesia has some critical issues related to land use and land rights. Firstly, in the mid-1980s, the Indonesian government has prompted the developing of non-oil and gas products, especially plantation or estate crops, to support domestic economy. This policy leads to pressure on the Indonesian forests, because of large scale forest conversion. In 2000, 2.4 million hectare is used for estate crops (Holmes, 2000; Pagiola, 2000); most of them took place in Riau, Jambi, South Sumatra, and West, Central and East Kalimantan. These activities change the trend of land use change in Indonesia (Pagiola, 2000). Secondly, inequality in the distribution of land ownership between large companies and small scale farmers, which leads to another serious obstacle: idle or abandon land which is estimated about 7.3 million hectare (Winoto, 2009). Idle land is granted land by Government that is not used or cultivated based on the purpose of entitlement (President Decree Number 11/2010). Indonesia government gives high attention in idle land problems. Some impacts of idle lands or land abandonments are economic loss from potential land; land cannot give the public welfare; access obstruction of socio economic; vulnerability of food security and national economy; environmental degradation; increase of land conflicts and disputes; delays in the provision of land for public infrastructure development (Badan Pertanahan Nasional, 2010). Therefore, in 2010 and also preceded in 2008, the President of the Republic of Indonesia, Susilo Bambang Yudhoyono gave instruction to the NLA to control idle lands in Indonesia accordance with laws and related rules. Idle land controlling becomes one of main focus areas for land policy reformation by the NLA (Winoto, 2009). Since 2010, the NLA started the idle land controlling in Indonesia based on Government Regulation Number 11/2010, Regulation of the Head Office of National Land Agency Number 4/2010, and Regulation of the Head Office of National Land Agency Number 9/2011.

1.2

Problem statement

Based on the indicated idle land inventory, the NLA targeted 1656 objects for idle land controlling in the period 2011 until 2013, which corresponded to about 3.5 million hectare indicated idle land all over Indonesia. However, until 2014, the idle land controlling still cannot reach maximum result. Until March 2014, only 104 land parcels out of the total of 1656 objects were determined as idle land, which corresponded to 68,969.85 hectare. Many obstacles appear during this process, both during the indicated idle land inventory as well as during the idle land controlling itself. Some of those obstacles occur because of lack of understanding about the existing regulations and inaccurate objects of idle land controlling. The indicated idle land inventory will be the source to determine the object of idle land controlling, which needs a good planning for both budget and human resources. Therefore, the NLA should put more attention in indicated idle land inventory. However, indicated idle land inventory has its own obstacle 2

because most of land rights, especially HGU, are located at inland areas far away from the capital city. This means that a lot of budget and human resources are necessary to perform the inventory. The new technologies have the possibility to solve these obstacles and remote sensing might be the answer to these obstacles. Remote sensing has been used successfully for many years in detection of changes (Lu et al., 2004) and time series analysis by using remote sensing techniques will give better understanding about land surface and land use changing (Herold et al., 2006; Yang and Lo, 2002). The new technologies in combination with the fast developments in the field of remote sensing also can be used to support oil palm plantation management. It is beneficial for oil palm planter to monitor characteristics of oil palm plantation for instance the phenology of oil palm (Ibrahim et al., 2001). Moreover, remote sensing is considered as a more reliable, timeliness and low cost data collection technique than sample field surveys for oil palm planning and management (Jusoff, 2009). Some studies have been performed to monitor oil palm plantations in Indonesia. Santoso et al. (2011) used QuickBird satellite images to detect the basal stem rot disease and its spatial pattern in oil palm plantations in North Sumatra. They used six vegetation indices derived from visible and near infrared bands (NIR) to identify palms infected by the disease. Sunaryathy et al. (2010) used Moderate Resolution Imaging Spectroradiometer (MODIS) to estimate Net Primary Productivity (NPP) for oil palm plantation in South Sulawesi. Another study is from Gandharum (2010), he used FORMOSAT-2 satellite images to classify growing stages of oil palms. For idle land controlling, the NLA needs high resolution satellite imagery in order to observe land use in very accurate and detailed level. However, the high price of high-resolution satellite imagery causes that the high-resolution satellite imagery purchasing cannot be implemented every year. Therefore, free satellite imagery, e.g., Landsat, and time series of satellite images are future possibilities to support land use change monitoring and idle land early detection. Landsat satellite images have some advantages, e.g. free download (Jia et al., 2014; USGS, 2008; Xian et al., 2009) and medium spatial resolution (Maxwell et al., 2010). However, using Landsat data has also some drawbacks, like: relatively low temporal frequency and cloud cover in tropical forest areas (Hansen et al., 2009; Zhu and Woodcock, 2014). Landsat has been used for many purposes, e.g. land cover detection and classification, forest disturbance monitoring, and idle land identification (De Vries et al., 2015; Maxwell et al., 2010; Zhu and Woodcock, 2014) which used (supervised) classification and time series approach to detect the changes. Regarding change detection, the classification method is considered as an old method and is less accurate when performing land cover classification. Meanwhile, the time series approach is possible to make better understanding of the event processes. Moreover, the time series approach is capable to detect changes in near real time (Verbesselt et al., 2012; Zhu and Woodcock, 2014). In order to perform idle land early detection, this study might be able to detect whether the HGU holder performs land clearance for areas larger than 20 hectare, in which most of HGU land cover will change from forest to bare soil. In general, this change in land cover will take place by a quick removal process. The detection of land clearance was followed by monitoring the land cover during a monitoring period. By using the 30 m spatial resolution of Landsat data, it might be difficult to detect the specific plants or utilization as mentioned in the purpose of entitlement; 3

it might need satellite images with a higher spatial resolution to support land cover detection. However, the presence of bare soil and grassland by the end of monitoring time indicates the HGU as an idle land. A study performed by Maxwell et al. (2010) has determined how to detect idle land represented as bare soil using the Normalized Difference Vegetation Index (NDVI) for one year Landsat images in 2000. Their research used a classification approach that intersected the Landsat data with the California Department of Water Resources land use map. Since the research only used one year Landsat images, there was no clear description about any phenomena and process that might take place in the study area. Moreover, the term idle land in their research has some differences with the term used by the NLA. The definition of idle land by the NLA emphasis on all uses that are not used or cultivated based on the purpose of entitlement, e.g. bare soil, grassland, different plantation between the entitlement and the reality. The publication year of idle land controlling regulation in 2010 was the starting point to monitor land cover of HGUs in Indonesia. For all HGU certificates which are published after 2010, three years after the publication of certificate is the key whether a HGU will be targeted in idle land controlling or not. Hence, all changes in this time are very important to be monitored especially the land clearance process and the existence of bare soil or grassland. Meanwhile, for all HGU certificates which were published before 2010, is still monitored by NLA and able to be targeted in idle land controlling. Therefore, land use history and existing land cover play an important role as a consideration for NLA to determine the idle land controlling target. One of time series methods that are possible to detect abrupt changes is Breaks For Additive Seasonal and Trend (BFAST) (Verbesselt et al., 2010a). BFAST has been used to detect phenological changes in vegetation indices, droughts, and deforestation (Verbesselt et al., 2010b; Verbesselt et al., 2012; Hutchinson et al., 2015; DeVries et al., 2015) but has not yet been validated to detect land clearance in idle land controlling. Therefore, in this thesis we evaluated the capacity of BFAST to detect sudden changes related to land clearance.

1.3

Research objectives and research questions

The overall objectives of this thesis are to test and evaluate a Landsat based change detection method to detect land clearance, in which most of the changes will be from forest-bare soil- oil palm. The specific objectives of this thesis are: - To find out whether there is any land clearance/deforestation on HGU land parcel or not. - To find what is the HGU land cover, 4 years after the publication of idle land controlling regulation and compare if it is the same as the purpose of the entitlement. This thesis aims to answer these following questions: - How sensitive is BFAST monitor in detecting land clearance? - How accurate is Landsat maximum likelihood supervised classification to distinguish land covers? - How accurate is Landsat maximum likelihood supervised classification to detect idle land? 4

1.4

Overview of report

This thesis investigates the possibilities of using Landsat time series datasets to support idle land controlling in Indonesia. We used Landsat time series and the BFAST monitor to detect land clearance, which can be considered as an indicator that the land owners already use their land as stated in the entitlement published by BPN. Furthermore we used Landsat images to perform a land cover classification and tested its ability in distinguishing land cover classes for idle land controlling. Chapter one provides the problems and issues related to idle land controlling in Indonesia. It also describes the problems due to idle land controlling and the techniques chosen. Moreover, it provides the research objectives, research questions, and an overview of report. Chapter two includes the theoretical background and literature review for this thesis. We describe the structure of the National Land Agency of the Republic of Indonesia, provide information about the current state of idle land controlling in Indonesia and the requirements of BPN and provide a theoretical background of oil palm plantation in Indonesia. Moreover we provide information about disturbance detection in Landsat time series using the BFAST monitor, about maximum likelihood supervised land cover classification, and about the accuracy assessment of the classification results. Chapter three gives an overview of the methodology used in this thesis. It briefly describes the study area, the data sets, and methodology used. It provides a flowchart of the whole process, which includes all steps, from the data obtained until the final result. For every process a description is given which includes both the data and methods used to answer all the research questions. Chapter four presents the main results obtained for answering the research questions. It presents the results of observing land clearance in HGU land parcels and the land cover classification results. Moreover, we describe in detail the Landsat multi temporal NDVI profile, land clearance detection using the BFAST monitor, and land clearance validation and accuracy assessment. Besides that, we present a land cover classification based on Landsat maximum likelihood supervised classification using SPOT 6 images as the ground truth. Furthermore, the results of a land cover classification accuracy assessment and validation are given by using the Kampar land cover map with a scale 1:50,000 and SPOT 6 images. Chapter five discusses the insights gathered from the main results in this thesis. Moreover, it also suggests possible alternatives to be investigated in future research. Finally, chapter six summarizes the main conclusions obtained with regard to the objectives.

5

2. Theoretical background 2.1

The National Land Agency of the Republic of Indonesia

Badan Pertanahan Nasional Republik Indonesia (BPN RI) or the National Land Agency of the Republic of Indonesia (NLA) is a governmental institution, non-ministry, who is responsible to the President of the Republic of Indonesia in accordance with Peraturan Presiden Nomor 6 Tahun 2006 Tentang Badan Pertanahan Nasional (Presidential Decree Number 10/2006). In September 2013, this regulation was replaced with Peraturan Presiden Nomor 63 Tahun 2013 Tentang Badan Pertanahan Nasional Republik Indonesia (Presidential Decree Number 63/2013). And again in January 2015, this regulation was replaced with Peraturan Presiden Nomor 20 Tahun 2015 Tentang Badan Pertanahan Nasional Republik Indonesia (Presidential Decree Number 20/2015). Generally, the main task of NLA is supporting the President in land administration and land management at national, regional and sectoral levels, which is about land surveying and mapping, land rights and registration, land planning, land controlling, land acquisitioning for development of public interest, and land disputes assisting. In October 2014, NLA’s structure had been changed to be merged with another ministry called Kementerian Agraria dan Tata Ruang/Badan Pertanahan Nasional (Ministry of Land and Spatial Planning/National Land Agency). NLA has a vertical organizational structure (Fig.1): central office, regional offices at province level, and local offices at district level. A vertical organizational structure means NLA has its office representation in central, province and district level. Each level has different tasks and functions with central office as the policies coordinator.

BPN Pusat/NLA central office (Central office)

Kantor Wilayah BPN/Regional Office (Provincial level, 33 offices)

Kantor Pertanahan Kabupaten/Kota/ Local Office (District level, 440 offices)

Figure 1. The organizational structure of NLA. Source: adapted from Presidential Decree Number 20/2015.

6

In NLA central office, the Directorate of Basic Mapping - Deputy Survey, Measurement, and Mapping is the division that responsible to provide both basic map and high resolution satellite images for land registration. Every year they will buy high resolution satellite images based on the budget planning of previous year. However, they cannot provide high resolution satellite images continuously because of budget limitation. If other divisions in NLA need the satellite images, they will ask directly to the Directorate of Basic Mapping. In 2012, the President of the Republic of Indonesia gave an instruction that only Lembaga Penerbangan dan Antariksa Nasional/LAPAN (National Institute of Aeronautics and Space) is responsible to the provision of satellite imagery for all government departments in Indonesia. Therefore, officially in January 2015, LAPAN gave the high resolution satellite images to 9 departments in Indonesia included NLA. Thus, the Directorate of Basic Mapping - Deputy Survey, Measurement, and Mapping is still the coordinator of the high resolution satellite images (Fig.2).

Figure 2. The management of high resolution images in Indonesia and the National Land Agency of the Republic of Indonesia. Source: adapted from www.bpn.go.id and www.lapan.go.id. 7

2.2

Idle land controlling in Indonesia

Idle land controlling is one of priority programs in NLA which aims to make land right holders use their land according to the purpose of the entitlement. Moreover, it intends the land right holders to use their land for the welfare of the other people. Since 2010, the NLA started the idle land controlling based on Government Regulation Number 11/2010, Regulation of the Head Office of National Land Agency Number 4/2010, and Regulation of the Head Office of National Land Agency Number 9/2011. The subject of idle land controlling is the right holder who abandons the land, either individual or corporate. Furthermore, the main object of idle land controlling is land that has been granted with a right by the State: Hak Milik (freehold title), Hak Guna Usaha (right to cultivate), Hak Guna Bangunan (right to build), Hak Pakai (right to use), Hak Pengelolaan (right to manage) that does not used as the purpose of the entitlement as regulated in article 2 Government Regulation Number 11/2010. Based on Regulation of the Head Office of National Land Agency Number 4/2010, the idle land controlling has 4 main steps (Fig.3).

Indicated idle land inventory

Indicated idle land identification

Rights holder warning

Idle land determination

Figure 3. The idle land controlling main steps. Source: Regulation of the Head Office of National Land Agency Number 4/2010.

8

Firstly, indicated idle land inventory aims to collect the data of indicated idle land rights and update the indicated idle land database. Regional and local offices carry out the inventory by observing in the field and check the existing use of the land parcel. Regional and local offices will receive budget to do this activity every year for different land parcel. The inventory is also possible from people or public those see if there is indicated abandon land. During the inventory, officers of NLA will monitor the existing land cover together with the associated land use history. For instance, if officers discover bare soil or grassland as the existing land cover, they have to check whether it is already bare soil or grassland for a while or it is just there for a short period of time. Hence, the existing land cover and the land use history should be combined as a consideration. The result of inventory will be reported to central office as a basic to make idle land controlling planning both activity and budget planning for upcoming year. Secondly, based on the updated indicated idle land inventory, central office will determine the object and target of idle land controlling, thus the regional office will perform identification, both juridical and physical data. Officers will check all related documents with the land parcel and perform field identification: location and land use identification using Global Positioning System (GPS) handheld. Some indications showing the land owner already use the land parcel: the existing crop is the same as mentioned in the purpose of the entitlement and land clearance of land parcel. The results of this activity are detail report and land use map that should be bound as one book. Thirdly, Panitia C, a group consists of some determined officers; will have official discussion based on identification report. They decide whether the land owner will get a warning or not. The warning aims to warn the owner to cultivate the land based on the purpose of the right. The regional office will give the warning to the owner every month for three times in total. And by the end of every warning, officers will go to the land parcel and observe the progress whether the owner start cultivating the land or not. All of this progress will be reported in an official report for further step. Finally, based on final observation, regional office will send all documents needed to determine an idle land determination by the Head of NLA. The officers in central office will prepare the summary of an idle land controlling sent by regional office for the Head of NLA. After idle land determination by the Head of NLA, the ex land owner should leave the land parcel.

2.3

Oil palm plantation in Indonesia

In recent decades, there has been a large conversion of agricultural lands, rain forest, and formerly forested areas into oil palm plantations in Indonesia (Germer and Sauerborn, 2008; McCarthy and Cramb, 2009). This phenomenon follows the increase of world demand for cooking oil, vegetable oil, and bio diesel (Donald, 2004; Ramdani and Hino, 2013). Other factors that the high interest in oil palm plantations are: 1) oil palm is the highest-yielding vegetable oil crop (Donald, 2004), 2) low cost labor and production, 3) the availability of area to be planted, and 4) the suitable climate and conditions of Southeast Asia (Basiron, 2007; McCarthy and Cramb, 2009). The basic habitat for oil palm (Elaeis) is tropical rainforest with 1780–2280 mm annual rainfall and a temperature range of 24–30°C (Sheil et al., 2009). Table 1 shows the summary of favorable soil conditions for oil palm. 9

Table 1. Favorable soil conditions for oil palm plantations. Source: adapted from Turner and Gillbanks (1974). Property Terrain (slope degree) Effective soil depth (cm) Soil texture pH Peat depth (m)

Suitable condition 12 - 20 50 - 75 Deep sandy to clay loam 4.0 - 6.0 0 - 0.6

Indonesia is the largest producer of palm oil, with 9 million ton produced in 2010. In 2011, oil palm plantation reached 8.7 million hectares distributed over whole Indonesia (Ramdani and Hino, 2013). Sumatra and Borneo contribute 50% of oil palm over the world. Indonesia’s government plans to increase the production and explore oil palm plantation in East part of Indonesia (Ramdani and Hino, 2013). Moreover, some provincial and district government support the oil palm plantations (Suara, 2007). Oil palm plantations need some steps before they can be harvested. Oil palm nursery needs at least one year before planting in the plantation area (Basiron, 2007). Next, the oil palm usually will be planted at a 9 m by 7.5 m spacing and the resulting 148 palms per ha. The palms grow approximately 80 cm per year and will reach 20 m in 25 years (Sheil et al., 2009). Moreover, the fruit of oil palms can be harvested 2 – 3 years after planting (Basiron, 2007) or even 4 years after planting. Oil palm produces more during the rainy season, when palms receive more water (Feintrenie et al., 2010). The most productive period for an oil palm tree is between an age of 9 – 15 years (Sheil et al., 2009). After 25 – 30 years, the trees should be replaced (Basiron, 2007). Oil palm plantation appearance can be recognized from aerial photo by its fronds, i.e., large compound leafs (Fig.4). Red boxes show the oil palm plantation because it has regular spacing as mentioned by Sheil et al (2009) and its fronds have unique appearance (Turner and Gillbanks, 1974). The canopy of oil palm looks different compared to forest canopy that is shown by yellow boxes. Meanwhile, the blue box shows either young oil palm or a used oil palm plantation, it shows a regular line pattern and the edges of one plantation from another one.

Figure 4. Oil palm plantation in Kampar-Riau Province. Source: https://www.google.nl/ maps//place/Kampar+Regency,+Riau,+Indonesia. This image comes from QuickBird satellite image with a spatial resolution of 61 cm (SIC, 2014). 10

According to Turner and Gillbanks (1974), the ripeness of oil palm fruit harvested is a correlation between the weather, the length of time between harvesting rounds, and the minimum ripeness standard accepted. Thus, it will affect the harvesting frequency. For 7 – 8 years oil palm with the good ripeness, the best harvesting frequency/interval is 10 – 12 days. Turner and Gillbanks (1974) describe the process from oil palm fruits to oil. First, oil palm fruits will be received in the factory; it should be weighted in a good way to calculate field production. Then, the fruits will be stored in temporary storage, in which it should be stored as short as possible to minimize the degradation of oil quality. The next step is fruit sterilization using steam injection which aims to kill any microorganism, to loosen the fruits on the bunch, to soften the pericarps, and to dehydrate the fruits and partially dehydrate the nuts, as well as to shrink the kernel slightly. After the sterilization, the fruits will be separated from the bunch called stripping or threshing using stripper drum. The stripped fruits will continue to next step and the empty bunches are able to use for fertilizer or another use. Next, fruits will pass to the digesters, a real process of oil and kernel begins. The aims of digestion are: 1) to release oil from the pericarps cells by mashing them, 2) to raise the temperature of the mash to facilitate subsequent pressing, and 3) to drain away free oil, reducing the volume to be pressed. Next step is the extraction of crude oil palm. The pericarps will be pressed, thus the crude oil will come out. Next, it will be purified since the crude oil might still contain small materials. The crude oil will pass to a heated crude oil tank. It will purify and clarify the oil from wet oily sludge and water. Finally, the palm oil leaves the purification system and goes to storage tanks which are on a temperature of 45°C to keep the oil in an ideal condition and prevent from cooling. After this one, oil is ready for shipment.

2.4

Disturbance detection in Landsat time series using BFAST Monitor

Break For Additive Seasonal and Trend (BFAST) algorithm has been used to decompose remote sensing NDVI time series data into additive components. “BFAST integrates the iterative decomposition of time series into trend, seasonal and noise components with methods for detecting and characterizing changes (i.e. breakpoints) within time series” (Verbesselt et al., 2010a). BFAST is able to detect abrupt changes in the seasonal and trend components of time series and characterize gradual and abrupt ecosystem changes by deriving the time, magnitude, and direction of the change within the trend component of the time series. It has shown the capability to detect breakpoints in the linear trend: forest stands development (Lambert et al., 2011), drought disturbance in Somalia (Verbesselt et al., 2012), phenological change in vegetation indices (Verbesselt et al., 2010b), vegetation change and trend in U.S. army training Fort Riley-Kansas (Hutchinson et al., 2015). However, BFAST change detection approach is not developed to detect disturbances in recently acquired data to assess the stability. In 2012, a method to detect changes in near real-time has been developed based on BFAST-type season-trend model to assess the stability of linear regression models (Verbesselt et al., 2012). This method is called the BFAST monitor. The BFAST monitor will detect disturbances in near real-time by comparing newly acquired data with a modeled stable period based on historic insitu or satellite time series data. The method allows the automated differentiation between normal and abnormal change in near real-time based on the BFAST concept (Verbesselt et al., 2012).

11

Figure 5. BFAST monitor plot output for single pixel time series. Source: http://dutri001. github.io/ bfastSpatial/.

Three steps for near real-time monitoring are (Fig.5): 1) identify a stable history period; 2) model the stable history period; and 3) monitoring if the new observations in the monitoring period are conform to the expected behavior of the history sample (Verbesselt et al., 2012). BFAST monitor as a near real-time disturbance monitoring approach has been tested in drought detection in Somalia by Verbesselt et al. (2012). It is also used for detecting forest disturbances in tropical forest of southern Ethiopia by DeVries et al. (2015). The method was applied to Landsat time series in tropical forests enabling a rapid response by detecting forest disturbances in near real-time. The advantages using BFAST are: there is no need to select a reference period, set a threshold, or define a change trajectory (Verbesselt et al., 2010a). Threshold values often produce misleading results due to different spectral and phenological characteristics of land cover types (Lu et al., 2004). Another advantage of using BFAST is that it is applicable to any time series data, e.g. Landsat and land cover types (Verbesselt et al., 2010a; DeVries et al., 2015). Furthermore, it is able to function as a change warning system and gives spatial and temporal information (Verbesselt et al., 2010a). Evaluating BFAST monitor that able to detect abrupt change, possible to apply in high temporal Landsat images and another advantages, BFAST monitor was proposed to detect land clearance in this study. However, BFAST monitor never been used to detect land clearance. Thus, detailed analyze was performed in this study.

12

Figure 6. The principle of maximum likelihood classifier. Source: Lillesand et al. (p. 556; 2008).

2.5

Land cover classification

Land cover on HGU land parcels are very important to be monitored periodically to make sure that the owner uses the land as mentioned in the purpose of the entitlement. Landsat images that are free, adequate spatial and high temporal resolution are possible to detect change (Zhu and Woodcock, 2014). Moreover, by the successful launch of Landsat 8 in February 2013 which had more new features, it is expected that Landsat 8 will give better contribution in land cover mapping (Jia et al., 2014). Therefore, Landsat is good choice to monitor land cover in idle land controlling, especially in indicated idle land inventory. However, Landsat is never been used by the NLA to monitor land cover of land parcels for idle land controlling purpose. Many applications use a land cover map and “land cover classification is one of the most studied topics in remote sensing” (Zhu and Woodcock, 2014). The objective of image classification is to assign all of pixels in the image to several classes or themes (water, forest, bare soil, etc) and to make a thematic map (Foody, 2002; Weng, 2010). Moreover, land cover map is able to predict the areal extent of various cover classes (Stehman and Czaplewski, 1998). Some methods to perform classifications are supervised classification, unsupervised classification, and time series based classification, e.g., Continuous Change Detection and Classification (CCDC) algorithm (Lillesand et al., 2008; Zhu and Woodcock, 2014). Supervised classification is classification in which the users use their knowledge of the class classification and this process starts by selecting training areas (Lillesand et al., 2008). One of the parametric rules in supervised classification and the common one is maximum likelihood (Hashim and Ahmad, 2014; Partoyo and Shrestha, 2013). “The maximum likelihood classifier quantitatively evaluates both the variance and covariance of the category spectral response patterns when classifying an unknown pixel” (Fig. 6) (Lillesand et al., 2008). On the other hand, the principle of time series based classification is that different land cover classes will have different shapes 13

for the estimated time series models. For instance, “CCDC algorithm uses the coefficients of time series models as the inputs for land cover classification. After the change detection process, each pixel will have its own time series models before and after any changes. By classifying the time series model coefficients, this algorithm can provide a land cover type for the entire time period for each time series model” (Zhu and Woodcock, 2014).

2.6

Accuracy assessment

The accuracy assessment is important after performing a digital land cover classification and is considered as an important part of spatial data (Congalton, 2001; Lillesand et al., 2008). “In thematic mapping from remotely sensed data, the term accuracy is used typically to express the degree of ‘correctness’ of a map or classification” (Foody, 2002). According to Stehman and Czaplewski (1998), three basic components of an accuracy assessment are: 1) the sampling design used to select the reference sample; 2) the response design used to obtain the reference land cover classification for each sampling unit; and 3) the estimation and analysis procedures. The primary reference sample/sampling units are pixels, polygon, and fixed area plots. Meanwhile, some sampling design are simple random, stratified random, cluster, and systematic sampling. Another sampling consideration is sample size. In general, 50 samples is minimum number to be included for each land cover category and the number might be increased by the increasing of area (Congalton, 2001; Lillesand et al., 2008). Some methods to calculate the accuracy are overall accuracy of a classification error matrix (confusion matrix) and Kappa coefficient (Congalton, 2001; Lillesand et al., 2008; Olofsson et al., 2014). The accuracies derived from error matrix are overall accuracy, user’s accuracy, and producer’s accuracy (Lillesand et al., 2008; Olofsson et al., 2014). “Overall accuracy is computed by dividing the total number of correctly classified pixels by the total number of reference pixels”, user’s accuracy is defined as the probability that a pixel labeled as a certain land cover class in the map is really this class, and producer’s accuracy is defined as the probability that a certain land cover of an area on the ground is classified as such (Lillesand et al., 2008). There are some studies which had different overall accuracy percentage to assess the quality of land cover classification. Thomlinson et al. (1999) set a target of an overall accuracy of 85% with each class is not less than 70% as an accurate. Meanwhile, Zhu et al. (2014) set 90% of overall accuracy as a high accuracy. However, many studies had overall accuracy below generally recommended in a value 85%, a large range between 40% and 100% (Foody, 2002). “This accuracy assessment is concerned as an informative one. However, it has a problem that some cases may have been allocated to the correct class purely by chance” (Pontius, 2000). Hence, to take into account the effects of chance agreement, there is a Kappa coefficient as a standard measure of classification (Foody, 2002). “Kappa analysis is a KHAT statistic (an estimate of Kappa), which is another measure of agreement or accuracy. It is a measure of the difference between the actual agreement between reference data and an automated classifier and the chance agreement between the reference data and a random classifier. The values can range from +1 to –1” (Congalton, 2001; Landis and Koch, 1977; Lillesand et al., 2008; Stehman, 1996; Stehman, 1997). 14

Kappa coefficient can be defined as,

Kappa =

observed accuracy − chance agreement 1 − chance agreement

Landis and Koch (1977) characterized the possible ranges for Kappa coefficient into six groups: 1) a value greater than 0.80 until 1 represents strong agreement/ almost perfect; 2) a value between 0.61 and 0.80 represents substantial agreement; 3) a value between 0.41 and 0.60 represents moderate agreement; 4) a value between 0.21 and 0.40 represents fair agreement; 5) a value between 0.00 and 0.20 represents slight agreement; 6) a value less than 0.00 represents poor agreement.

15

3. Materials and Methods 3.1

Study area

The study areas of this thesis are located at the Riau Province - Indonesia (Fig.7). The geographical coordinates of the study areas are between 02°25’00” N - 01°05’00” S and 100°00’00” E - 105°05’00” E. The total area of Riau Province is 87,023.66 km2 (Kementerian Dalam Negeri Republik Indonesia -Ministry of Home Affairs of the Republic of Indonesia, 2014b). Riau has a wet tropical climate, with an average rainfall ranging from 1,000 – 3,000 mm per year, with a dry season (April–September) and a rainy season (October–March) (Ramdani and Hino, 2013; Dephut, 2014).

: HGU land parcel

Figure 7. Study area of five HGU land parcels (red) in Riau Province – Indonesia. Source: NLA.

16

3.2

Materials

3.2.1 HGU land parcels and Kampar regency land cover map In this study, we used 5 HGU land parcels located in Riau Province (Table 2). All HGU land parcels were obtained from the National Land Agency of the Republic of Indonesia in shapefile (.shp) format with WGS 1984 Universal Transverse Mercator (UTM) zone 47 N as projected coordinate system. Table 2 provides a description about those land parcels. All HGU have validity period for 25 until 35 year and for confidential reasons, the true identities were not mentioned. We used alphabetic identity to name all land parcels instead of used numerical since NLA using numerical system to give identity to HGU land parcels. Until 2015, HGU A until HGU E is not targeted in idle land controlling. Therefore, we chose HGU with various years of publication date and did not pay attention with the rule of three year after the publication of land certificate in which these land parcels are still potential to be targeted in idle land controlling. We also used Peta Rupa Bumi of Kampar regency (Kampar land cover map) year 2014 with a scale of 1:50,000 from Badan Informasi Geospasial (BIG – The Geospatial Information Agency of the Republic of Indonesia) for validation when performing land cover classification. The map was in shapefile (.shp) format with World Geodetic System (WGS) 1984 UTM zone 47 N as projected coordinate system. 3.2.2 Satellite images datasets We used Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper plus (ETM+) images for BFAST monitor input, one Landsat 8 Operational Land Imager (OLI) for land cover classification, and SPOT 6 images for ground truth and validation. Landsat TM, ETM+, and OLI have a 30 meter spatial resolution and a 16 day revisit time (Lillesand et al., 2008; Roy et al., 2014). The TM has 7 spectral bands, ETM + has 8 bands including a panchromatic band, and OLI has 11 bands (Table3) (Lillesand et al., 2008; Chander et al., 2009; Roy et al., 2014). All Landsat images had WGS 1984 UTM zone 47 N as projected coordinate system.

Table 2. General overview of HGU land parcels. Source: NLA. HGU code HGU A HGU B HGU C HGU D HGU E

Publication year August 2003 March 2005 June 2007 May 2008 November 2003

Purpose Oil palm Oil palm Oil palm Oil palm Oil palm

Area (Ha) 23,318 226,863 501,251 309,946 22, 232

Location Kampar-Riau Kampar-Riau Kampar-Riau Kampar-Riau Kampar-Riau

17

Table 3. The overview of Landsat TM, ETM+, Landsat 8 spectral bands. Source: adapted from Chander et al. (2009), Lillesand et al. (p. 411; 2008), Roy et al. (2014), USGS (2014). Satellite

Sensor

Band

Nominal spectral location

Spectral range (µm)

Spatial resolution (m)

Landsat 4

TM

Band 1

Blue

0.45 – 0.52

30

Band 2

Green

0.52 – 0.60

30

Band 3

Red

0.63 – 0.69

30

Band 4

Near IR

0.76 – 0.90

30

Band 5

Mid IR

1.55 – 1.75

30

Band 6

Thermal IR

10.4 – 12.5

120

Band 7

Mid IR

2.08 – 2.35

30

Band 1

Blue

0.452 – 0.514

30

Band 2

Green

0.519 – 0.601

30

Band 3

Red

0.631 – 0.692

30

Band 4

Near IR

0.772 – 0.898

30

Band 5

Mid IR

1.547 – 1.748

30

Band 6

Thermal IR

10.31 – 12.36

60

Band 7

Mid IR

2.065 – 2.346

30

Band 8

Panchromatic

0.515 – 0.896

15

OLI

Band 1

Blue

0.43 – 0.45

30

TIRS

Band 2

Blue

0.45 – 0.51

30

Band 3

Green

0.53 – 0.59

30

Band 4

Red

0.64 – 0.67

30

Band 5

Near IR

0.85 – 0.88

30

Band 6

Shortwave IR

1.57 – 1.65

30

Band 7

Shortwave IR

2.11 – 2.29

30

Band 8

Panchromatic

0.50 – 0.68

15

Band 9

Cirrus

1.36 – 1.38

30

Band 10

Thermal IR

10.60 – 11.19

100

Band 11

Thermal IR

11.50 – 12.51

100

Landsat 5

Landsat 7

Landsat 8

ETM

18

Table 4. Monthly frequency of available Landsat data from 2008 to 2014 for path 127 row 60. Month Number of images

Jan.

Feb.

Mar.

Apr.

May

June

July

Aug.

Sep.

Oct.

Nov.

Dec.

5

11

8

10

11

12

10

10

9

12

11

5

Since October 2014, USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) provided Demand Interface to order bulk processed satellite images (https://espa.cr.usgs.gov/). All available Level 1 Terrain (corrected) (1T) Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper plus (ETM+) images with path 127 and row 60 and a cloud cover less than 100% and Normalized Difference Vegetation Index (NDVI) as the spectral indices ranging from 2008 until 2014 were chosen and downloaded for free from http://glovis.usgs.gov/. Since BFAST monitor was applied to each pixel of the Landsat images (DeVries et al., 2015), we took into account the geographic navigation accuracy of a particular pixel because it is important for determining the changes in time series datasets (Masek, 2001). Thus, we chose X and Y root mean square error (RMSE) less than 15 m (0.5 pixels) for all images, to get a precise geometric match among Landsat images (Vicente et al., 2008). After the RMSE selection, a total of 114 images were used: 26 images from Landsat TM 5 and 88 images from Landsat ETM+ (Fig.8). From the monthly distribution of the images in the selected period of 7 years, we can observe that the number of images during the rainy season is considerably lower, especially in January and December with only 5 images. The frequency of images for the other months is relatively the same (Table 4). One Landsat 8 for land cover classification was acquired on 12 June 2014 with 13.38% of cloud cover. We chose this date because our SPOT 6 image for ground truth/reference and validation of existing land cover classification was acquired on 17 June 2014. Hence, our Landsat image had the least different time with our SPOT image compared to the other Landsat images. Moreover, this Landsat image had the least cloud cover during 2014 for our study area.

Figure 8. The number of Landsat TM/ETM+ images from 2008 to 2014 for path 127 row 60 (Riau Province). 19

Table 5. The overview of SPOT 6 spectral bands. Source: adapted from Hopwood (2014). Band

Nominal spectral location

Spectral range (µm)

Spatial resolution (m)

Band 1

Blue

0.45 – 0.52

6

Band 2

Green

0.53 – 0.59

6

Band 3

Red

0.63 – 0.70

6

Band 4

Near IR

0.76 – 0.89

6

Band 5

Panchromatic

0.45 – 0.75

1.5

Besides Landsat images, we also used Satellite Pour l’Observation de la Terre (SPOT) 6 images for validation. SPOT 6 was launched on 12 September 2012. It has a daily revisit time, 4 multispectral bands with a 6 meter of spatial resolution and 1 panchromatic band with a 1.5 meter spatial resolution (Table 5) (Hopwood, 2014). We used ortho rectified SPOT 6 images which were acquired on 7 January 2013 and 17 June 2014 (Fig.9). All SPOT 6 images had WGS 1984 Universal Transverse Mercator (UTM) zone 47 N as projected coordinate system. These images were obtained from LAPAN (the National Institute of Aeronautics and Space of the Republic of Indonesia).

(a)

(b)

Figure 9. Oil palm plantation from SPOT 6 satellite images: (a) panchromatic band, (b) multispectral bands 1-2-3 composites. Source: LAPAN, acquired on 17 June 2014.

20

3.2.3

Software

Most pre processing and analysis related to NDVI time series, BFAST monitor change magnitude and breakpoint detection and some graphs were produced using free software R 3.0.3 version. Most activities related to Landsat 8 maximum likelihood supervised classification and analysis was performed using Erdas Imagine 2014. We used Microsoft Excel for accuracy calculation and graphical plotting. For validation we used both Erdas Imagine 2014 and ArcGIS 10.2 software package (ESRI, Inc.), for cartography and map production, we used ArcGIS 10.2.

3.3

Data analysis methods

This part describes the data processing and analysis to answer all research questions. The flow chart (Fig. 10) shows the steps from the acquisition of satellite images of the study area to the extraction of the required information. By these procedures, we have answered the possibility of using Landsat based change detection method to detect land clearance and the accuracy of Landsat 8 to distinguish land cover 4 years after the publication of idle land controlling regulation.

21

Figure 10. Flowchart of the method for developing a method to characterize land use history using Landsat time series as an idle land early detection in Indonesia 22

3.3.1

Preprocessing

Preprocessing of satellite images is important to get good images which can be compared in time. Some of the Landsat preprocessing steps are radiance calibration, atmospheric correction, mosaic, and subset (Jia et al., 2014). In this study, the preprocessing stages were: 1) extraction of processed NDVI; 2) application of CFmask; 3) subset/cropped clean NDVI based on study area; and 4) creation of clean NDVI time series stack. All Landsat images for path 127 and row 60 in the period 2008 until 2014 have been processed for atmospheric correction using Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Huang et al., 2009). Atmospheric correction was needed for change detection and classification of multi temporal images. Without atmospheric correction those uncorrected images will lead to misinterpretations and wrong classifications (Song et al., 2001). Moreover, clouds and cloud shadows were detected in every image using the object-oriented Fmask algorithm and subsequently masked out of the images (Huang et al., 2010).

3.3.2

The determination of forest mask

In order to focus only to the pixels to be included in the disturbance monitoring algorithm, we created forest masks for all the study areas. DeVries et al. (2015) used a supervised classification of Landsat images in the beginning of the monitoring period to create the forest mask. In this study we created forest masks by classifying multispectral SPOT 6 images acquired on 7 January 2013 as the beginning of the monitoring period using a maximum likelihood supervised classification method. We made some initial classes: forest, oil palm, bare soil, and grassland based on visual interpretation as described in section 3.3.5 (Table 8), thereafter we aggregated those classes into two final classes: forest and non-forest. We resampled the 6 m SPOT spatial resolution forest mask to 30 m of Landsat spatial resolution. All non forest pixels according to the 2013 forest mask were masked out of the images and produced a clean forest NDVI time series stack. 3.3.3

Land clearance detection in Landsat time series using BFAST Monitor

In order to detect the disturbances, we applied the BFAST monitor to each pixel of a Landsat NDVI time stack. BFAST monitor has been applied for drought detection in Somalia by Verbesselt et al. (2012) and forest disturbances in the tropical forest of southern Ethiopia by DeVries et al. (2015). According to Verbesselt et al. (2012), the three steps for near real-time monitoring are: 1. Defining the history period: data that has been acquired and which will be analyzed for stability in order to model normal vegetation dynamics. The length of history period should be more than 2 years (>2years) (Verbesselt et al., 2012). Thus, we defined Landsat time series data from 2008 until 2012 as our stable history period to determine normal prediction response in the monitoring period (2013- 2014). This monitoring time was corresponding with the publication year of idle land controlling regulation in 2010 as the starting point to monitor land cover of HGUs in Indonesia.

23

2. Monitoring period, the period representing new observed data that has been monitored for existence of disturbance. Our monitoring time in our algorithm started in 1st January 2013 and ended in 1st July 2014. This monitoring time was chosen due to availability of our SPOT 6 images for validation. 3. The season-trend model is estimated based on a stable history period and predicts the normal data variation. The most important BFAST monitor outputs are the disturbance time and the change magnitude per pixel. “The magnitude (ΔNDVI) is defined as the median of the residuals within a monitoring period between observations and expected values (based on the history period model)” (DeVries et al., 2015). It means that all pixels are assigned a value whether a breakpoint has been detected or not (Fig. 11). By observing the disturbance in time series, it is possible to estimate the time and location of the breakpoints (Verbesselt et al., 2010a). However, the breakpoints are sometimes considered with positives magnitudes or near zero. Those positive breakpoints do not represent deforestation, but an increase in NDVI that might represent regrowth (DeVries et al., 2015). Thus, breakpoints are not good enough to determine changes. Breakpoints in combination with the magnitude value are possible solutions to determine the changes processes. Moreover, a short monitoring period (e.g., one year) could avoid the effect of an increased number of observations before and after the change occurred to magnitude calculation. “Change magnitude in longer monitoring periods would be affected by follow-up land use” (DeVries et al., 2015). In this study, we combined the breakpoints and magnitudes value to determine land clearance and no land clearance class. We used breakpoints and set the threshold to 0, -0.025, -0.05, -0.075, -0.1 for the magnitude value to evaluate which threshold would obtain the best accuracy in detecting land clearance for five study areas. Meanwhile, the no land clearance class was determined by non-breakpoint pixels with zero or positive magnitude value

Figure 11. Demonstration of the sequential monitoring approach used by DeVries et al., 2015. A breakpoint (dotted red line) and magnitude (M) are shown on each plot. The residuals, defined as the difference between observed values (black dots) and expected values (blue line) are shown as broken grey lines in each monitoring period. The solid black line shows both the end of the history period and the beginning of the monitoring period. Source: DeVries et al., 2015, http://dx.doi.org/10.1016/j.rse.2015.02.012. 24

3.3.4

Land clearance validation and accuracy assessment

It was very important to validate the land clearance detected by the BFAST monitor to know how sensitive the BFAST monitor is in detecting land clearance. We performed a stratified random sampling to validate both classes following DeVries et al. (2015). In case of the land clearance class, we validated all pixels as the result of different magnitude thresholds (5 strata) by using SPOT images acquired on 7 January 2013 and 17 June 2014. Meanwhile, for no land clearance class, we sampled between 19 and 100 pixels for each study area and then validated those pixels using SPOT images of years 2013 and 2014. For both classes, we made a confusion matrix to assess the accuracy of the BFAST monitor in detecting land clearance and no land clearance. For every sample design protocol, the response needs a specific agreement when comparing the result and the references data (Stehman et al., 2003). Hence, we chose the cross point between user’s and producer’s accuracy following DeVries et al. (2015) in order to get best threshold for detecting land clearance. We considered the cross point as an unbiased estimator variance which would minimize the biased (underestimated) standard errors in the accuracy assessment (Olofsson et al., 2014; Stehman et al., 2003; Stehman et al., 2012). 3.3.5

Land cover classification and analysis

Multiple methods have been used to perform land cover classification as mentioned in section 2.5. We performed maximum likelihood supervised classification for the available Landsat 8 image to evaluate the land cover of HGU land parcels 4 years after the publication of idle land controlling regulation (Regulation of the Head Office of National Land Agency Number 4/2010). According to Vicente et al. (2008), a precise geometric match in Landsat image could be obtained by choosing an image with X and Y RMSE less than 15 m (0.5 pixels). Hence, we chose one Landsat 8 image that was acquired on 12 June 2014 with image quality 9, a cloud cover percentage of 13.38, corresponding geometric RMSE model of 8.365 m, geometric RMSE model Y of 5.508 m, and geometric RMSE model X of 6.296 m. According to Lillesand et al. (2008), the choice in band combination is depending on personal preferences. All bands were considered useful for land cover mapping, for instance band 1 for coastal, band 2 for distinguishing soil from vegetation, band 5 for biomass content, band 9 for cirrus cloud contamination, etc (Hurd and Civco, 2009; Jia et al, 2014; USGS, 2013). Hence, we stacked six surface reflectance bands: band 2, 3, 4, 5, 6, and 7, as input for the classification. The corresponding spatial resolution is 30 m. Band 2 has a bandwidth of 0.45 µm - 0.51 µm and belongs to the green part of the electromagnetic spectrum and is used for bathymetric mapping, to distinguish soil form vegetation and deciduous forest from coniferous vegetation. Band 3 with a bandwidth of 0.53 µm - 0.59 µm is used to emphasize peak vegetation which is useful for assessing plant vigor. Band 4, 0.64 µm - 0.67 µm in the red part of the electromagnetic spectrum is used to discriminate vegetation slopes. Band 5 belongs to the near infrared (NIR) part of the spectrum and has a band width of 0.85 µm - 0.88 µm. It emphasizes the biomass content and shorelines. Band 6 we used to discriminate the moisture content of soil and vegetation. This band is included in the shore-wave Infrared and has a band width of 1.57 µm - 1.65 µm. This band has as additional property that it penetrates through thin clouds. Band 7 is also included in the short wave infra red (SWIR) and has a band width of 2.11 µm - 2.29 µm. It has been used to improve the estimation of the moisture content of the soil and vegetation. Just like band 6, this band is able to penetrate through thin clouds. 25

Table 6. The descriptions of land cover classes for Landsat classification. Source: adapted from Badan Standarisasi Nasional – The National Standardization Agency of the Republic of Indonesia (2010). Land cover class Oil palm Forest

Grassland Bare soil

Description Homogeneous plantations planted with oil palm, which is not replaced by the other types of plants for at least two years. Forest that grows on dry land habitats, which can be low lands, hills, and mountainous forest or highlands of tropical forests. The forest is still compacted and no human intervention have taken place e.g., forest logged. Open area that is dominated by a variety of heterogeneous grass. Land area without any vegetation.

Furthermore, related to atmosphere correction, Song et al. (2001) and Huang et al. (2002) mentioned that it is not necessary to perform atmosphere correction when performing classification for only one single date Landsat image. However, our surface reflectance bands have been processed for atmospheric correction using the L8SR algorithm which involves the solar zenith and view zenith angles in the calculation (USGS, 2015a; USGS, 2015b). One step in supervised classification is the determination of the reference data used for training and accuracy assessment (Jia et al., 2014; Yuan et al., 2005). Yuan et al. (2005) used aerial photos and Jia et al. (2014) used Google Earth for those purposes. Meanwhile, we used the SPOT 6 image acquired on 17 June 2014 as the ground truth/guidance to determine the training areas for classification and accuracy assessment. First, we made reference classes based on SPOT 6 image interpretation. According to Olson (1960), the basic characteristics to interpret images are shape, size, pattern, tone (hue), texture, shadows, site, association, and resolution (Table 7). Furthermore, textural information is very important in visual image interpretation (Zhang et al., 2003) and panchromatic images are able to give good topographic information (Shaker et al., 2012). Hence, in order to make reference classes, we referred to land cover published by The National Standardization Agency of the Republic of Indonesia as described in Table 6, the basic characteristics of image interpretation as described in Table 7, and literatures that mentioned about how to identify land cover from high resolution images (Table 8).

Table 7. The basic characteristics of image interpretation. Source: adapted from Olson (1960), Green (2000), and Lillesand et al. (p. 191-195; 2008). Characteristic Shape Size Pattern Tone (hue) Texture

Description refers to general form, configuration, or outline of individual objects refers to objects scale relates to the spatial arrangement of objects refers to the relatives brightness or color of objects the frequency of tonal change, produced by an aggregation of unit features 26

Table 7. (Continued) Characteristic Shadow Site Association Resolution

Description refers to the an area where light from a light source is obstructed by an object refers to topographic or geographic location refers to the occurrence of certain features in relation to the others refers to the smallest size of objects that still can be distinguished from the others

Table 8. The determination of reference classes based on SPOT 6 images visual assessment. Source: adapted from BSN– The National Standardization Agency of the Republic of Indonesia (2010), CRISP (2001), and Mongabay (2012). Land cover Visual assessment Visual interpretation class Oil palm

- the shape of the leafs has a regular pattern - the surface texture has a medium density - the plants are inside the rectangle, associated with oil palm fields - dark color indicates vegetation (low reflectance)

Forest

- the surface texture is high density/rough - irregular shape - dark color indicates vegetation (low reflectance)

27

Table 8. (Continued) Land cover class

Visual assessment

Grassland

- the surface texture has a soft density - homogeneous color in open areas - the color is brighter than for forest and oil palm

Bare soil

- bright color due to high reflectance

Visual interpretation

In total, we used four land cover classes based on our visual interpretation of SPOT image and the Indonesian National Standardization for land cover classes scale 1:50,000 (BSN- Badan Standarisasi Nasional – The National Standardization Agency of the Republic of Indonesia, 2010). The four classes were: oil palm, forest, grassland, and bare soil (Table 6 and Table 8). Yuan et al. (2005) made five until twenty training areas for every class in their Landsat land cover classification in Minnesota. Based on that, we made for every class two until ten training areas through polygons of the AOI’s (area of interest). The number of pixels for every AOI varied, ranged from 10 until 35 pixels. Thereafter, we assessed the signature plot and merged it into one class. 3.3.6

Land cover classification accuracy assessment

In order to assess the land cover classification accuracy, we performed a visual observation and quantitative classification accuracy indicators. We validated the classification result for the five study areas with SPOT 6 images and the Kampar land cover map, which has a scale of 1: 50,000. The purpose of this validation was to check whether the land cover map of the five study areas corresponded to the correct land cover class or not. For the first accuracy assessment, we performed random sampling and used pixel as the sampling unit following Jia et al. (2014), DeVries et al. (2015), and Zhu et al. (2014). Simple random sampling is more easily modified during the accuracy assessment and “in general, the simple random selection protocol will better satisfy the desirable design criteria and is the recommended option” (Olofsson et al., 2014; Stehman et al., 2012). Furthermore, 50 samples is the minimum number to be included for each land cover category, this number might increase by an increasing of the area. The number sample size can also be determined by the importance of a 28

certain land cover (Congalton, 2001; Lillesand et al., 2008). In the validation performed by Zhu et al. (2014), 500 pixels were chosen to assess the change detection accuracy. Meanwhile, Jia et al. (2014) used 50% of the number of pixels in their training areas as their sample for each land cover class. Hence, we used 800 pixels in total for all classes in the five HGU land parcels to assess the land cover classification accuracy following the procedure of Jia et al. (2014). In order to check the real land cover based on SPOT 6 images acquired on 17 June 2014, we made visual observation based on Table 8 and a double check: first, we overlaid the land cover map with SPOT 6 images in ArcGIS; second, we linked three images: a panchromatic image, a multispectral image with 4-3-2 composite and the classification result. This process has been done in ERDAS imagine 2014. Instead of making similar spatial resolution for both datasets, we used point’s coordinate to determine the sample and the real land covers. Therefore, we could keep the quality of SPOT images for validation. Next, we made a confusion matrix to analyze the accuracy between Landsat land cover classification and SPOT 6 images for five HGU land parcels. For the second validation, we overlaid the Kampar land cover map (shapefile format) and land cover classification of the five study areas. We performed visual validation using Arcgis and used transparency in display to make the analysis easier.

29

4. Results 4.1

Observing land clearance in HGU land parcels

4.1.1 Landsat multi-temporal NDVI profile We produced and analyzed the NDVI time series stack to obtain important information which was necessary to support idle land controlling especially the land cover history of HGU land parcels. The publication year of idle land controlling regulation in 2010 was the starting point to monitor land cover of HGUs in Indonesia. For all HGU certificates which were published after 2010, three years after the publication of the certificate, the HGU parcels will be monitored whether the owner has acted as mentioned in the purpose of the entitlement. Meanwhile, for all HGU certificates which were published before 2010, NLA will still monitor whether the owner used the land like mentioned in the purpose of the entitlement and those parcels can still be targeted in idle land controlling programs. Hence, the land use history plays an important role for NLA to determine idle land controlling targets. We plotted the annual median NDVI for the years 2008 until 2014 to observe the yearly median NDVI value of the whole study area. It gave the opportunities to analyze what happens in the study area based on this plot (Fig. 12; Appendix I Fig. 22, 23, 24, and 25). For instance, we can observe a general trend of NDVI value for HGU parcel D (Fig. 12). We focused on the area included in the blue circle in the middle of the plot for the purpose of our analysis. The certificate of HGU D was published in May 2008. The NDVI value from year to year in some of the areas was quite fluctuating (blue circle). In the picture of 2008, in which the certificate was published, we observed that the land parcel had a rectangular pattern (indicated by the arrows). It was more clearly shown in the picture of 2009. The value of the NDVI at the start of the history period in 2008 was quite low with values around 0.2 - 0.3. For the years 2009 and 2010 an increase in the NDVI can be observed with corresponding NDVI values of 0.4 and 0.5 respectively. Unfortunately, the quality of the image for 2011 was not good enough to analyze this picture. However, if we follow the general trend in the NDVI values, the value for the year 2011 would be around 0.5 - 0.6. In 2012, we can observe an increase in the NDVI to values around 0.6 - 0.8. For 2013 and 2014, a decrease in the NDVI values for those years can be observed to values around 0.3 - 0.4. A possible explanation of this decrease in the NDVI values can be harvesting of a certain crop or land clearance. This general pattern can also be observed for other areas in this HGU parcel, i.e. the area in the north-east of the parcel had the same trend, however an increase in the NDVI values for 2014 can be observed. Overall, we can see that the NDVI values stay relatively high with values around 0.8 for the other areas in HGU parcel D (Fig. 12). In HGU A (Appendix I Fig. 22), we cannot find any rectangle pattern like in HGU D. However, we can observe a fluctuative NDVI (shown in red circle). The NDVI had a relatively low value around 0.55 in 2008. In 2009, the same area had an NDVI of 0.7 and for 2010, it increased further to values around 0.8. In 2011, the NDVI values dropped to an average value around 0.6 for the area shown in the red circle. In 2012, due to the fact that SLC was off for Landsat 7, we did not have a full picture of the area. In 2013 and 2014, the NDVI in the selected red circle was quite stable with values in the range of 0.6-0.7. This pattern with fluctuative NDVI values we also see in other areas within HGU A, e.g., the area in the southeast corner of the picture. 30

0.8 0.7 0.6 0.5 0.4 0.3 0.2

0.8 0.7 0.6 0.5 0.4

0.8 0.7 0.6 0.5 0.4

0.8 0.7 0.6 0.5 0.4 0.3

0.8 0.6 0.4 0.2 0

0.8 0.6 0.4 0.2

0.8 0.7 0.6 0.5 0.4 0.3

Figure 12. Annual median NDVI of HGU D for period 2008 until 2014. In HGU B (Appendix I Fig. 23), we also do not recognize any rectangle pattern and the NDVI values were high for almost the whole time for the period 2008 until 2014. However, in 2011, it had a decrease NDVI (shown by red circle). Moreover we recognize some small areas, which had a low NDVI value throughout the whole time period in this pixel. Those low NDVI values can be caused by bare soil or buildings. In HGU C (Appendix I Fig. 24), it was almost the same with HGU B, the NDVI value was high for the whole time in the period 2008 until 2014. However, we can observe that in 2011 (shown by red circle) there was a decrease in NDVI value. The NDVI for almost the whole area is constant with values around 0.7-0.8. In HGU E (Appendix I Fig.25), the most significant changes were observed in the east part of this land parcel. In the period 2012 to 2013, there was a decrease NDVI value (shown by red circle). In the period 2008-2012, the NDVI values for this area were quite stable with NDVI values around 0.7. After 2012 there was a strong decrease in NDVI values, this can be possibly caused by land clearance. The NDVI values for 2013 and 2014 in the red circle were around 0.4-0.5. However, in the western part of this HGU parcel an opposite trend can be observed. In this area, lower NDVI values were observed for the years 2008 until 2011, with values around 0.6. For 2012 we 31

did not have any information and for 2013 and 2014 an increase in the NDVI values was observed with values around 0.75-0.8. 4.1.2 Land clearance detection using BFAST monitor The history and monitoring period are two of the three steps for near real time monitoring used to model the normal vegetation dynamics and represents new observed data that has been monitored for existence of disturbance (Verbesselt et al., 2012). Table 9 shows the mean and standard deviation clear sky pixel observation of the forest NDVI time series stack for the history (1st January 2008 to 31st December 2008) and monitoring period (1st January 2013 to 1st July 2014), corresponding to 81 and 26 Landsat images respectively. In the history period, we can observe that HGU C had the highest mean clear sky pixel observation, corresponding to 25.32. Meanwhile, HGU A had the lowest value, with a value of 14.78. For HGU B, D, and E, the mean clear sky pixel observations were 24.11, 23.89, and 16.49 respectively (Table 9). The results in the monitoring period were different from the history period. HGU B was the area which had the highest mean clear sky pixel observation, corresponding to 6.37. Meanwhile, HGU E had the lowest one, corresponding to 5.68. For HGU D, A and C the mean clear sky pixel observations were 6.32, 6.22, and 6.21 respectively (Table 9). We applied BFAST monitor to clean forest NDVI time series stack for five study areas resulting three raster layers for each study area: 1) breakpoint timing, 2) change magnitude, and 3) error flag (Appendix III Fig. 28, 29, 30, 31, and 32). Breakpoint timing is shown in decimal year format. Meanwhile, “the magnitude is defined as the median of the residuals within a monitoring period between observations and expected values (based on the history period model)” (DeVries et al., 2015). It means that all pixels are assigned a value whether a breakpoint has been detected or not. Meanwhile, “error flag is possible values of 1 if an error has been encountered in a particular pixel, and NA where the algorithm was successful (or a mask had already been applied previously)” (Dutrieux et al., 2014).

Table 9. The mean and standard deviation clear sky pixel observation of the forest NDVI time series stack for history and monitoring period for five HGU land parcels.

HGU A HGU B HGU C HGU D HGU E

History period (01-01-2008 to 31-12-2012)

Monitoring period (01-01-2013 to 01-07-2014)

(Total Landsat scenes: 81)

(Total Landsat scenes: 26)

Mean observation

Standard deviation observation

Mean observation

Standard deviation observation

14.78 24.11 25.32 23.89 16.49

1.61 1.92 2.13 2.72 2.48

6.22 6.37 6.21 6.32 5.68

1.35 1.28 1.08 1.35 0.95

32

Table 10. The number of pixels for land clearance class for five different magnitude thresholds. Magnitude threshold 0 -0.025 -0.050 -0.075 -0.1

HGU A

HGU B

HGU C

HGU D

HGU E

43 36 31 27 22

138 82 50 24 19

104 93 83 31 15

38 35 32 27 22

25 21 20 12 9

After we got five outputs of BFAST monitor for five study areas, next we used breakpoints and negative magnitude value to determine land clearance. We set the threshold of 0, -0.025, -0.05, 0.075, -0.1 to magnitude value in order to evaluate which threshold would obtain the best accuracy in detecting land clearance for five study areas (Fig. 13, Fig 14; Appendix IV Fig. 33, 34, 35, and 36). Furthermore, Table 10 shows the number of pixels for land clearance class for five different magnitude thresholds in which lower thresholds result in a lower number of pixels identified as land clearance for five HGU land parcels.

33

2013

MT = 0

MT = -0.050

MT = -0.025

MT = -0.1

MT = -0.075

2014 2014

Figure 13. Land clearance detected by BFAST monitor for HGU B. The upper five maps display the land clearance class determined by combining the breakpoints and five different magnitude thresholds (MT). The SPOT 6 image for the five upper maps was acquired on 7 January 2013. The lower map shows the SPOT 6 image acquired on 17 June 2014. This image shows the land cover of HGU B in 2014 and was used to validate the detected land clearance class of the upper five maps. 34

Figure 14. Example of a land clearance detected by BFAST monitor for single pixel time series in HGU B. The upper graph shows the output of the BFAST monitor for a single pixel shown by the red arrow. The breakpoint time indicates that there was a disturbance in October 2013 (296). The lower two pictures reveals that land clearance occurred for that specific pixel in the period 2013 and 2014.

4.1.3 Land clearance validation and accuracy assessment In order to assess the accuracy of BFAST monitor in land clearance detection, we sampled all pixels resulted in section 4.1.2 Table 10. Meanwhile, for no land clearance class, we sampled non breakpoint pixels with zero or positives magnitudes value which means this class was not affected by the magnitude threshold for land clearance class. Thus, we sampled once for every HGU resulting the number of pixels of 19, 100, 100, 60, and 56 for HGU A, B, C, D, and E respectively. After determining the land clearance and no land clearance class, we formed confusion matrix based on visual assessment validation using SPOT 6 images for the year 2013 and 2014 (Table 11; Appendix IV Table 15, 16, 17, and 18). Furthermore, Figure 15 shows the influence of the magnitude threshold setting on the user’s, producer’s and overall accuracy for detecting land clearance in HGU A, B, C, D, and E.

35

Table 11. The BFAST monitor accuracy assessment in detecting land clearance for HGU B. Land clearance = breakpoint and magnitude < 0 Land clearance No land clearance Total 30 108 138 Land clearance 6 94 100 No land clearance 36 202 238 Total Producer’s (%) 83.33 46.53 Overall (%) Land clearance = breakpoint and magnitude < -0.025 Land clearance No land clearance Total 30 52 82 Land clearance 6 94 100 No land clearance 36 146 182 Total Producer’s (%) 83.33 64.38 Overall (%) Land clearance = breakpoint and magnitude < -0.05 Land clearance No land clearance Total 30 20 50 Land clearance 6 94 100 No land clearance 36 114 150 Total Producer’s (%) 83.33 82.46 Overall (%) Land clearance = breakpoint and magnitude < -0.075 Land clearance No land clearance Total 20 4 24 Land clearance 6 94 100 No land clearance 26 98 124 Total Producer’s (%) 76.92 95.92 Overall (%) Land clearance = breakpoint and magnitude < -0.1 Land clearance No land clearance Total 17 2 19 Land clearance 6 94 100 No land clearance 23 96 119 Total Producer’s (%) 73.91 97.91 Overall (%)

User’s (%) 21.74 94 52.10 User’s (%) 36.58 94.00 68.13 User’s (%) 60.00 94.00 82.67 User’s (%) 83.33 94.00 91.94 User’s (%) 89.47 94.00 93.27

In general, we had the same pattern for five HGU land parcels (Fig. 15): 1) a decrease in the overall accuracy and the user’s accuracy with the magnitude threshold approaching 0, thus the maximum accuracies were reached for a magnitude threshold value of -0.1; and 2) the producer’s accuracy had an opposite trend in which increased with the magnitude threshold approaching 0. We chose the cross point of user’s and producer’s accuracy as an unbiased estimator for the variance which minimizes the biased (underestimated) standard errors in the accuracy assessment. HGU A had a cross point with magnitude -0.075. In this case, the user’s and producer’s accuracy were 89.15% and the overall accuracy was 82.61% (Fig.15). In HGU B, the user’s and producer’s accuracy crossed at a point with magnitude of -0.070 corresponding to 79% and an 89% of overall accuracy (Fig. 15). In HGU C, those accuracies crossed at a point with magnitude of -0.080, corresponding to 53% and an 81% for overall accuracy (Fig. 15). In HGU D, the cross point was found at a point with a magnitude of -0.060. The user’s and producer’s accuracy were 67% and the corresponding overall accuracy was 78% (Fig. 15). Meanwhile, for HGU E, the user’s and producer’s accuracy crossed at a point with magnitude of -0.075. At this cross point, user’s and producer’s accuracy were 66.67% and overall accuracy was 88.24% (Fig. 15).

36

Figure 15. Overall accuracy (OA), producer’s accuracy (PA), and user’s accuracy (UA) of land clearance detection as a function of magnitude for HGU A, B, C, D, and E. 37

4.2

Observing land cover classification

4.2.1 Land cover based on Landsat maximum likelihood supervised classification It is very important to monitor the land parcels in Indonesia periodically, to make sure that the owner uses the land as mentioned in the purpose of the entitlement. This monitoring will support idle land inventory in all provinces in Indonesia. We performed a maximum likelihood supervised classification for Landsat 8. This has been done to evaluate the land cover of HGU land parcels 4 years after the publication of the idle land controlling regulation. We used six surface reflectance bands of Landsat 8: band 2, 3, 4, 5, 6, and 7. In order to remove cloud contamination, we applied the CFmask to each band then stacked them together as an input for classification (Fig. 16). The maximum likelihood classification has been done to evaluate the land cover of HGU land parcels 4 years after the publication of the idle land controlling regulation. Figure 16 shows a comparison for the Landsat image before and after applying the CFmask. It is clearly shown by applying the CFmask were able to produce a clean image of the area. We cropped our Landsat image for our five HGU land parcels and made training areas based on a SPOT 6 image, which was acquired on 17 June 2014. This resulted to the fact that we could use in total 1597 pixels as our reference data (Table 12). The corresponding numbers of pixel and the area in hectare for all the study areas are shown in Table 13. The land cover classes recognized after performing a maximum likelihood supervised classification were oil palm, forest, grassland, and bare soil (Fig.17). After the classification, we observed that not all the HGUs had all mentioned land cover classes in their area. There were only two plots, HGU A and HGU B, which have all the four land cover classes (Table 12, Fig.17). For HGU A 9.36% of the area was classified as oil palm, which corresponded with 28 pixels. In this plot 10.71% of the area was classified as forest, 58.20% as grassland, and 21.73% as bare soil (Table 13). Whereas for HGU B, in which also four land cover classes have been recognized, the area was covered for 73.77% by oil palm, 8% by forest, 10.19% by grassland, and 8.04% by bare soil (Table 13).

Figure 16. The comparison of Landsat before (left) and after (right) applying CFmask in 4-3-2 composites. Both pictures show the same area (red arrow shows the same road). Left picture was a Landsat image before applying the CFmask and the right picture shows the one after applying the Cfmask.

38

Table 12. The number of pixels in each class type as reference data for training and validating the Landsat classification for the five study areas. HGU A HGU B HGU C HGU D HGU E Total

Oil palm 22 236 350 129 57 794

Forest 12 65 0 0 0 77

Grassland 65 36 35 289 38 463

Bare soil 35 40 21 151 16 263

Total 134 377 406 569 111 1597

For HGU C three land cover classes have been recognized, namely: oil palm, grassland, and bare soil. In this plot, oil palm was the dominant one with 72.23%. The other land cover classes cover 18.97% of the area for grassland and only 8.80% for the bare soil (Table 13). In HGU D, only three land cover classes have been recognized, those were: oil palm, grassland, and bare soil. In this plot, grassland and bare soil covered the largest part of the area with 67.37% and 25.52% respectively. Oil palm covered only 7.11% of the total area in plot D (Table 13). In HGU plot E, the area was in particular covered with grassland and oil palm. Those land cover classes cover respectively 49.65% and 42.50% of the area. Bare soil had the least area with only 7.85% (Table 13). From the five selected study areas, the oil palm class had the largest area with 583.830 ha and forest had the smallest area with 21.690 ha in total. Meanwhile, the area of grassland and bare soil were 368.550 ha and 154.890 ha respectively (Table 13). However, those calculated areas were the results of the classification, which still consisted of misclassified land cover classes in which one land cover class was mixed with another land cover.

Table 13. The number of pixels and the area (ha) after performing a land cover classification for the five selected study areas. HGU A HGU B HGU C HGU D HGU E Total pixel Total in ha

Oil palm 28 1926 4158 256 119 6487 583.830

Forest 32 209 0 0 0 241 21.690

Grassland 174 266 1092 2424 139 4095 368.550

Bare soil 65 210 506 918 22 1721 154.890

Total pixel 299 2611 5756 3598 280 12,544 1,128.960

Total in ha 26.910 234.990 518.040 323.820 25.200

39

Figure 17. Land cover classification of HGU A, B, C, D, and E from Landsat 8 image of 12 June 2014. 40

4.2.2 Land cover classification validation and accuracy assessment 4.2.2.1 Land cover validation using SPOT 6 images We used a SPOT 6 image, which was acquired on 17 June 2014 for our first validation. The purpose of this validation was to check whether the land cover map of the five study areas felt into the classes we determined based on SPOT 6 visual interpretation or not. The number of samples we used for each HGU land parcel was determined by more or less half of the number of pixels used for the training areas for every class (Table 12). The total sample pixels used in the random sample design for the classification accuracy assessment were 397 pixels for oil palm, 39 pixels for forest, 232 pixels for grassland, and 132 pixels for bare soil. We made a confusion matrix from this validation and calculated the producer’s accuracy, user’s accuracy, overall accuracy, and Kappa coefficient (Table 14). The results of the accuracy assessment showed producer’s accuracies of 88.80%, 86.21%, 75.39%, and 100% for oil palm, forest, grassland, and bare soil respectively (Table 14). Meanwhile, the user’s accuracy were 87.90%, 64.10%, 83.19%, and 92.42% for oil palm, forest, grassland, and bare soil respectively (Table 14). We calculated the overall accuracy and the Kappa coefficient resulting in 86.13% and 0.78 respectively, which showed a substantial agreement. The producer’s accuracy is defined as the probability that a certain land cover of an area on the ground is classified as such (classified pixels) (Lillesand et al., 2008). From the confusion matrix in Table 14, we can observe that 2.29% of the oil palm was wrongly classified as forest and 8.8% as grassland by the producer. For the classification of forest, the error occurs when making a difference between grassland and forest. This can be seen in the percentage of 13.79% that has been classified as grassland even though they should have been classified as forest. For the classification of grassland, oil palm had the highest correspondence in which 18.75% of the grassland was misclassified as oil palm. The percentages classified as forest and bare soil was negligible in the classification of grassland. The classification for bare soil had a producer’s accuracy of 100%. The user’s accuracy is defined as the probability that a pixel labeled as a certain land cover class in the map is really belonging to this class (observed pixels) (Lillesand et al., 2008). The user’s accuracy of oil palm was 87.90% and was only influenced by grassland with 12.10% (Table 14). The user’s accuracy of forest was strongly influenced by oil palm and grassland with 23.07% and 12.82% respectively. Because of these misclassifications, the user accuracy of forest was low with only 64.10%. The user’s accuracy of grassland was 83.19%. Grassland was labeled as oil palm in 15.08% of the cases and only 1.72% as forest. The user’s accuracy of bare soil was quite high with 92.42% and was only influenced by grassland with 7.57%. Overall, grassland had the lowest producer’s accuracy and forest had the lowest user’s accuracy. The largest confusions were between oil palm and grassland, between forest and grassland and also between oil palm and forest.

41

Table 14. The accuracy assessment for the classification using SPOT 6 images for the five HGU land parcels. Oil palm Forest Grassland Bare soil Total Producer’s (%)

Oil palm

Forest

Grassland

Bare soil

Total

349 9 35 0 393

0 25 4 0 29

48 5 193 10 256

0 0 0 122 122

397 39 232 132 800

88.80

86.21

75.39

100

User’s (%) 87.90 64.10 83.19 92.42

Overall (%) 86.13

4.2.2.2 Land cover validation using Kampar land cover map scale 1:50,000 After assessing the accuracy using SPOT 6 images, we validated the classification result for five study areas with the Kampar land cover map of the year 2014 with a scale 1: 50,000. The purpose of this validation was to check whether the land cover map of the five study areas corresponded with the same land cover as on the Kampar land cover map or not. This has been done by overlaying the Kampar land cover map (shapefile format) with the land cover classification of the five study areas. We adjusted the land use legend in Kampar land cover map from Bahasa to English to make it easier to understand (Fig. 18). The complete area of the HGU plots A, C, and E completely overlapped with the plantation class of the Kampar land cover map (Fig. 18). Our oil palm class as the result of Landsat classification was considered as plantation. However, if we compare this with the percentages of oil palm we got from the supervised classification this did not completely correspond. For plot A, 9.36% of the area was classified as palm oil and for plot C and E these percentages were 72.23% and 42.50% respectively (Table 13). The classes of forest, grassland, and bare soil were not recognized by the Kampar land cover map. In HGU plot B, according to the Kampar land cover map, almost the complete area of the plot was covered with oil palm and only partly with forest and grassland (Fig. 18). This corresponds also with the classification results for HGU plot B. The classification showed a distribution of 73.76% oil palm, 8% forest and 10.18% grassland. However, the bare soil class was not recognized in the Kampar land cover map. HGU plot D belongs partly to grassland, partly to forest and partly to oil palm plantation according to the Kampar land cover map. However the percentages we found after the classification showed the following distribution for this plot: 67.37% grassland and only 7.11% of the area was covered by oil palm. The bare soil class was not recognized in the Kampar land cover map (Fig. 18).

42

Figure 18. Overlay between HGU A, B, C, D, E and Kampar land cover map of 2014 with a scale 1:50,000. 43

5. Discussions 5.1

Landsat time series for idle land early detection

The land use history and existing land cover play an important role in idle land controlling especially in indicated idle land inventory. However, the National Land Agency of the Republic of Indonesia has many obstacles in this stage. Hence, it is necessary to develop methods in which they are able to overcome those obstacles in the monitoring system of land parcels all over Indonesia. Landsat images that are free and have an adequate spatial and high temporal resolution are possible to detect changes and able to give better contribution in land cover mapping (DeVries et al., 2015; Jia et al., 2014; Zhu and Woodcock, 2014). 5.1.1 Observing NDVI time series for land use history NDVI as one of vegetation indices provided by Landsat has been used for many purposes since it has strong relationship with vegetation (Pettorelli et al., 2005). An example of this relationship is that high NDVI values correspond to dense green vegetation (Bhalli et al., 2013; Myneni et al., 1997). Dutrieux et al. (2014) gives an example of an annual median NDVI time series representing a dynamic area in which dramatic forest changes occurred. In this thesis we produced annual median NDVI time series and analyzed the NDVI value to obtain important information related to land use history. Figure 12 reveals that HGU D has been planted at the same year when it got the land certificate in May 2008. The rectangular pattern shows there is a plantation, which is possibly associated with oil palm as stated by CRISP (2001) and Gandharum (2010). Furthermore, by observing the relative high NDVI values, around 0.7 to 0.8, in the rectangular pattern and considering the characteristics of mature oil palm trees that are ready for harvesting 2 until 4 years after planting (Basiron, 2007), it might be concluded that the oil palm already existed before the land owner registered the HGU certificate. From legal perspective, based on the regulation of right to cultivate (Peraturan Pemerintah Republik Indonesia Nomor 40 Tahun 1996 Tentang Hak Guna Usaha, Hak Guna Bangunan dan Hak Pakai Atas Tanah - Government Regulation No. 40 Year 1996), as long as the land owner has the forest conversion certificate from Ministry of Forestry, this plantation is legal. Otherwise, the person who plants this area is indicated as breaking the law by performing a plantation on forest land. The annual dynamic median NDVI value in some area in Figure 12 also might be considered there is an activity in this land parcel. Meanwhile, Figure 23 and Figure 24 (Appendix I) shows that both HGU B and C which were published in March 2005 and June 2007 had relative high stable NDVI values from the year 2008 until 2014, with values ranging from 0.6 to 0.8 for this period of seven years. We might assume that 10 years after the publication of land certificate, both land parcels have quite dense green vegetation in which can be assumed as mature oil palm, forest, or another mature plantations. Meanwhile, HGU A and HGU E (Appendix I Fig. 22 and 25), which were published in August and November 2003 respectively, had also a dynamic NDVI value. We are able to assume that there are plantations activities, those activities involve tree cut, deforestation or growth on both land parcels.

44

5.1.2 BFAST monitor in detecting land clearance In this study, we present a land clearance detection approach by the BFAST monitor for the five HGU land parcels. We combined breakpoints and negative magnitudes to determine the land clearance. We set five different magnitude thresholds for the magnitude value. The importance to set the magnitude threshold in combination with the breakpoints is: 1) the breakpoints do sometimes have positives magnitudes or magnitude values near zero, which do not represent deforestation (DeVries et al., 2015), and 2) to get the best combination of breakpoints and magnitude thresholds to determine land clearance; pixels in which land cover will change from forest to bare soil. We had different breakpoints in combination with different best thresholds when determining the land clearance for the five HGU land parcels. The best accuracies for HGU A, B, C, D, and E were achieved with breakpoints and a magnitude of -0.075, -0.070, -0.080, -0.060, and -0.075 respectively (Fig. 15). From those five different best thresholds for each study area, in general the magnitude threshold of -0.075, as the mode of the thresholds, is the best threshold to be used for detecting land clearance in this particular study area. The best threshold value for determining the land clearance class in this study was different than the breakpoint and a magnitude threshold values determined by DeVries et al. (2015). In his study, the breakpoint and magnitude combination led to a value of -0.18 for determining discrete forest land clearance in southern Ethiopia. In general, the changes in our study areas typically had a higher magnitude value than the changes in Ethiopia. It means the decrease in NDVI values in our study areas were less than in Ethiopia. Although both Riau-Sumatra and Ethiopia have a humid tropical climate and have a humid tropical rain forest (Margono et al., 2012; Schmitt et al., 2010), the location of the study areas for this research is closer to the equator line, which might affect the vegetation characteristics, e.g. vegetation height, canopy crown, etc. Hence, the magnitude of land clearance might be different between our study and the study performed by DeVries et al. (2015). 5.1.3 Potential errors and limitations from BFAST monitor analysis Land clearance detection is important in land parcel monitoring by NLA. The purpose is to check whether the land owner already uses the land as mentioned in the right entitlement or not. We have detected land clearance using the BFAST monitor algorithm by combining breakpoints and magnitude value. However, we still have potential errors and limitations in the present study, which will be discussed below. 5.1.3.1 Historical and monitoring data availability The history and monitoring period are important for near real time monitoring (Verbesselt et al., 2012). However, it was a challenge to get dense temporal Landsat images for our time series because our five HGU land parcels are located in tropical forest area and cloud cover is a main issue for tropical areas (Hansen et al., 2009; Margono et al., 2012; Mitchard et al., 2011). Furthermore, gaps due to Landsat 7 ETM+ Scan line corrector (SLC-off) might be the source of data losses in our study areas (Fig. 19). Table 4 shows us that during the rainy season, especially in January and December, the number of Landsat images is very low. The data availability in the history and monitoring period might cause errors in our study. The history period is an important aspect for change detection in near real-time, which is important to 45

define a stable reference model (DeVries et al., 2015; Verbesselt et al., 2012; Zhu and Woodcock, 2014). In the five years period, from 1st of January 2008 until the 31st of December 2012, the mean clear sky observations for the five selected study areas ranged from 14.78 to 25.32 out of the 81 available Landsat scenes (Table 9). It means that the history data availability for the five study areas was quite low, in average; there were 55 to 66 loss observations during that period. To determine our history period, we followed the paper of Verbesselt et al. (2012) which states that

the length of the history period should be more than 2 years (>2years). However, during the data analysis in this study we found out that a history period of more than 2 years was insufficient to obtain a good result in land clearance detection. We observed two pixels on HGU A and found that both pixels had large no data gaps for the period between the middle of 2011 until the end of 2012, which was resulting in the false breakpoints (Fig. 19). On the other hand, an example in which the fewer gaps in the history period resulted in true breakpoints is shown for some pixels in Fig. 14 and Appendix V Fig. 37. The observation data gaps might be caused by the high cloud cover and the fact that SLC was off (Fig.19). With the area of 23.318 ha, it was very challenging to get full observation due to the Landsat striping. Moreover, Fig. 19 shows the cloud cover for observation and not available data per pixel. It is clear that in the period of 2008, 2010, 2011, and 2012, the cloud cover was relatively high, with values above 50% and some scenes were almost completely covered by clouds, 90% of cloud cover. Those images with a high amount of cloud cover were especially found at the beginning and the end of the year, due to the rain season. Our low temporal frequency of time series data for the history period is in line with what is stated in the paper of DeVries et al. (2015). This low history data leads to a problem when the time series method forms a stable model (DeVries et al., 2015) and might increase the commission errors (Zhu and Woodcock, 2014), which means that the algorithm detects changes at locations in which in reality there is no changes. To avoid these commission errors, we should use as many observations as possible (Zhu and Woodcock, 2014). The data availability during our monitoring period might be the cause of our potential errors. The data during our monitoring period represents the new observed data that has been monitored for existence of disturbance. We used 26 Landsat images in our monitoring period, which was the period between the 1st of January 2013 and the 1st of July 2014. The mean clear sky observations for the five selected study areas ranged from 5.68 to 6.37 out of the 26 available Landsat scenes (Table 9). It means our monitoring data availability for our study areas was quite low, in average; there were around 20 loss observations during this period. As the BFAST monitor works only for every available pixel in the time series, the low temporal data in the monitoring period might lead to high omission errors, which will result to the fact that the BFAST monitor does not detect the real changes.

46

A

2014 2014

B

C

D

Figure 19. Example of the commission errors arising from low data availability in the temporal history period. Both pixels (red point (A)) had large gaps of no data in the period from the middle of 2011 until the end of 2012, which led to false break detection (B). The Landsat 3-2-1 composites show the four periods from 2011 until 2012 with cloud cover and SLC off (C). The last graph (D) shows the percentage of cloud cover per observation and the not available (gaps) data for both pixels. 47

5.1.3.2 Forest mask determination The possible error in BFAST monitor used in this study was the determination of the forest mask at the beginning of the monitoring period. Our forest mask was only determined by supervised classification of SPOT 6 images for the year 2013 (Appendix II Fig. 26 and 27). We did not refine our forest mask with a tree cover from Moderate Resolution Imaging Spectroradiometer Vegetation Continuous Fields (MODIS VCF) product as DeVries et al. (2015). The misclassification might occur especially when the HGU land parcels were already planted with oil palm (Appendix II Fig. 26). Hence, the omission error might increase in which the BFAST monitor would not detect the land clearance of misclassified pixels. 5.1.3.3 Temporal accuracy assessment This study has a limitation in change detection accuracy in the temporal domain which is in line with Verbesselt et al. (2012). It means we did not assess the accuracy of the changes time detected by BFAST monitor because of a lack of information when the land clearance took place and limited amount of SPOT 6 images to validate the changes. We only used two SPOT 6 images which were acquired on 7 January 2013 and 17 June 2014. Hence, we can only observe the changes during that period and not in detailed detected breakpoints within this time interval. However, our results confirmed all changes took place in the period between January and June 2014, as our starting and end of the monitoring time (Appendix III Fig. 28, 29, 30, 31, and 32), from which can be assessed that all the changes took place during the monitoring time. 5.1.3.4 The cause of land clearance One limitation of the BFAST monitor is that it does not directly provide information on the causes of forest disturbance (Verbesselt et al., 2012). Hence, we cannot make any analysis what is the source of the land clearance in the five HGU land parcels, whether it is because a proper land clearance or illegal land clearance, e.g., forest fire. Many farmers in Riau-Sumatra use fire in land clearance for oil palm plantation (Suyanto, et al., 2004) which is against the Government Regulation No. 11/2001 article 11 that prohibit people from engaging in forest fire activities.

5.2

Landsat images classification to support for idle land controlling

Relatively lower values in the producer’s, user’s accuracy, and also largest confusions have been found in the classification of the land cover classes oil palm, forest, and grassland. The largest confusions were between oil palm and grassland, oil palm and forest, and also between forest and grassland (Table 14). As can be seen in Fig. 20, the spectral signatures of oil palm, forest, and grassland in some area are overlapping. Even in the layer 4 as the representation of band 5 (NIR) which is suitable for vegetation detection (Jia et al, 2014; USGS, 2013b), these three land covers are not well distinguished. Our result is in line with Jia et al. (2014) in which forest and grassland has the large confusion. This might occur, due to the fact that grasslands are often present close to the boundary of forest or in some small areas even in the middle of the forest (Jia et al., 2014). Moreover, in our case, grassland might be present near oil palm and oil palm was sometimes adjoined with forest. These conditions led to some misclassifications. 48

Figure 20. Spectral signature plot of land cover classes in Landsat 8 classification.

The accuracy assessment showed that the bare soil class had the highest accuracy among the other land covers (Table 14). This is in accordance with the result of Jia et al. (2014), in which the bare soil land cover class also had the highest accuracy among the vegetation classes. The overall accuracy of this Landsat land cover classification was relative high with a percentage of 86.13% and a Kappa coefficient of 0.78 (a substantial agreement) (Table 14). Several comparable studies in which Landsat supervised maximum likelihood land cover classification has been performed resulted in high overall accuracies as well, like the studies of Yuan et al. (2005), Wu et al. (2006), and Jia et al. (2014), which have an overall accuracy of 93.98%, 89.32% and 90.4% respectively. Our approach is in line with Jia et al. (2014), in which Landsat 8 supervised maximum likelihood classification of surface reflectance bands 2 until 7 has been performed. Therefore, we can conclude that Landsat 8 OLI data has a satisfactory performance in land cover classification. Apparently, although Landsat 8 maximum likelihood supervised classification gave relative high total accuracies and Kappa coefficient, some classes had a low accuracy, which resulted in a large confusion (Table 14). Hence, it is quite difficult to assess the accuracy of those classes and distinguish between the different land covers. Bare soil seems to be the only land cover that might give a high accuracy, which can be used to detect idle land. The possible sources of errors in our classification are: 1) The determination of the reference classes for training and accuracy assessment. Our determination of reference classes was only based on visual interpretation of SPOT 6 images. We might have made some mistakes when determining the training areas and validating the classification result. Hence, it might cause a lower accuracy for some land cover classes. 2) The determination of the training areas. The training areas used for our classification were based on visual interpretation of our SPOT image. Although our method in training 49

areas determination follows the approach used by Yuan et al. (2005), there is a possibility that our training areas were not evenly distributed. Hence, again it might cause a lower accuracy for some land cover classes. 3) The quality of Kampar land cover map as the validation data. According to the Geospatial Information Agency of the Republic of Indonesia, in general land cover maps are made based on aerial photographs using ground control points of the Control Network Center of Geodesy and Geodynamics for correction purpose. The level of detail of every single map is adjusted to the needs of the users in which every map scale has its own purposes and refers to Land Cover Classification System United Nation – Food and Agriculture Organization (LCCS-UNFAO) and ISO 19144-1 (BSN, 2010). Moreover, land cover maps with a scale of 1:50,000 have generalized land cover classes and are not very detailed (BSN, 2010). This lack in detail can be shown by the fact that the area of 1 cm2 on the map corresponds with 500 x 500 m = 250,000 m2 = 25 ha. HGU plot A and plot E have in total a size of around 25 ha, which is only 1 cm2 on the map. Besides this lack in detail we were also missing information about the classification techniques or inventory techniques used to determine the Kampar land use classes. An example of how the classification technique or inventory technique influences the classification of the landscape units is discussed in the paper of Hubert-Moy et al. (2001). In this paper a description is given of how the choice of a classification technique can significantly influence the results of the classification and what the effect is on the landscape units classified or identified within a study area. This in combination with the fact that we did not have any information about the quality of our land cover map, makes it important that we keep into consideration that the identified land use classes in the Kampar land use map also has some inaccuracies, which influence the accuracy assessment.

5.3

The combination of NDVI time series profile, BFAST monitor, and Landsat classification as idle land early detection

Idle land controlling is one of the priority programs in the National Land Agency of the Republic of Indonesia which aims to make land right holders use their land according to the purpose of the entitlement. Since this program will abolish the legal relationship between the land owner and their land, every stage in idle controlling should be carried out based on the existing regulations. As mentioned in the explanation of idle land inventory stage (Fig. 3), the NLA officers should put attention to existing land cover and its land use history, for instance if officers discover bare soil or grassland as the existing land cover, they have to check whether it is already bare soil or grassland for a while or it is just there for a short period of time. In this study, we combine the results of the BFAST monitor in detecting land clearance and Landsat supervised maximum likelihood in detecting existing land cover to support idle land controlling in Indonesia. Figure 21 illustrates an example of how we are able to obtain some findings which can be taken into consideration during indicated idle land inventory. It shows the combination of the BFAST monitor and Landsat 8 classification as idle land early detection (showed in red circle): (a) Land clearance detected by BFAST monitor (SPOT image in January 2013), (b) SPOT 6 image in June 2014, (c) Landsat surface reflectance bands in 5-4-3 composites in June 2014, (d) land cover map resulted from Landsat 8 classification. All these 50

data sources are available for the NLA. In particular after the newest development that the National Institute of Aeronautics and Space provides high resolution satellite images, including SPOT images, to 9 departments in Indonesia including NLA from January 2015 onwards. However, the temporal resolution of SPOT images is lower than Landsat images which are freely downloadable. If we look at Landsat 8 image in 5-4-3 composites (Fig. 21 (c)), we cannot recognize what the land covers are in HGU E. With the infrared false color, we only might recognize some parts of the vegetation (red) and some parts of non vegetation (non red). After we classify the Landsat image, we get the land cover map (Fig. 21 (d)). From the confusion matrix in Table 14, we can conclude that bare soil, as early indication of idle land, had the highest accuracy compared to other land cover classes. Red circle in Fig. 21 (d) shows that there is bare soil as the existing land cover. However, we cannot just make final conclusion based on this land cover, we still should look at what might happen in the past. Hence, we look at our BFAST monitor result to check what might happen in the past. Red circle in Fig. 21 (a) shows that for the period January 2013 and June 2014, there is a land clearance/activity on HGU E which means the land owner uses the land. Based on the BFAST monitor result, this area is eight pixels which equals to 0.72 ha.

Magnitude

07-01-2013

17-06-2014

(a)

(b)

12-06-2014 12-06-2014

(c) (c)

(d)

Figure 21. The combination of BFAST monitor and Landsat classification as idle land early detection. 51

Next, we look at grassland, which is indicated with the black arrow (Fig. 21 (d)). The existence of grassland in the existing land cover is also sometimes indicated as idle land. During 2013 until 2014, we do not have any changes, because this area is not considered as forest class. Therefore the BFAST monitor algorithm does not take into account this area in the calculation. Hence, we look at our annual median NDVI to see what might happen in this area. Fig. 25 in Appendix I shows there is a decrease in NDVI values during 2012 and 2013, which might lead to the conclusion that there is an activity in that area (black arrow). These observations and findings would help NLA to decide whether HGU E will be targeted for idle controlling or not.

5.4

Future research needs

Although the BFAST monitor in our study was able to detect land clearance in our land parcels, some limitations in our study should be considered as the input for future research. This study only used five HGU land parcels. This is, due to limitations in available high resolution images for validation. Further improvements of our methodology should include a higher number of HGU land parcels, equally spread over all provinces in Indonesia. Hence, we would be able to assess the sensitivity of BFAST monitor in detecting land clearance with various conditions and circumstances. The low data availability in our study areas was the major problem which led to high commission and omission errors. According to Zhu and Woodcock (2014), as many observations as possible is needed in change detection. Hence, we should use a longer time in our history period in future research. The low data availability might be overcome by combining the time series approach with another remote sensing method, e.g., data integration of multi-temporal Landsat and Phased Array type L-band Synthetic Aperture Radar in forest monitoring (Eberenz, 2014; Morel, et al., 2012). This integration is a topic for future research. Furthermore, because “NDVI is known to be less sensitive to the forest structure than other spectral indices” (DeVries et al., 2015), it will be more interesting if we used other vegetation indices provided by Landsat in detecting forest disturbances/land clearance, e.g. Normalized Difference Moisture Index (Hais, et al., 2009; Jin and Sader, 2005; Wilson and Sader, 2002) and Soil Adjusted Vegetation Index (Matricardi et al., 2010). In addition, we might use medium-high resolution satellite images such as Landsat 8 (15 and 30 m spatial resolution) and SPOT 6 (1.5 and 6 m spatial resolution) in order to better detect land clearance. Another further improvement of our methodology should include multiple classifier approaches for pixel based classification, not only maximum likelihood, but also approaches like e.g. minimum distance classifier (Hubert et al., 2001) or support vector machine (Jia et al., 2014; Pal and Foody, 2012). We also might use object based classification as a non-pixel based classification to detect existing land cover classes in HGU land parcels (Desclée et al., 2006; Im et al., 2008; Stow et al, 2007). Government Regulation No.40/1996 article 12 (1-f) states that every HGU rights owner should provide a written report, by the end of every year related to their own specific HGU land parcels. Since idle land controlling is one of priority programs in NLA, it might be a further improvement for NLA to build a system in which HGU land owners have to publish online reports regularly, 52

about the progress on their lands, including the land clearance. Hence, NLA is able to use the progress reports to support idle land controlling and other programs related to land management in Indonesia.

53

6. Conclusions In this section, an overview of conclusions about the main findings related to idle land early detection is provided. It also addresses the results for each research questions formulated at the beginning of this study and concludes about the feasibility of BFAST monitor as a method for idle land early detection. Furthermore, a conclusion will be made about suitability and accuracy of Landsat maximum likelihood supervised classification as an approach for idle land early detection. For the National Land Agency of Indonesia, the land use history and the existing land cover plays an important role in the idle land controlling, especially for the indicated idle land inventory. In our study, it was shown that the annual NDVI time series and the BFAST monitor were able to show the land use history of our five selected HGU land parcels. The annual median NDVI was able to show the plantation pattern and from its value, we could derive the possible activities on the HGU land parcels. Furthermore, land clearance practices were evidenced from the Landsat time series using the BFAST monitor by combining the breakpoints and different negative magnitude threshold values. Overall, the magnitude threshold between -0.060 and 0.080 was the best threshold in detecting land clearance for each study area. However, in general the magnitude threshold of -0.075, as the mode of the thresholds, was the best threshold to be used in land clearance detection for this particular study area. For assessing existing land cover on the HGU land parcels, Landsat 8 supervised maximum likelihood classification of surface reflectance bands 2 until 7 has been used. This approach had a good performance in land cover classification, with an overall accuracy of 86.13% and a Kappa coefficient of 0.78, which shows a substantial agreement. We concluded that detecting bare soil, as one of the indications of idle land, is more accurate than vegetation i.e. oil palm, forest, and grassland. A major aspect which limited the usability of our methodology was the high cloud cover in tropical areas. This situation led to a low data availability, which results in high commission and omission errors. By bringing together all our results in this study, it gives many opportunities for NLA to have an idle land early detection method to decide whether a land parcel will be targeted in the idle land controlling or not. A continuously development in the idle land controlling stages is necessary for NLA to gain a good land management in Indonesia.

54

References Badan Pertanahan Nasional. (2010). Penertiban Dan Pendayagunaan Tanah Terlantar Berdasarkan Peraturan Pemerintah Nomor 11 Tahun 2010. A presentation in front of Ministry of Home Affairs of the Republic of Indonesia. Badan Pusat Statistik (National Statistic Agency of the Republic of Indonesia). (2014). Retrieved September 19, 2014, from http://www.bps.go.id/. Basiron, Y. (2007). Palm oil production through sustainable plantations. European Journal of Lipid Science and Technology, 109(4), 289-295. Bell, K.C. & Srinivas, S. (2013). A BIG Change : land management Indonesia. Geospatial World, June 2013, 38 – 41. Bhalli, M. N., Ghaffar, A., Shirazi, S. A., & Parveen, N. (2013). An analysis of the Normalized Difference Vegetation Index (NDVI) and its relationship with the population distribution of Faisalabad-Pakistan. Pakistan Journal of Science, 65(4). BSN - Badan Standarisasi Nasional (National Standardization Agency of the Republic of Indonesia). (2010). SNI: klasifikasi penutup lahan. Retrieved on March 2, 2015 from http://www.bakosurtanal.go.id/assets/download/sni/SNI/15.%20SNI%207645-2010% 20Klasifikasi% 20penutup% 20lahan.pdf Chander, G., Markham, B. L., & Helder, D. L. (2009). Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sensing of Environment, 113(5), 893-903. Congalton, R. G. (2001). Accuracy assessment and validation of remotely sensed and other spatial information. International Journal of Wild land Fire, 10(4), 321-328. CRISP. (2001). The Centre for Remote Imaging, Sensing and Processing: interpreting optical remote sensing images. Retrieved February 23, 2015 from http://www.crisp.nus.edu.sg/~research/tutorial/opt_int.htm. Dephut. (2014). Profil kehutanan 33 provinsi: Provinsi Riau. Retrieved in December 3, 2014 from http://www.dephut.go.id/uploads/files/76333af5b0c4474a6498f7d3d1303470.pdf. Desclée, B., Bogaert, P., & Defourny, P. (2006). Forest change detection by statistical objectbased method. Remote Sensing of Environment, 102(1), 1-11. DeVries, B., Verbesselt, J., Kooistra, L., & Herold, M. (2015). Robust monitoring of small-scale forest disturbances in a tropical montane forest using Landsat time series. Remote Sensing of Environment. Donald, P. F. (2004). Biodiversity impacts of some agricultural commodity production systems. Conservation Biology, 18(1), 17-38. Dutrieux, L., DeVries, B., & Verbesselt, J. (2014). Introduction to bfastSpatial. Retrieved from http://dutri001.github.io/bfastSpatial/. Eberenz, J. (2014). Integration of Landsat and SAR time series for near real-time deforestation monitoring (Master thesis). Wageningen University, Wageningen. Feintrenie, L., Chong, W. K., & Levang, P. (2010). Why do farmers prefer oil palm? Lessons learnt from Bungo district, Indonesia. Small-Scale Forestry, 9(3), 379-396. Foody, G. M. (2002). Status of land cover classification accuracy assessment. Remote Sensing of Environment, 80(1), 185-201. Gandharum, L. (2010). Classification of Oil Palm in Indonesia Using FORMOSAT-2 Satellite Image (Doctoral dissertation, Master thesis). National Central University (NCU), Taiwan. 55

Germer, J., & Sauerborn, J. (2008). Estimation of the impact of oil palm plantation establishment on greenhouse gas balance. Environment, Development and Sustainability, 10(6), 697716. Green, K. (2000). Selecting and interpreting high-resolution images. Journal of Forestry, 98(6), 37-40. Hansen, M. C., Stehman, S. V., Potapov, P. V., Arunarwati, B., Stolle, F., & Pittman, K. (2009). Quantifying changes in the rates of forest clearing in Indonesia from 1990 to 2005 using remotely sensed data sets. Environmental Research Letters, 4(3), 034001. Hashim, U. M., & Ahmad, A. (2014). The effects of training set size on the accuracy of maximum likelihood, neural network and support vector machine classification. Science International, 26(4), 1477-1481. Hais, M., Jonášová, M., Langhammer, J., & Kučera, T. (2009). Comparison of two types of forest disturbance using multitemporal Landsat TM/ETM+ imagery and field vegetation data. Remote Sensing of Environment, 113(4), 835-845. Herold, M., Woodcock, C. E., Di Gregorio, A., Mayaux, P., Belward, A. S., Latham, J., & Schmullius, C. C. (2006). A joint initiative for harmonization and validation of land cover datasets. Geoscience and Remote Sensing, IEEE Transactions, 44(7), 1719-1727. Holmes, D. (2000). Deforestation in Indonesia: a review of the situation in Sumatra, Kalimantan and Sulawesi. Jakarta: World Bank (processed). Hopwood, D. (2014). Status of Astrium GEO-Information services, EO satellite constellation. Retrieved February 7, 2015 from https://calval.cr.usgs.gov/wordpress/wp-content/ uploads/Astrium-Constellation-Status.pdf. Huang, C., Davis, L. S., & Townshend, J. R. G. (2002). An assessment of support vector machines for land cover classification. International Journal of Remote Sensing, 23(4), 725-749. Huang, C., Goward, S. N., Masek, J. G., Gao, F., Vermote, E. F., Thomas, N., Townshend, J. R. et al. (2009). Development of time series stacks of Landsat images for reconstructing forest disturbance history. International Journal of Digital Earth, 2(3), 195-218. Huang, C., Thomas, N., Goward, S. N., Masek, J. G., Zhu, Z., Townshend, J. R., & Vogelmann, J. E. (2010). Automated masking of cloud and cloud shadow for forest change analysis using Landsat images. International Journal of Remote Sensing, 31(20), 5449-5464. Hubert-Moy, L., Cotonnec, A., Le Du, L., Chardin, A., & Perez, P. (2001). A comparison of parametric classification procedures of remotely sensed data applied on different landscape units. Remote Sensing of Environment, 75(2), 174-187. Hurd, J. D., & Civco, D. L. (2009). Creating an image dataset to meet your classification needs: a proof-of-concept study. In Baltimore, Maryland: ASPRS 2009 Annual Conference. Hutchinson, J. M. S., Jacquin, A., Hutchinson, S. L., & Verbesselt, J. (2015). Monitoring vegetation change and dynamics on US Army training lands using satellite image time series analysis. Journal of Environmental Management, 150, 355-366. Ibrahim, S., Hasan, Z. A., & Khalid, M. (2001). Application of optical remote sensing technology for oil palm management. In Cutting-edge technologies for sustained competitiveness: Proceedings of the 2001 PIPOC International Palm Oil Congress, Agriculture Conference, Kuala Lumpur, Malaysia, 20-22 August 2001. (499-504). Malaysian Palm Oil Board (MPOB). Im, J., Jensen, J. R., & Hodgson, M. E. (2008). Object-based land cover classification using high-posting-density LiDAR data. GIScience & Remote Sensing, 45(2), 209-228. 56

Jia, K., Wei, X., Gu, X., Yao, Y., Xie, X., & Li, B. (2014). Land cover classification using Landsat 8 Operational Land Imager data in Beijing, China. Geocarto International, 29(8), 941-951. Jin, S., & Sader, S. A. (2005). Comparison of time series tasseled cap wetness and the normalized difference moisture index in detecting forest disturbances. Remote Sensing of Environment, 94(3), 364-372. Jusoff, K. (2009). Sustainable management of a matured oil palm plantation in UPM campus, Malaysia using airborne remote sensing. Journal of Sustainable Development, 2(3), 195. Kementerian Dalam Negeri (Ministry of Home Affairs of the Republic of Indonesia). (2014a). Retrieved September 19, 2014, from http://otda.kemendagri.go.id/images/file/data2014/ file_konten/jumlah_daerah_otonom_ri.pdf Kementerian Dalam Negeri (Ministry of Home Affairs of the Republic of Indonesia). (2014b). Profil Provinsi Riau. Retrieved February 10, 2015, from http://www.kemendagri.go.id/ pages/profil-daerah/provinsi/detail/14/riau Lambert, J., Jacquin, A., Denux, J.-P., & Chéret, V. (2011). Comparison of two remote sensing time series analysis methods for monitoring forest decline. In Analysis of Multi-temporal Remote Sensing Images (Multi-Temp), 2011 6th International Workshop. IEEE, 93-96. Landis, J. R., & Koch, G. G. (1977). The measurement of observer agreement for categorical data. Biometrics, 159-174. Lillesand, T. M., Kiefer, R. W., & Chipman, J. W. (2008). Remote sensing and image interpretation (No. Ed. 6). Hoboken: John Wiley & Sons Ltd. Lu, D., Mausel, P., Brondizio, E., & Moran, E. (2004). Change detection techniques. International Journal of Remote Sensing, 25(12), 2365-2401. Margono, B. A., Turubanova, S., Zhuravleva, I., Potapov, P., Tyukavina, A., Baccini, A., Hansen, M. C., et al. (2012). Mapping and monitoring deforestation and forest degradation in Sumatra (Indonesia) using Landsat time series data sets from 1990 to 2010. Environmental Research Letters, 7(3), 034010. Masek, J. G., Honzak, M., Goward, S. N., Liu, P., & Pak, E. (2001). Landsat-7 ETM+ as an observatory for land cover: Initial radiometric and geometric comparisons with Landsat-5 Thematic Mapper. Remote Sensing of Environment, 78(1), 118-130. Matricardi, E. A., Skole, D. L., Pedlowski, M. A., Chomentowski, W., & Fernandes, L. C. (2010). Assessment of tropical forest degradation by selective logging and fire using Landsat imagery. Remote Sensing of Environment, 114(5), 1117-1129. Maxwell, S. K., Airola, M., & Nuckols, J. R. (2010). Using Landsat satellite data to support pesticide exposure assessment in California. International Journal of Health Geographics, 9, 46-46. McCarthy, J. F., & Cramb, R. A. (2009). Policy narratives, landholder engagement, and oil palm expansion on the Malaysian and Indonesian frontiers. The Geographical Journal, 175(2), 112-123. Mitchard, E. T. A., Saatchi, S. S., White, L. J. T., Abernethy, K. A., Jeffery, K. J., Lewis, S. L., et al. (2011). Mapping tropical forest biomass with radar and spaceborne LiDAR: overcoming problems of high biomass and persistent cloud. Biogeosciences Discussions, 8(4), 8781. Mongabay. (2012). In pictures: Rainforests to palm oil. Retrieved February 23, 2015 from http://news.mongabay.com/2012/0702-rainforests-to-palm-oil-photos.html.

57

Morel, A. C., Fisher, J. B., & Malhi, Y. (2012). Evaluating the potential to monitor above ground biomass in forest and oil palm in Sabah, Malaysia, for 2000–2008 with Landsat ETM+ and ALOS-PALSAR. International Journal of Remote Sensing, 33(11), 3614-3639. Myneni, R. B., Keeling, C. D., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). Increased plant growth in the northern high latitudes from 1981 to 1991. Nature, 386(6626), 698-702. Olofsson, P., Foody, G. M., Herold, M., Stehman, S. V., Woodcock, C. E., & Wulder, M. A. (2014). Good practices for estimating area and assessing accuracy of land change. Remote Sensing of Environment, 148, 42-57. Olson, C. E. (1960). Elements of photographic interpretation common to several sensors. Photogrammetric Engineering, 26(4), 651-656. Pagiola, S. (2000). Land use change in Indonesia. Background paper prepared for the Environment Department, World Bank, Washington, DC. Pal, M., & Foody, G. M. (2012). Evaluation of SVM, RVM and SMLR for accurate image classification with limited ground data. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 5(5), 1344-1355. Partoyo & Shrestha, R. P. (2013). Monitoring farmland loss and projecting the future land use of an urbanized watershed in Yogyakarta, Indonesia. Journal of Land Use Science, 8(1), 5984. Pettorelli, N., Vik, J. O., Mysterud, A., Gaillard, J. M., Tucker, C. J., & Stenseth, N. C. (2005). Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends in Ecology and Evolution, 20(9), 503-510. Pontius, R. G. (2000). Quantification error versus location error in comparison of categorical maps. Photogrammetric Engineering and Remote Sensing, 66(8), 1011-1016. Ramdani, F., & Hino, M. (2013). Land use changes and GHG emissions from tropical forest conversion by oil palm plantations in Riau Province, Indonesia.PloS one, 8(7), e70323. Roy, D. P., Wulder, M. A., Loveland, T. R., Woodcock, C. E., Allen, R. G., Anderson, M. C., Zhu, Z., et al. (2014). Landsat-8: Science and product vision for terrestrial global change research. Remote Sensing of Environment, 145, 154-172. Santoso, H., Gunawan, T., Jatmiko, R. H., Darmosarkoro, W., & Minasny, B. (2011). Mapping and identifying basal stem rot disease in oil palms in North Sumatra with QuickBird imagery. Precision Agriculture, 12(2), 233-248. Schmitt, C. B., Denich, M., Demissew, S., Friis, I., & Boehmer, H. J. (2010). Floristic diversity in fragmented Afromontane rainforests: altitudinal variation and conservation importance. Applied Vegetation Science, 13(3), 291-304. Shaker, A., Yan, W. Y., & El-Ashmawy, N. (2012). Panchromatic Satellite Image Classification for Flood Hazard Assessment. Journal of Applied Research and Technology, 10(6), 902911. Sheil, D., Casson, A., Meijaard, E., Van Noordwijk, M., Gaskell, J., Sunderland-Groves, J., Kanninen, M., et al. (2009). The impacts and opportunities of oil palm in Southeast Asia: What do we know and what do we need to know?. Center for International Forestry Research. SIC. (2014). Satellite Imaging Corporation: Google Earth vs. Custom Satellite Imagery. Retrieved May 18, 2015 from http://www.satimagingcorp.com/contact/google_earth/. Song, C., Woodcock, C. E., Seto, K. C., Lenney, M. P., & Macomber, S. A. (2001). Classification and change detection using Landsat TM data: when and how to correct atmospheric effects?. Remote sensing of Environment, 75(2), 230-244. 58

Stehman, S. V. (1996). Estimating the kappa coefficient and its variance under stratified random sampling. Photogrammetric Engineering and Remote Sensing, 62(4), 401-407. Stehman, S. V. (1997). Selecting and interpreting measures of thematic classification accuracy. Remote sensing of Environment, 62(1), 77-89. Stehman, S. V., Wickham, J. D., Smith, J. H., & Yang, L. (2003). Thematic accuracy of the 1992 National Land-Cover Data for the eastern United States: Statistical methodology and regional results. Remote Sensing of Environment, 86(4), 500-516. Stehman, S. V., Olofsson, P., Woodcock, C. E., Herold, M., & Friedl, M. A. (2012). A global land-cover validation data set, II: augmenting a stratified sampling design to estimate accuracy by region and land-cover class. International Journal of Remote Sensing, 33(22), 6975-6993. Stehman, S. V., & Czaplewski, R. L. (1998). Design and analysis for thematic map accuracy assessment: fundamental principles. Remote Sensing of Environment, 64(3), 331-344. Stow, D., Lopez, A., Lippitt, C., Hinton, S., & Weeks, J. (2007). Object‐based classification of residential land use within Accra, Ghana based on QuickBird satellite data. International Journal of Remote Sensing, 28(22), 5167-5173. Suara, P. (2007). Izin 28 Perusahaan Kebun Sawit di Jambi Akan Dicabut. [Partnership or power relationship] Studi Kemitraan di PT Sari Aditya Loka/Astra Group di Provinsi Jambi, Setara Jambi/Yayasan Keadilan Rakyat , p. 9. Sunaryathy, P. I., Busu, I., & Rasib, A. W. (2010). Evaluation of net primary productivity of oil palm plantation in South Sulawesi Indonesia. In 33rd Asian Conference on Remote Sensing 2012. ACRS 2012, (2), 1830-1834. Suyanto, S., Applegate, G., Permana, R. P., Khususiyah, N., & Kurniawan, I. (2004). The role of fire in changing land use and livelihoods in Riau-Sumatra. Ecology and Society, 9(1), 15. Thomlinson, J. R., Bolstad, P. V., & Cohen, W. B. (1999). Coordinating methodologies for scaling landcover classifications from site-specific to global: Steps toward validating global map products. Remote Sensing of Environment, 70(1), 16-28. Turner, P. D. and R. Gillbanks (1974). Oil palm cultivation and management. Kuala Lumpur Malaysia Incorporated Society of Planters. USGS. (2008). Opening the Landsat Archive: Imagery for everyone. Retrieved September 24, 2014 from http://landsat.usgs.gov/documents/USGS_Landsat_Imagery_Release.pdf. USGS. (2013). The best spectral bands for some studies. Retrieved February 27, 2015 from http://landsat.usgs.gov/best_spectral_bands_to_use.php. USGS. (2014). The band designations for the Landsat satellites. Retrieved February 7, 2015 from http://landsat.usgs.gov/band_designations_landsat_satellites.php. USGS. (2015a). Product Guide Provisional Landsat 8 Surface Reflectance Product. Retrieved February 15, 2015 from http://landsat.usgs.gov/ documents/provisional_l8sr_ product_guide.pdf. USGS. (2015b). Provisional Landsat 8 Surface Reflectance Data Available. Retrieved February 15, 2015 from http://landsat.usgs.gov/about_LU_Special_Issue_3.php. Verbesselt, J., Hyndman, R., Newnham, G., & Culvenor, D. (2010a). Detecting trend and seasonal changes in satellite image time series. Remote Sensing of Environment, 114(1), 106-115. Verbesselt, J., Hyndman, R., Zeileis, A., & Culvenor, D. (2010b). Phenological change detection while accounting for abrupt and gradual trends in satellite image time series. Remote Sensing of Environment, 114, 2970-2980. 59

Verbesselt, J., Zeileis, A., & Herold, M. (2012). Near real-time disturbance detection using satellite image time series. Remote Sensing of Environment, 123, 98-108. Vicente, S. M., Pérez-Cabello, F., & Lasanta, T. (2008). Assessment of radiometric correction techniques in analyzing vegetation variability and change using time series of Landsat images. Remote Sensing of Environment, 112 (10), 3916-3934. Weng, Q. (2010). Remote sensing and GIS integration: theories, methods, and applications (p. 416). New York: McGraw-Hill. Wilson, E. H., & Sader, S. A. (2002). Detection of forest harvest type using multiple dates of Landsat TM imagery. Remote Sensing of Environment, 80(3), 385-396. Winoto, J. (2009). Taking land policy and administration in Indonesia to the next stage and National Land Agency’s strategic plan. International Federation of Surveyors Forum. Wu, Q., Li, H., Q., Wang, R., S., Paulussen, J., He., Y., Wang, M., Wang, B., H., Wang, Z. (2006). Monitoring and predicting land use change in Beijing using remote sensing and GIS. Landscape and Urban Planning, 78(4): 322-333. Xian, G., Homer, C., & Fry, J. (2009). Updating the 2001 National Land Cover Database land cover classification to 2006 by using Landsat imagery change detection methods. Remote Sensing of Environment, 113(6), 1133-1147. Yang, X., & Lo, C. P. (2002). Using a time series of satellite imagery to detect land use and land cover changes in the Atlanta, Georgia metropolitan area. International Journal of Remote Sensing, 23(9), 1775-1798. Yuan, F., Sawaya, K. E., Loeffelholz, B. C., & Bauer, M. E. (2005). Land cover classification and change analysis of the Twin Cities (Minnesota) Metropolitan Area by multitemporal Landsat remote sensing. Remote Sensing of Environment, 98(2), 317-328. Zhang, Q., Wang, J., Gong, P., & Shi, P. (2003). Study of urban spatial patterns from SPOT panchromatic imagery using textural analysis. International Journal of Remote Sensing, 24(21), 4137-4160. Zhu, Z., & Woodcock, C. E. (2014). Continuous change detection and classification of land cover using all available Landsat data. Remote Sensing of Environment, 144, 152-171. List of Laws and Regulations: Undang-undang Dasar Republik Indonesia 1945 (the Constitution 1945). Undang-undang Nomor 5 Tahun 1960 Tentang Peraturan Dasar Pokok-pokok Agraria (the Basic Agrarian Law No. 5 Year 1960). Peraturan Pemerintah Republik Indonesia Nomor 40 Tahun 1996 Tentang Hak Guna Usaha, Hak Guna Bangunan dan Hak Pakai Atas Tanah (Government Regulation No. 40 Year 1996) Peraturan Pemerintah Republik Indonesia Nomor 4 Tahun 2001 Tentang Pengendalian Kerusakan dan atau Pencemaran Lingkungan Hidup yang Berkaitan dengan Kebakaran Hutan dan atau Lahan (Government Regulation No. 4 Year 2001). Peraturan Pemerintah Republik Indonesia Nomor 11 Tahun 2010 Tentang Penertiban dan Pendayagunaan Tanah Terlantar (Government Regulation No. 11 Year 2010).

60

Peraturan Kepala Badan Pertanahan Nasional Republik Indonesia Nomor 4 Tahun 2010 Tentang Tata Cara Penertiban Tanah Terlantar (Regulation of the Head Office of National Land Agency No. 4 Year 2010). Peraturan Kepala Badan Pertanahan Nasional Republik Indonesia Nomor 9 Tahun 2011 Tentang Perubahan Atas Peraturan Kepala Badan Pertanahan Nasional RI Nomor 4 Tahun 2010 Tentang Tata Cara Penertiban Tanah Terlantar (Regulation of the Head Office of National Land Agency No. 9 Year 2011).

61

Appendixes Appendix I. The annual median NDVI of HGU A, B, C, and E for period 2008 until 2014.

0.8 0.7 0.6 0.5

0.8 0.7 0.6 0.5

0.8 0.75 0.7 0.65 0.6 0.55

0.9 0.8 0.7 0.6 0.5

0.8 0.75 0.7 0.65 0.6

0.8 0.7 0.6 0.5 0.4

0.8 0.7 0.6 0.5 0.4 0.3

Figure 22. Annual median NDVI of HGU A for period 2008 until 2014.

62

0.8 0.7 0.6 0.5 0.4

0.8 0.7 0.6 0.5

0.8 0.7 0.6 0.5

0.8 0.7 0.6 0.5

0.9 0.8 0.7 0.6 0.5 0.4

0.8 0.7 0.6 0.5

0.8 0.7 0.6 0.5

Figure 23. Annual median NDVI of HGU B for period 2008 until 2014.

63

0.8 0.7 0.6 0.5 0.4 0.3

0.8 0.7 0.6 0.5 0.4 0.3

0.8 0.7 0.6 0.5 0.4

0.9 0.8 0.7 0.6 0.5 0.4 0.3

0.8 0.7 0.6 0.5 0.4 0.3

0.8 0.7 0.6 0.5 0.4 0.3

0.8 0.7 0.6 0.5 0.4 0.3

Figure 24. Annual median NDVI of HGU C for period 2008 until 2014.

64

0.8 0.75 0.7 0.65 0.6

0.75 0.7 0.65 0.6 0.55

0.8 0.75 0.7 0.65 0.6 0.55 0.5

0.8 0.7 0.6 0.5 0.4 0.3

0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55

0.8 0.7 0.6 0.5

0.7 0.6 0.5 0.4

Figure 25. Annual median NDVI of HGU E for period 2008 until 2014.

65

Appendix II. The forest mask of HGU A, B, C, D, and E.

Non-forest Forest

(a) (b) (c) Figure 26. The determination of HGU A’s forest mask: (a) SPOT 6 image acquired on 7 January 2013, (b) non-forest and forest class, (c) final forest mask.

Non forest

Figure 27. Forest mask of HGU B, C, D, and E resulted from classification of SPOT 6 images acquired on 7 January 2013.

66

Appendix III. The BFAST monitor output of HGU A, B, C, D, and E.

Figure 28. The output of BFAST monitor for HGU A.

Figure 29. The output of BFAST monitor for HGU B.

Figure 30. The output of BFAST monitor for HGU C.

67

Figure 31. The output of BFAST monitor for HGU D.

Figure 32. The output of BFAST monitor for HGU E.

68

Appendix IV. The accuracy assessment of BFAST monitor in detecting land clearance. 2013

MT = -0.050

MT = -0.025

MT = 0

MT = -0.075

MT = -0.1

2014 2014

Figure 33. Land clearance detected by BFAST monitor for HGU A. The upper five maps display the land clearance class determined by combining the breakpoints and five different magnitude thresholds (MT). The SPOT 6 image for the five upper maps was acquired on 7 January 2013. The lower map shows the SPOT 6 image acquired on 17 June 2014. This image shows the land cover of HGU A in 2014 and was used to validate the detected land clearance class of the upper five maps. 69

Table 15. The BFAST monitor accuracy assessment in detecting land clearance for HGU A. Land clearance = breakpoint and magnitude < 0 Land clearance No land clearance Total 26 17 43 Land clearance 4 15 19 No land clearance 30 32 62 Total Producer’s (%) 86.67 46.87 Overall (%) Land clearance = breakpoint and magnitude < -0.025 Land clearance No land clearance Total 25 11 36 Land clearance 4 15 19 No land clearance 29 26 55 Total Producer’s (%) 86.21 57.69 Overall (%) Land clearance = breakpoint and magnitude < -0.05 Land clearance No land clearance Total 24 7 31 Land clearance 4 15 19 No land clearance 28 22 50 Total Producer’s (%) 85.71 68.18 Overall (%) Land clearance = breakpoint and magnitude < -0.075 Land clearance No land clearance Total 23 4 27 Land clearance 4 15 19 No land clearance 27 19 46 Total Producer’s (%) 85.19 78.94 Overall (%) Land clearance = breakpoint and magnitude < -0.1 Land clearance No land clearance Total 19 3 22 Land clearance 4 15 19 No land clearance 23 18 41 Total Producer’s (%) 82.61 83.33 Overall (%)

User’s (%) 60.47 78.94 66.13 User’s (%) 69.44 78.94 72.73 User’s (%) 77.42 78.94 78.00 User’s (%) 85.19 78.94 82.61 User’s (%) 86.36 78.94 82.93

70

2013

MT = -0.025

MT = 0

MT = -0.050

MT = -0.075

MT = -0.1

2014 2014

Figure 34. Land clearance detected by BFAST monitor for HGU C. The upper five maps display the land clearance class determined by combining the breakpoints and five different magnitude thresholds (MT). The SPOT 6 image for the five upper maps was acquired on 7 January 2013. The lower map shows the SPOT 6 image acquired on 17 June 2014. This image shows the land cover of HGU C in 2014 and was used to validate the detected land clearance class of the upper five maps. 71

Table 16. The BFAST monitor accuracy assessment in detecting land clearance for HGU C. Land clearance = breakpoint and magnitude < 0 Land clearance No land clearance Total 20 84 104 Land clearance 12 88 100 No land clearance 32 172 204 Total Producer’s (%) 62.50 51.16 Overall (%) Land clearance = breakpoint and magnitude < -0.025 Land clearance No land clearance Total 20 73 93 Land clearance 12 88 100 No land clearance 32 161 193 Total Producer’s (%) 62.50 54.66 Overall (%) Land clearance = breakpoint and magnitude < -0.05 Land clearance No land clearance Total 20 63 83 Land clearance 12 88 100 No land clearance 32 151 183 Total Producer’s (%) 62.50 58.28 Overall (%) Land clearance = breakpoint and magnitude < -0.075 Land clearance No land clearance Total 15 16 31 Land clearance 12 88 100 No land clearance 27 104 131 Total Producer’s (%) 55.56 84.62 Overall (%) Land clearance = breakpoint and magnitude < -0.1 Land clearance No land clearance Total 10 5 15 Land clearance 12 88 100 No land clearance 22 93 115 Total Producer’s (%) 45.45 94.62 Overall (%)

User’s (%) 19.23 88.00 52.94 User’s (%) 21.51 88.00 55.96 User’s (%) 24.09 88.00 59.02 User’s (%) 48.39 88.00 78.63 User’s (%) 66.67 88.00 85.21

72

MT = -0.025

MT = 0

MT = -0.075

MT = -0.050

MT = -0.1

2014

Figure 35. Land clearance detected by BFAST monitor for HGU D. The upper five maps display the land clearance class determined by combining the breakpoints and five different magnitude thresholds (MT). The SPOT 6 image for the five upper maps was acquired on 7 January 2013. The lower map shows the SPOT 6 image acquired on 17 June 2014. This image shows the land cover of HGU D in 2014 and was used to validate the detected land clearance class of the upper five maps. 73

Table 17. The BFAST monitor accuracy assessment in detecting land clearance for HGU D. Land clearance = breakpoint and magnitude < 0 Land clearance No land clearance Total 22 16 38 Land clearance 10 50 60 No land clearance 32 66 98 Total Producer’s (%) 68.75 75.75 Overall (%) Land clearance = breakpoint and magnitude < -0.025 Land clearance No land clearance Total 21 14 35 Land clearance 10 50 60 No land clearance 31 64 95 Total Producer’s (%) 67.74 78.12 Overall (%) Land clearance = breakpoint and magnitude < -0.05 Land clearance No land clearance Total 20 12 32 Land clearance 10 50 60 No land clearance 30 62 92 Total Producer’s (%) 66.67 80.64 Overall (%) Land clearance = breakpoint and magnitude < -0.075 Land clearance No land clearance Total 20 7 27 Land clearance 10 50 60 No land clearance 30 57 87 Total Producer’s (%) 66.67 87.71 Overall (%) Land clearance = breakpoint and magnitude < -0.1 Land clearance No land clearance Total 19 3 22 Land clearance 10 50 60 No land clearance 29 53 82 Total Producer’s (%) 65.51 94.33 Overall (%)

User’s (%) 57.89 83.33 73.46 User’s (%) 60.00 83.33 74.73 User’s (%) 62.50 83.33 76.08 User’s (%) 74.07 83.33 80.45 User’s (%) 86.36 83.33 84.14

74

2013

MT = -0.025

MT = 0

MT = -0.075

MT = -0.050

MT = -0.1

2014

2014

Figure 36. Land clearance detected by BFAST monitor for HGU E. The upper five maps display the land clearance class determined by combining the breakpoints and five different magnitude thresholds (MT). The SPOT 6 image for the five upper maps was acquired on 7 January 2013. The lower map shows the SPOT 6 image acquired on 17 June 2014. This image shows the land cover of HGU E in 2014 and was used to validate the detected land clearance class of the upper five maps. 75

Table 18. The BFAST monitor accuracy assessment in detecting land clearance for HGU E. Land clearance = breakpoint and magnitude < 0 Land clearance No land clearance Total 13 12 25 Land clearance 4 52 56 No land clearance 17 64 81 Total Producer’s (%) 76.47 81.25 Overall (%) Land clearance = breakpoint and magnitude < -0.025 Land clearance No land clearance Total 11 10 21 Land clearance 4 52 56 No land clearance 15 62 77 Total Producer’s (%) 73.33 83.87 Overall (%) Land clearance = breakpoint and magnitude < -0.05 Land clearance No land clearance Total 11 9 20 Land clearance 4 52 56 No land clearance 15 61 76 Total Producer’s (%) 73.33 85.25 Overall (%) Land clearance = breakpoint and magnitude < -0.075 Land clearance No land clearance Total 8 4 12 Land clearance 4 52 56 No land clearance 12 56 68 Total Producer’s (%) 66.67 92.86 Overall (%) Land clearance = breakpoint and magnitude < -0.1 Land clearance No land clearance Total 7 2 9 Land clearance 4 52 56 No land clearance 11 54 65 Total Producer’s (%) 63.64 96.29 Overall (%)

User’s (%) 52.00 92.86 80.25 User’s (%) 52.38 92.86 81.82 User’s (%) 55.00 92.86 82.89 User’s (%) 66.67 92.86 88.24 User’s (%) 77.78 92.86 90.77

76

Appendix V. The BFAST monitor plots for single pixel time series.

Figure 37. Some examples of BFAST monitor plot for single pixel time series in detecting disturbances/land clearances.

77