Modeling Future Urban Sprawl and Landscape Change in the Laguna de Bay Area, Philippines

: This study uses a spatially-explicit land-use/land-cover (LULC) modeling approach to model and map the future (2016–2030) LULC of the area surrounding the Laguna de Bay of Philippines under three different scenarios: ‘business-as-usual’, ‘compact development’


Introduction
Urban sprawl is occurring at an accelerated pace in many developing countries worldwide due to rapid global economic and population growth coupled with globalization.Currently, 54% of the world's population lives in urban areas, and the United Nations has predicted that by 2050, 66% of the world's population will live in urban areas [1].This rapid increase in urban population has forced nations to meet the changing demands for necessities such as food, energy, land, and water.A major concern related to this urban sprawl is land-use (LU)/land-cover (LC) change, which can dramatically alter the landscape in areas with high rates of urban expansion [2].These LULC changes are often based on the plans of local governments to increase economic development and to support their growing populations.However, such plans may fail to consider other factors, including climate conditions, water resources, and food security [3,4].Since population increases are expected to Land 2017, 6, 26 3 of 21 categorized into six types of approaches [31]; machine learning and statistical, cellular automata, sector-based economic, spatially disaggregated economic, agent-based, and hybrid approaches.In this work, we focused on the machine learning category, namely the Multilayer-Perceptron (MLP) Artificial Neural Network (ANN) approach, which is gaining attention for modeling LULC changes.ANN are powerful tools for modeling complex behaviors [32] (e.g, relationships between land transitions and their driving forces), and the usage of ANN for LULC change modeling has increased in recent years.As one example, Grekousis et al. [33] demonstrated the use of ANN to model future urban growth based on demographic time-series data.Triantakonstantis and Stathakis [32] used MLP for modeling future LULC transition probabilities based on information on past LULC changes and geomorphic drivers such as elevation, slope, and distance variables from specific land features, etc.Similar methods can be seen in a number of studies, although the number and types of driver variables utilized vary [27,[34][35][36][37][38].
Past studies have been conducted at various geographical locations and spatial extents, e.g., the city level or provincial level.However, not many works have used the MLP method for modeling complex larger regions.In our case, we have focused on a larger ecological scale that encompasses multiple municipalities and provinces, each of which have different development policies, infrastructure, topography, climatic conditions, etc.Therefore, the MLP method for modeling the LULC transitions throughout the region becomes a strong decision tool, even when prior knowledge is lacking [31,37], and this is one of its advantages compared to other methods like SLEUTH (Slope, Land use map, Excluded area, Urban area, Transportation map, Hillside area), which require coefficient values to be set [39], and other cellular automata methods that require a suitability map [40] based on prior knowledge of change behavior.MLP can also handle a large number of data sets, so it could be a good predictor for recognizing the patterns of the changes in the area.
Our main objective in this study is to model the future LULC changes in the river basins that drain into the Laguna de Bay up to year 2030.This work implements the MLP with the Markov Chain method embedded in the Land Change Modeler (LCM) of the TerrSet software package [41].This model is based on transition probabilities calculated using historical LULC change data and other freely/openly available geospatial datasets.

Study Area
The Philippines is one of the most rapidly developing countries in Asia.In particular, the Metropolitan Manila area has experienced large and rapid LULC changes owing to urban area expansion [2,6].The Laguna de Bay, located just southeast of Manila, is the largest and most important and dynamic lake in the Philippines owing to its vital economic, political, and socio-cultural significance.With a surface area of 900 km 2 , this lake is also one of the largest in Southeast Asia.From the 21 major river systems, more than 100 rivers that traverse the 292,000 ha watershed flow into the lake.The study area includes the municipalities surrounding the Laguna de Bay, located southeast of Metro Manila.Figure 1 shows that Metro Manila, situated to the northwest, has a denser urban concentration; while smaller, less densely populated cities and municipalities can be seen in the west and southwest areas surrounding the Laguna de Bay.The population of this area was estimated to be about 15 million in 2010 [42].Agricultural lands are distributed mainly along the southwestern, southeastern, and northeastern shores of the lake.Large forest areas are located on the east side of the lake with mountain ranges (300-600 m), the highest peak of which, Mt.Banahao (2170 m), is located southeast of the lake.Just south of Laguna de Bay, Mt.Makilling (1090 m) is seen.The area is broad and the climate varies along different provinces of the area.Average annual temperatures are cooler in the mountainous areas (23 • C) than in the lower altitude plains and cities (25 • C and 27 • C respectively).Annual precipitation ranges from over 3000 mm in the east mountainous areas to 1900 mm in the western area around Manila bay.The area in focus is approximately 60 km × 80 km in the north-south and east-west directions, respectively, and it was chosen to visualize how urban sprawl will affect the LULC in the areas surrounding the Laguna de Bay.Therefore, we have ignored the center and northern areas of Metro Manila, even though it also experienced urban sprawl.

Overview
As shown in Figure 2, the overall flowchart of the study details three main steps in generating maps of future LULC.The first step is to gather evidence of past LULC transitions by identifying the historical LULC changes in the area.For this, a map of recent LULC change from 2007 to 2015 was generated by using optical and synthetic aperture radar (SAR) satellite images from 2007 and 2015 and automated image classification techniques.Four LULC classes, including Built-up, Forests, Crop-Grass, and Water Bodies, were mapped for each year, and the LULC changes between 2007 and 2015 were identified by overlaying the two maps.The full details of the LULC change mapping methodology employed in this study are given in Johnson et al [44], although brief information will be stated about the methods and the result of the developed LULC change map.In the second step, drivers of these historical LULC changes were identified by using various ancillary spatial datasets containing demographic, topographic, and climate information.In the third and final step, the future LULC of the area (2030) was modeled and mapped by using Markov Chain analysis.
In addition to the main procedure, the validity of the LULC change model was examined by comparing a 'simulated 2015 LULC map' with a 'reference 2015 LULC map'.The 'simulated 2015 LULC map' was generated using similar methods to those mentioned above.First, LULC maps from the years 2007 and 2010 were utilized along with the ancillary spatial datasets to model the LULC conditions in the year 2015 LULC (i.e., simulated 2015 LULC map).This map was then compared with the actual (i.e.reference) 2015 LULC map.

Overview
As shown in Figure 2, the overall flowchart of the study details three main steps in generating maps of future LULC.The first step is to gather evidence of past LULC transitions by identifying the historical LULC changes in the area.For this, a map of recent LULC change from 2007 to 2015 was generated by using optical and synthetic aperture radar (SAR) satellite images from 2007 and 2015 and automated image classification techniques.Four LULC classes, including Built-up, Forests, Crop-Grass, and Water Bodies, were mapped for each year, and the LULC changes between 2007 and 2015 were identified by overlaying the two maps.The full details of the LULC change mapping methodology employed in this study are given in Johnson et al. [44], although brief information will be stated about the methods and the result of the developed LULC change map.In the second step, drivers of these historical LULC changes were identified by using various ancillary spatial datasets containing demographic, topographic, and climate information.In the third and final step, the future LULC of the area (2030) was modeled and mapped by using Markov Chain analysis.
In addition to the main procedure, the validity of the LULC change model was examined by comparing a 'simulated 2015 LULC map' with a 'reference 2015 LULC map'.The 'simulated 2015 LULC map' was generated using similar methods to those mentioned above.First, LULC maps from the years 2007 and 2010 were utilized along with the ancillary spatial datasets to model the LULC conditions in the year 2015 LULC (i.e., simulated 2015 LULC map).This map was then compared with the actual (i.e., reference) 2015 LULC map.

LULC Maps of 2007, 2010 and 2015
To develop a categorical LULC change map of the study area, we utilized optical (Landsat 5 and Landsat 8) and synthetic aperture radar (ALOS PALSAR-1 and PALSAR-2) satellite images from the years 2007, 2010, and 2015 and classified the pixels in the images from each year into one of four LULC classes (Built-up, Crop-Grass, Trees, and Water) using a semi-unsupervised classification approach [44].Crop-Grass includes cropland, paddy fields, grassland, and pasture (however the majority lies within paddy and other agriculture).Trees includes forest and agroforestry plantations (e.g.coconut, banana, etc.).This four class LULC classification system, although relatively simple, was chosen because it was representative of the area and allowed us to maintain a relatively high LULC change mapping accuracy (adding more specific LULC classes usually decreases LULC mapping accuracy, and LULC change mapping accuracy even more so due to error propagation).The overall accuracy of the 2007-2015 LULC change map was estimated as 90.2%.

Evidence of LULC Change Transition
By using the two different periods of the LULC change maps, the net change in area of each LULC class was calculated, and the spatial distributions of all of the LULC changes were analyzed.The areas of transition and persistence of each LULC type within the 2007-2015 analysis was for both training and validation data with all of the driver variables when performing transition sub-modeling.

Collection and Processing of Data on Potential Driver Variables
Spatial data related to various potential drivers of LULC change were collected via the Internet.Only datasets which were openly available online were used, so our modeling approach can easily be replicated by other researchers.The drivers influencing LULC change processes are extremely diverse as well as highly variable from one location to another [45,46].What is known is that the changes are typically the results of the local population's responses to economic opportunities [45], which gives relevance to various contextual features such as the distance from a location to nearby infrastructural features like major roads, town centers, and so forth, and many works show these kind of factors to be used as potential drivers for calibrating the change probability.These context features are also considered in this work for the model calibration.However, as mentioned before, we are not so aware of all the LULC change drivers in the area, so we want to shed some light on the question of 'What factors influence the LULC transitions?' Thus, data for a large number of possible

LULC Maps of 2007, 2010 and 2015
To develop a categorical LULC change map of the study area, we utilized optical (Landsat 5 and Landsat 8) and synthetic aperture radar (ALOS PALSAR-1 and PALSAR-2) satellite images from the years 2007, 2010, and 2015 and classified the pixels in the images from each year into one of four LULC classes (Built-up, Crop-Grass, Trees, and Water) using a semi-unsupervised classification approach [44].Crop-Grass includes cropland, paddy fields, grassland, and pasture (however the majority lies within paddy and other agriculture).Trees includes forest and agroforestry plantations (e.g., coconut, banana, etc.).This four class LULC classification system, although relatively simple, was chosen because it was representative of the area and allowed us to maintain a relatively high LULC change mapping accuracy (adding more specific LULC classes usually decreases LULC mapping accuracy, and LULC change mapping accuracy even more so due to error propagation).The overall accuracy of the 2007-2015 LULC change map was estimated as 90.2%.

Evidence of LULC Change Transition
By using the two different periods of the LULC change maps, the net change in area of each LULC class was calculated, and the spatial distributions of all of the LULC changes were analyzed.The areas of transition and persistence of each LULC type within the 2007-2015 analysis was for both training and validation data with all of the driver variables when performing transition sub-modeling.

Collection and Processing of Data on Potential Driver Variables
Spatial data related to various potential drivers of LULC change were collected via the Internet.Only datasets which were openly available online were used, so our modeling approach can easily be replicated by other researchers.The drivers influencing LULC change processes are extremely diverse as well as highly variable from one location to another [45,46].What is known is that the changes are typically the results of the local population's responses to economic opportunities [45], which gives relevance to various contextual features such as the distance from a location to nearby infrastructural features like major roads, town centers, and so forth, and many works show these kind of factors to be used as potential drivers for calibrating the change probability.These context features are also considered in this work for the model calibration.However, as mentioned before, we are not so aware of all the LULC change drivers in the area, so we want to shed some light on the question of 'What factors influence the LULC transitions?' Thus, data for a large number of possible drivers was collected to test which variables had the greatest levels of influence on LULC changes (a brief explanation is made later for why each variable was considered as potentially relevant).Table 1 shows the complete list of collected data related to these drivers; not all of these datasets were selected for the final LULC modeling.Gridded climate data with a 1 km resolution was obtained from the Philippine GIS Data Clearinghouse (PhilGIS) [47] (the original global dataset is distributed at WorldClim [52]).We included all of the climate data that was processed in the form of bioclimatic variables [53,54].These climate variables can be considered drivers of LULC change because the agro-climatic zones with different climatic conditions can affect the suitability of agricultural lands for its productivity [55]; thus, no patterns of change in the preferred area might be identified relative to this factor.This is one challenge in our work since not many related studies implement climatic information in their calibration.The reason can be considered that, depending on the scale of the study area, the spatial variation of climatic factors might be too small to have any influence on LULC changes.However, our study area encompasses areas with different climatic conditions, which may allow us to identify if climate is affecting LULC change.Topographic data including elevation, slope, and aspect were obtained from the Shuttle Radar Topography Mission (SRTM) 1 s (30 m) digital elevation model (DEM) [48].Topography is often a significant driver of LULC change [56] because areas with steep slopes are typically more difficult and thus less likely than flatter areas to be converted to built-up land or cropland.The SRTM DEM contains only elevation information; therefore, gridded slope and aspect data were generated by using the TerrSet software package [41].Road and waterway data through March 03, 2016, were collected from OpenStreetMap (OSM) [49], and a 25 m grid map containing the Euclidean distance of each pixel to the nearest road/waterway was calculated by using TerrSet software [41].The distance from roads in addition to other various LULC context features are often drivers of LULC change because more developed road networks are found with a greater rate of conversion [57].This study focuses not only on the roads in general but also on different types/functions of roads by separating roads into detailed classes.It is expected that the importance of different road types could further distinguish the patterns of transitions throughout the study region.The case is similar for the waterways.Nightlight intensity information in the form of monthly average radiance composite images was obtained from the Visible Infrared Imaging Radiometer Suite (VIIRS) and Global Defense Meteorological Satellite Program-Operational Linescan System (DMSP-OLS) nighttime lights time series dataset [50].The nightlight changes from 2007 to 2015 were computed by performing radiometric normalization of the images to ensure that the 2007 data had a radiance range similar to that of the 2015 data.The differences in radiance at each pixel location were then calculated.Nightlight intensity was considered a driver variable because it is strongly correlated with economic activity and the gross domestic product (GDP) [58,59].Changes in nightlight intensity over time can also be considered an indicator of economic growth.Polygon data on the locations of protected areas were collected from the PhilGIS [47].By using this dataset, the distance from each pixel to the nearest protected area was calculated in TerrSet [41], and all pixels located within a protected area were assumed to experience no LULC conversions in the LULC change modeling process.Unique LULC data such as golf course information were also collected and used to generate the distance information from those features.Compared to other unique LULC features such as markets and town centers, the location of a golf course would mostly not change over time, and there would be a lesser chance for new land areas to be developed as golf courses, keeping the consistency of the patterns; therefore we have used those data.Gridded population data were obtained from WorldPop [51].This data is based on census population counts at the Barangay level but were downscaled to a 100 m × 100 m grid level utilizing various other spatial datasets, as outlined in Stevens et al. [60] and Linard et al. [61].Population growth was calculated on the basis of population change between 2007 and 2015.All of these gridded datasets were resampled from their original resolutions to a 25 m resolution by using a cubic convolution resampling approach to match the resolution of the LULC change map (25 m) of the area.For the visual interpretation of the variables used in the model, summarized images are given in Figure S1 of the supplementary materials.
The explanatory power of all of these variables in relation to different LULC transitions was computed and examined by using Cramer's V [38].Also known as Cramer's Coefficient (V), this method is used for quantifying the explanatory power of each variable, which is an optional quick test used to determine whether the variables are worthy of consideration in the model [29].The final variables used for the modeling will be considered on certain criteria of this value.

Processing Transition Sub-Models (MLP)
MLP, a type of ANN method widely used for modeling complex behaviors and patterns, uses the back propagation algorithm to learn the characteristics of all the factors influencing the LULC transitions.Several studies show the advantages of MLP compared to logistic regression and other empirical models [62,63].Further details of the MLP algorithm can be found in Riccioli et al. [37].In our study, MLP's ability to handle a large number of input variables (some of which may be irrelevant and/or highly correlated with one another) in the model calibration process was very useful, as it allowed us to investigate over 20 explanatory variables.We focused on modeling the changes of three LULC classes; Built-up, Crop-Grass, and Trees.For all of the variables measuring 'distance from' a pixel to some geographic feature (e.g., road, built-up area, etc.), the distances were recalculated at a one year interval using that year's modeled LULC map.A random sample of 10,000 pixels from the 2007-2015 LULC change map was used for the building model.Of these, 50% were used for training and 50% were used for testing through a cross-validation process.

Change Modeling (Three Scenarios)
The probability of changes occurring in different years in the future was calculated using Markov Chain analysis, a technique for predictive change modeling that is able to model future changes based on past changes.On the basis of the observed data between the two periods (2007 and 2015 in our case), the Markov Chain computes the probability that a pixel will change from one LULC type to another within a specified period [64].Table 2 shows the matrix of the probability that each LULC category will change to every other category (base rate), which is known as the transition probability.In this method, the probability is determined by the actual changes shown in the developed LULC map; further details have been reported by Takada et al. [65].The target year of the modeling in the present study was set to 2030; the transition for each year was also produced to review the continuous dynamic changes in the study area.A total of three scenarios were output for the comparison.The base scenario considers the transition rate to be the same as the 2007 to 2015 change rate (i.e., a 'business-as-usual scenario').The second scenario is if the development policy changes and the rate of LULC change reduces to half of the 2007-2015 transition rate (i.e., a 'compact development scenario').The third scenario is if the rate of LULC change further accelerates to twice the 2007-2015 rate (i.e., a 'high sprawl scenario').The deceleration and acceleration of the transition rates are controlled by simply half and double the values for each changing class in the transition probability matrix, respectively.
During the process of the simulation, the data of protected areas are used as constrain maps to control the process of transitions.The transition potentials associated with each transition are multiplied by the constraints map [64], so a value of 1 means unconstrained, while a values near to 0 acts as a disincentive and above 1 acts as an incentive.The protected area is given the value of very near 0 (i.e., 0.01).

Validation of Modeled Map
We assessed the validity of the model by comparing the simulated 2015 LULC (i.e., derived from LULC modeling using the 2007-2010 LULC change data) with the reference (i.e., mapped) 2015 LULC.Utilizing the 2007 and 2010 LULC maps, a LULC change model was calibrated using the same driver variables as used in the 2007-2015 calibration.Using the 2007-2010 model, the year 2015 LULC was simulated, and the projected 2015 map was compared with the reference 2015 map.Two statistical indexes were calculated for the validation; Figure of Merit (FoM) [66] and Kappa Index [67].FoM determines the accuracy of LULC hits (model predicted change and actually observed change) compared with the sum of hits, misses (model predicted persistence but is observed change), and false alarms (model predicted change but it observed persistence), giving 0% for no match between the modeled and the reference map and 100% for a perfect match.Kappa Index is widely used in the remote sensing society for assessing the reliability of the map.Other than the standard Kappa (K standard ), Kappa for no ability (K no ) and Kappa for location (K location ) are computed.K no considers Land 2017, 6, 26 9 of 21 and fixes the major problems of the standard Kappa, wherein it penalizes for large quantification error and fails to reward the simulation for specifying quantity [67].K location indicates how well the grid cells are located on the landscape [64,67].Kappa values range from −1 (no agreement) to 1 (perfect agreement).The water bodies were masked out for the computation.

2007-2015 LULC Changes
By using the Cross-Tabulation module, the transitions of each class from 2007 to 2015 were computed, as shown in Figure 3a.The gains and the losses for each LULC type, in ha, are also shown in Figure 3b.The transition map confirmed that the majority of changes in the Built-Up class are attributed to decreases in the Crop-Grass area, although a significant number of Trees areas are expected to become Built-Up land outside the protected forest areas.
LULC.Utilizing the 2007 and 2010 LULC maps, a LULC change model was calibrated using the same driver variables as used in the 2007-2015 calibration.Using the 2007-2010 model, the year 2015 LULC was simulated, and the projected 2015 map was compared with the reference 2015 map.Two statistical indexes were calculated for the validation; Figure of Merit (FoM) [66] and Kappa Index [67].FoM determines the accuracy of LULC hits (model predicted change and actually observed change) compared with the sum of hits, misses (model predicted persistence but is observed change), and false alarms (model predicted change but it observed persistence), giving 0% for no match between the modeled and the reference map and 100% for a perfect match.Kappa Index is widely used in the remote sensing society for assessing the reliability of the map.Other than the standard Kappa (Kstandard), Kappa for no ability (Kno) and Kappa for location (Klocation) are computed.Kno considers and fixes the major problems of the standard Kappa, wherein it penalizes for large quantification error and fails to reward the simulation for specifying quantity [67].Klocation indicates how well the grid cells are located on the landscape [64,67].Kappa values range from −1 (no agreement) to 1 (perfect agreement).The water bodies were masked out for the computation.however, smaller but important signs were detected at different parts of the surrounding environment, which possibly has relationships with the smaller cities located nearby.The changes of Crop-Grass to Trees were located mostly at the southern side of the study area, which is near protected areas.Thus, the Grass-Shrub-type LC is slowly changing into denser and taller vegetation due to the absence of effects from human activities.
Crop-Grass to Built-Up, Trees to Built-Up, Trees to Crop-Grass, and Crop-Grass to Trees.This method enables identification and understanding of the spatial trends of the transition, which can provide a better comprehension of the sites of different changes at different spatial locations.The assumption shows the change patterns that occurred in the area from 2007 to 2015.For the Built-Up class, significant changes were detected from north to south along the west side of the Laguna de Bay.High probabilities of transition were detected from the Trees (Crop-Grass) class to Built-Up class north (west) of the center of the study area.The trends for the Trees class changing into Crop-Grass were most dominant north of the center of the study area, followed secondly by the southwest area and thirdly by the south to southeast area.The largest changes were recognized at the southwestern side; however, smaller but important signs were detected at different parts of the surrounding environment, which possibly has relationships with the smaller cities located nearby.The changes of Crop-Grass to Trees were located mostly at the southern side of the study area, which is near protected areas.Thus, the Grass-Shrub-type LC is slowly changing into denser and taller vegetation due to the absence of effects from human activities.

Potential Explanatory Power of Driver Variables
Table 3 shows the driver variables selected for modeling of the transition potential.Each variable has a potential explanatory power that describes the strength of its relationship to the actual transition of the classes and is computed by contingency table analysis.V takes values from 0 to 1.

Potential Explanatory Power of Driver Variables
Table 3 shows the driver variables selected for modeling of the transition potential.Each variable has a potential explanatory power that describes the strength of its relationship to the actual transition of the classes and is computed by contingency table analysis.V takes values from 0 to 1. Values near 0 show little association between variables, whereas those near 1 indicate strong association.A value of about 0.15 contains little information for explanation; more than 0.4 is considered to be a good variable [64].Only the final variables selected for the model are listed in the table.The selection criteria was that a variable had either (i) an overall V value greater than 0.3 or (ii) V values greater than 0.15 for all of the individual LULC classes.For example, although for Dist_Protect the overall V value is below 0.3, each class shows values above 0.15.
Table 4 shows the priority of the drivers compared with other variables, showing which are the most and least influential.For model (b), the most (least) influential was DEM (NL_2015).For model (c), the most (least) influential was Pop_2007 (Road_Dist2).For model (d), the most (least) influential was Slope (NL_ch).The accuracy of each transition model in explaining its overall power for detecting the correct changes are Crop-Grass to Built-Up, 74.21%; Crop-Grass to Trees, 70.37%; Trees to Built-Up, 90.91%; and Trees to Crop-Grass, 74.31%.

LULC Change Modeling and Its Landscape
For visual interpretation of the dynamic changes in LULC classes for each stage of the modeling (each year of 2016-2030), Video S1 is provided in the supplementary material (for business-as-usual scenario).Here, we discuss only the beginning (2015), middle (2023), and end (2030) years of the LULC map. Figure 5 shows that changes occurred from the modeling at the Laguna de Bay region.A few significant characteristic trends of change depending on each LULC class can be identified.For the Built-Up classes, the first large change is the expansion of urban areas spreading more southward and expansion at the west side of the Laguna de Bay.This type of trend has occurred in the past.Google Earth images from the 1980s in those regions [68] show strong evidence of rapid LULC change in the southern part of Metro Manila and at the west side of the Laguna de Bay.Built-Up areas have expanded in the southwestern part of the study area near a smaller lake and show development along road infrastructures at the east side of this lake.At the east side of the Laguna de Bay, significantly fewer LULC changes have occurred compared with the west side because the east side of the lake consists mainly of rural areas with low population density and less infrastructure that are thus not affected by rapid development.

Road_Dist2
Var. 14    For the Crop-Grass class, the main changes were detected in three areas.The first includes changes along with the development of Built-Up areas at the west side of Laguna de Bay.In addition, the Crop-Grass class area has expanded and has ensured a considerable amount of area by invading the Trees class.The second site is at the east coast of the Laguna de Bay, where Built-Up areas have not changed significantly, although the Crop-Grass class has again ensured its area by invading large amounts of Trees class areas.The third location is at the center north of the Laguna de Bay, where a similar trend of Crop-Grass invading the Trees class was noted.Future scenarios can change depending on whether countermeasures for protecting the forests were planned and implemented.Similar to using the protected areas as constraining images for unchanges, probability maps for specific regions can be applied to reduce the transition probabilities of those areas.These types of planning scenarios for the forested areas can be considered a strong decision tool for planning REDD+ (Reducing Emissions from Deforestation and forest Degradation in developing countries) actions [69] because total changes in forests can be easily calculated quantitatively, enabling sound estimation of the CO2 uptakes in those regions.
For the Trees class, the mountainous area at the east side of the Laguna de Bay and the protected areas show increments of Trees Class transition from the Crop-Grass.This occurrence depends on the ages of the people living in those rural areas, which have higher elevations and rugged terrain.When the residents become older, the land becomes more difficult to access.Agricultural lands would then begin to change into abandoned areas; thus, the trend in Crop-Grass LULC type would decrease.Other possible factors include the impact of the National Greening Program of the government [70].In other regions, the majority of the Trees class trend shows a decrease in total area owing to its transition to the Built-Up and Crop-Grass classes.

LULC Change Statistics
Figure 6 shows the accumulated increase in land area of each transition class from the starting year of 2015 to 2030 in yearly increments for the business-as-usual scenario, compact development scenario, and high sprawl scenario.In the business-as-usual scenario, most of the LULC change For the Crop-Grass class, the main changes were detected in three areas.The first includes changes along with the development of Built-Up areas at the west side of Laguna de Bay.In addition, the Crop-Grass class area has expanded and has ensured a considerable amount of area by invading the Trees class.The second site is at the east coast of the Laguna de Bay, where Built-Up areas have not changed significantly, although the Crop-Grass class has again ensured its area by invading large amounts of Trees class areas.The third location is at the center north of the Laguna de Bay, where a similar trend of Crop-Grass invading the Trees class was noted.Future scenarios can change depending on whether countermeasures for protecting the forests were planned and implemented.Similar to using the protected areas as constraining images for unchanges, probability maps for specific regions can be applied to reduce the transition probabilities of those areas.These types of planning scenarios for the forested areas can be considered a strong decision tool for planning REDD+ (Reducing Emissions from Deforestation and forest Degradation in developing countries) actions [69] because total changes in forests can be easily calculated quantitatively, enabling sound estimation of the CO 2 uptakes in those regions.
For the Trees class, the mountainous area at the east side of the Laguna de Bay and the protected areas show increments of Trees Class transition from the Crop-Grass.This occurrence depends on th ages of the people living in those rural areas, which have higher elevations and rugged terrain.When the residents become older, the land becomes more difficult to access.Agricultural lands would then begin to change into abandoned areas; thus, the trend in Crop-Grass LULC type would decrease.Other possible factors include the impact of the National Greening Program of the government [70].In other regions, the majority of the Trees class trend shows a decrease in total area owing to its transition to the Built-Up and Crop-Grass classes.

LULC Change Statistics
Figure 6 shows the accumulated increase in land area of each transition class from the starting year of 2015 to 2030 in yearly increments for the business-as-usual scenario, compact development scenario, and high sprawl scenario.In the business-as-usual scenario, most of the LULC change classes showed increases in total area and different growing rates; however, the Crop-Grass to Trees class showed a decrease in area beginning in 2026.The Trees to Crop-Grass and Crop-Grass to Trees classes showed a polynomial-like increase, resulting in increases of 26,409 ha and 21,166 ha, respectively, in 2030.The Crop-Grass to Built-Up and Trees to Built-Up classes showed a linear increase of land area, resulting in increases of 14,137 ha and 3946 ha, respectively, in 2030.These changes are limited to the study area.More meaningful values might be extracted when areas are divided according to administrative boundaries.In-depth information on the LULC changes for each municipal boundary is given in Spreadsheet S1 in the supplementary material.For the compact development scenario, all of the classes show a linear-like increase compared to the business-as-usual scenario.The Trees to Crop-Grass and Crop-Grass to Trees classes shows increases of 18,274 and 16,036 ha respectively, which is about 70% of the business-as-usual scenario changes.The Crop-Grass to Built-up and Trees to Built-up showed 5783 and 2000 ha increases, indicating slightly less than 50% of their areas modeled in the business-as-usual scenario.These values are similar to the modeled 2021 LULC in the business-as-usual scenario.The high sprawl scenario shows a similar pattern to the business-as-usual scenario but with a steeper increase and faster point of decrease in increments for the Crop-Grass to Trees class.The Trees to Crop-Grass and Crop-Grass to Trees classes show increases of 30,021 and 20,182 ha, respectively.The former class shows a 113% increase, but the latter remains at 95% of that of the business-as-usual scenario, showing that more Crop-Grass conversions are occurring.Crop-Grass to Built-up and Trees to Built-up classes also show an increase compared to business-as-usual scenario (26,446 and 7836 ha respectively, approximate double the area of business-as-usual scenario).
Land 2017, 6, 26 14 of 21 maps for specific regions can be applied to reduce the transition probabilities of those areas.These types of planning scenarios for the forested areas can be considered a strong decision tool for planning REDD+ (Reducing Emissions from Deforestation and forest Degradation in developing countries) actions [69] because total changes in forests can be easily calculated quantitatively, enabling sound estimation of the CO2 uptakes in those regions.
For the Trees class, the mountainous area at the east side of the Laguna de Bay and the protected areas show increments of Trees Class transition from the Crop-Grass.This occurrence depends on the ages of the people living in those rural areas, which have higher elevations and rugged terrain.When the residents become older, the land becomes more difficult to access.Agricultural lands would then begin to change into abandoned areas; thus, the trend in Crop-Grass LULC type would decrease.Other possible factors include the impact of the National Greening Program of the government [70].In other regions, the majority of the Trees class trend shows a decrease in total area owing to its transition to the Built-Up and Crop-Grass classes.

2007-2010 Model Validation
Firstly, FoM was carried out by computing hits, misses, and false alarms between the modeled and reference LULC maps of the year 2015.The total number of pixels with hits, misses, and false alarms was 517,671, 642,064, and 964,127 respectively, resulting in a FoM of 24.37%.Secondly, Kappa statistics were computed using the modeled and reference LULC maps.The calculated Kappa statistics were: Kstandard = 0.5825; Kno = 0.6620; and Klocation = 0.6217.A Kappa value of 1 indicates total agreement and 0 indicates totally by chance.This can be interpreted as for example using the Kstandard value that the modeled 2015 map is 58% better than a chance agreement.

Influences of Driver Variables Overview
Generally, the topographical drivers such as DEM and SLOPE influenced the LULC transitions for all classes.As expected, areas of higher elevations and areas with steep slopes tended to experience higher rates of transition to trees, while lower and flatter areas were more likely to convert from tree-covered areas to built-up or agricultural lands.Context drivers showed a strong influence for all LULC change classes other than the Crop-Grass to Built-up class, with the Dist_Canal variable being the top candidate for the other three transition classes.These transitions were thus affected by water availability.Road infrastructure has shown also importance (and also road type), which is logical because they are highly related to people's mobility.The nightlight and population data showed a straightforward result of influencing changes to Built-up land.For example, Crop-Grass to Built-up transitions occurred more frequently in areas with higher populations, while Trees to Built-up transitions occurred more frequently in lower populated areas.Biophysical drivers (climatic variables) did not show as much influence as expected, although some relations were shown for variables such as BIO15 (precipitation seasonality).Our primary expectation was that LULC changes to Crop-Grass might have some relationship with climate, but this was not the case in our study (the slope and population of the nearby area were much more significant drivers).

Scenarios of the Future
Comparing the three different scenarios of future changes that we modeled in this study, although the locations of the transitions did not change significantly (because all three scenarios were based on the same transition model), the quantity and rate of LULC changes differed.If we look at the compact development scenario, where the transition rate is half of business-as-usual

2007-2010 Model Validation
Firstly, FoM was carried out by computing hits, misses, and false alarms between the modeled and reference LULC maps of the year 2015.The total number of pixels with hits, misses, and false alarms was 517,671, 642,064, and 964,127 respectively, resulting in a FoM of 24.37%.Secondly, Kappa statistics were computed using the modeled and reference LULC maps.The calculated Kappa statistics were: K standard = 0.5825; K no = 0.6620; and K location = 0.6217.A Kappa value of 1 indicates total agreement and 0 indicates totally by chance.This can be interpreted as for example using the K standard value that the modeled 2015 map is 58% better than a chance agreement.

Influences of Driver Variables Overview
Generally, the topographical drivers such as DEM and SLOPE influenced the LULC transitions for all classes.As expected, areas of higher elevations and areas with steep slopes tended to experience higher rates of transition to trees, while lower and flatter areas were more likely to convert from tree-covered areas to built-up or agricultural lands.Context drivers showed a strong influence for all LULC change classes other than the Crop-Grass to Built-up class, with the Dist_Canal variable being the top candidate for the other three transition classes.These transitions were thus affected by water availability.Road infrastructure has shown also importance (and also road type), which is logical because they are highly related to people's mobility.The nightlight and population data showed a straightforward result of influencing changes to Built-up land.For example, Crop-Grass to Built-up transitions occurred more frequently in areas with higher populations, while Trees to Built-up transitions occurred more frequently in lower populated areas.Biophysical drivers (climatic variables) did not show as much influence as expected, although some relations were shown for variables such as BIO15 (precipitation seasonality).Our primary expectation was that LULC changes to Crop-Grass might have some relationship with climate, but this was not the case in our study (the slope and population of the nearby area were much more significant drivers).

Scenarios of the Future
Comparing the three different scenarios of future changes that we modeled in this study, although the locations of the transitions did not change significantly (because all three scenarios were based on the same transition model), the quantity and rate of LULC changes differed.If we look at the compact development scenario, where the transition rate is half of business-as-usual scenario, the conversion to Built-up in 2030 does not catch up even to that of 2022 in the business-as-usual scenario.On the other hand, if we double the rate, the development of built-up area would be mostly completed at 2023 compared to 2030 for the business-as-usual scenario.Studies on the potential impacts of these development scenarios on the local environment are still under progress, although we have concerns related to crop production, biodiversity loss (including losses of patches and corridors for wildlife), changes in local climate due to increasing heat fluxes, flood vulnerability, and river and lake water quality.A work by Wijesekara et al. [71] shows good practice of how modeling results can be used as decision/planning tools.
To achieve higher accuracy for the modeling results, two possible factors can be investigated; improving the accuracy of the LULC change map of the region and/or incorporating the zoning policies and development plans of the local governments into the modeling process.The development plans can help determine better scenarios and model areas with higher chances for transitions and the manner of development.A limitation of future modeling, especially using Markov Chains that possess stationary distribution, is that the development rate is considered from past evidence of the changes.If this rate remains constant, the trend is closer to actual trend of future change.However, past information does not always explain future modeling, meaning this can vary depending on decisions by the government or local authorities.

Other Relating Works
In this section, we compare the results of our study with those of other similar studies.First of all, how does our LULC change modeling accuracy compare to that of other similar studies?As stated in the results section, the FoM was calculated in our study as 24.37%.If we look at other similar studies for comparison, the FoM values range from 1% to 59% depending on the spatial resolution, spatial extent, and number of LULC classes mapped [40].To compare with our study, the Twin Cities or Detroit case [40] would be the most appropriate, as these studies used maps with a similar spatial resolution, number of pixels, and number of LULC classes to our study.These two studies had FoM values of 11% and 15% respectively; slightly lower than our FoM, even though we included one additional LULC class.The studies with much higher FoM values typically had access to site-specific driver variable information (local land allocation plans) and/or aggregated modeling results to coarser spatial resolutions (i.e., larger pixel sizes) to reduce the number of locational errors [40].Although the LULC models implemented differ among studies and have different reference characteristics, this kind of comparison can give a general idea of the modeling accuracy of our study.
Comparing our study to others that used the same MLP Markov Chain modeling method that we used [32,[34][35][36][37][38]72], we can find a few important aspects for discussion.One is that many previous studies did not provide reasoning for why certain driver variables were considered or discuss why certain factors were heavily influencing different types of LULC conversions.Although the number and types of driver variables investigated differ among studies, the level of importance of each variable and the reasons why the variable is affecting LULC changes also likely varies due to different social demands and issues in each study area.By understanding the factors driving LULC conversions, local governments can have a better idea of the issues that region that is facing.We addressed this in our study by first hypothesizing why each variable may be relevant and then by measuring the influence of each variable and discussing why certain variables were relevant for different types of LULC transitions.A second point that we would like to discuss is related to model behavior.Olmedo et al. [36] showed an important aspect of changes, which could not be simulated due to the acceleration of changes that occurred in the reference year, which did not show during the calibration time period.The MLP method is a stationary model, so it is determined that it will keep the same transition rate.This means, if the development scenario changes in the future, it could either increase or decrease its changing rates, resulting in an inaccurate projection.This issue is one reason why we simulated three scenarios of the future LULC in our study area.However, most other works only simulate one future scenario.Due to this model behavior issue and also to provide useful options to government agencies in the area, we recommend considering at least a few different scenarios in addition to the business-as-usual scenario.A third point we would like to emphasize is that there should be some general consensus for presenting model validation results.For instance, one issue with Kappa and other traditional LULC accuracy metrics is that they typically give high accuracy values in study areas with few LULC changes and much lower accuracy values in areas with many LULC changes [37,38].On the other hand, because the FoM is a ratio metric, it is unaffected by the quantity of LULC changes, so it is a good practice to also report FoM in LULC change modeling studies (as we have done in this study).

Accomplished Tasks and Future Works
This work follows similar processes to other MLP Markov Chain studies to model the future landscape of Laguna de Bay.Using Landsat and SAR images, a spatial resolution of 25 m was achieved, showing finer information compared to using Landsat data alone.Developing LULC maps in tropical regions such as the Philippines, Indonesia, or other Southeast Asian countries is a difficult task due to frequent cloud cover (which obstructs the view of the land surface for optical satellite images).Therefore utilizing SAR data was an advantage in developing the historical LULC change data.We have worked on finding the characteristics of the changing pattern and compared these with other related studies.The 27 variables we considered was an enormous amount, around double the number of variables from the study that showed the highest number (14 drivers [72]).We have found explainable relations with the drivers used, and the study gave a clear idea of how LULC is likely developing in this area.Although we did use a large number of driver variables, we could not take into account the local zoning policies of the cities and towns in the study area (due to a lack of data availability), which also have a significant effect on LULC change.It is a very time-consuming and challenging task to visit all of the municipalities in the area to collect (and possibly digitize) their zoning and development plans, but this would probably further enhance the modeling results.
We found that urban sprawl, which was the focus of this study, is expected to continue occurring throughout the future timeframe that we considered (2016-2030).Looking at the surrounding area as a continuous landscape, we were concerned especially with how the existing agricultural and forested lands could be maintained in the future.We plan to continue the work to present a scientific standard for how these lands should be preserved according to various social-environmental impacts, which could be caused by the future LULC changes.

Conclusions
The objectives of this study were to model and map the future LULC of the Laguna de Bay area of the Philippines, which has a significant impact on the lake's water quantity and water quality.The study area is significant mainly because of the importance of the lake to the population of Metro Manila and the surrounding cities and towns.The future LULC was modeled and mapped by using Markov Chain analysis, and the transition probabilities were calculated by using historical LULC change data and freely available data related to the drivers of the LULC changes such as demographic, economic, topographic, infrastructure, and climate variables.Three questions are addressed in this work: (1) Where are different LULC changes taking place?; (2) What are the variables that best explain the changes (e.g., what are the drivers)?; and (3) What is the rate at which these changes are likely to take place?The major LULC changes in the area included an increase in Built-Up areas at the west side of the Laguna de Bay, south of Metro Manila, and changes of many areas between the Crop-Grass and Trees classes, owing to logging, cropland abandonment, reforestation, and other factors.In total, approximately 7800 to 44,000 ha of land within the study area are modeled to be converted to the Built-Up class by 2030, depending on the development scenario.We tested three scenarios: 'business-as-usual', 'compact development', and 'high urban sprawl'.The study has shown the extent of the changed areas in addition to the patterns and the locations of these future changes and identified the variables considered to be significant drivers of these changes (which varied for different types of LULC transitions).This information can be used by decision makers in deciding the

Land 2017, 6 , 26 4 of 21 Figure 1 .
Figure 1.Overview of the study area, the Laguna de Bay district, and the surrounding environment in the Philippines.Metro Manila is shown in the northwest corner.Referenced from Google Earth [43].

Figure 1 .
Figure 1.Overview of the study area, the Laguna de Bay district, and the surrounding environment in the Philippines.Metro Manila is shown in the northwest corner.Referenced from Google Earth [43].

Figure 2 .
Figure 2. Overall flowchart of the methodology used in this study.

Figure 3 .
Figure 3. (a) Area of transitions occurring from 2007 to 2015 for each land-use/land-cover (LULC) type; (b) Total area of gains and losses computed within the study area for each LULC type.

Figure 3 .
Figure 3. (a) Area of transitions occurring from 2007 to 2015 for each land-use/land-cover (LULC) type; (b) Total area of gains and losses computed within the study area for each LULC type.

Figure 4 .
Figure 4. Patterns or trends of transition through space for (a) Crop-Grass to Built-up; (b) Crop-Grass to Trees; (c) Trees to Built-up; and (d) Trees to Crop-Grass.Red (blue) color indicates the higher (lower) chances of transition.The data contain no specific values because no mean is represented.

Figure 4 .
Figure 4. Patterns or trends of transition through space for (a) Crop-Grass to Built-up; (b) Crop-Grass to Trees; (c) Trees to Built-up; and (d) Trees to Crop-Grass.Red (blue) color indicates the higher (lower) chances of transition.The data contain no specific values because no mean is represented. Table

4
(a) to (d) show the accuracy and rankings of each variable that has the largest effect on the skill of the model for (a) Crop-Grass to Built-Up; (b) Crop-Grass to Trees; (c) Trees to Built-Up; and (d) Trees to Crop-Grass.The accuracy is based on the results of the 5000 testing pixels.The most (least) influential variable for model (a) was NL_2015 (BIO12).

Figure 5 .
Figure 5. Land-use/land-cover (LULC) modeling of the Laguna de Bay environment for (a) 2015; (b) 2023; and (c) 2030.(a-c) shows the hard classified map; (d-f) shows the percentage of each LULC class within the study area.The numbers indicated after the alphabet represent the different modeled scenarios.Thus (a-1) is the business-as-usual scenario; (a-2) is the compact development scenario; and (a-3) is the high sprawl scenario.Legends in the pie-chart correspond to the colors of the LULC class in the hard classified map.

Figure 5 .
Figure 5. Land-use/land-cover (LULC) modeling of the Laguna de Bay environment for (a) 2015; (b) 2023; and (c) 2030.(a-c) shows the hard classified map; (d-f) shows the percentage of each LULC class within the study area.The numbers indicated after the alphabet represent the different modeled scenarios.Thus (a-1) is the business-as-usual scenario; (a-2) is the compact development scenario; and (a-3) is the high sprawl scenario.Legends in the pie-chart correspond to the colors of the LULC class in the hard classified map.

Figure 6 Figure 6 .
Figure 6 shows the accumulated increase in land area of each transition class from the starting year of 2015 to 2030 in yearly increments for the business-as-usual scenario, compact development scenario, and high sprawl scenario.In the business-as-usual scenario, most of the LULC change classes showed increases in total area and different growing rates; however, the Crop-Grass to Trees class showed a decrease in area beginning in 2026.The Trees to Crop-Grass and Crop-Grass to Trees classes showed a polynomial-like increase, resulting in increases of 26,409 ha and 21,166 ha, respectively, in 2030.The Crop-Grass to Built-Up and Trees to Built-Up classes showed a linear increase of land area, resulting in increases of 14,137 ha and 3946 ha, respectively, in 2030.These changes are limited to the study area.More meaningful values might be extracted when areas are

Figure 6 .
Figure 6.Trends of LULC class and net area of increase at each stage of the modeling.The base year is 2015.C2B, Crop-Grass to Built-up; T2B, Trees to Built-up; T2C, Trees to Crop-Grass; C2T, Crop-Grass to Trees.(a) Business-as-usual scenario; (b) compact development scenario; (c) high sprawl scenario.
Overall flowchart of the methodology used in this study.

Table 1 .
Complete list of all variables collected in this study.

Table 3 .
Test of the explanatory power (Cramer's V) of each variable.

Table 4 .
The sensitivity of the model in maintaining selected inputs.The output shows the accuracy of the case in which all combinations of variables were used except for one to remain constant.Together it shows the ranking of variables from most to least influential given to the models for (a) Crop-Grass to Built-Up; (b) Crop-Grass to Trees; (c) Trees to Built-Up; and (d) Trees to Crop-Grass.