Changes in Gross Primary Production ( GPP ) over the Past Two Decades Due to Land Use Conversion in a Tourism City

Understanding the changes in gross primary production (GPP), which is the total carbon fixation by terrestrial ecosystems through vegetation photosynthesis, due to land use conversion in a tourism city is important for carbon cycle studies. Satellite data from Landsat 5, Landsat 7 and Landsat 8 and meteorological data are used to calculate annual GPP for 1995, 2003 and 2014, respectively, using the vegetation production model (VPM) in the tourism city Denpasar, Bali, Indonesia. Five land use types generated from topographic maps in three different years over the past two decades are used to quantify the impacts of land use changes on GPP estimation values. Analysis was performed for two periods to determine changes in land use and GPP value as well as their speed. The results demonstrated that urban land development, namely, the increase of settlement areas due to tourism activity, had overall negative effects on terrestrial GPP. The total GPP of the whole area decreased by 7793.96 tC year−1 (12.65%) during the study period. The decline is due to the conversion of agriculture and grassland area into settlements, which caused the city to lose half of its ability to uptake carbon through vegetation. However, although forest area is declining, forest maintenance and restoration by making them protection areas has been helpful in preventing a drastic decline in GPP value over the past two decades. This study provides information that is useful for carbon resource management, tourism, policy making and scholars concerned about carbon reduction in a tourism city.


Introduction
Around the globe, urbanization is expected to increase significantly in the coming decades and is projected to continue [1,2].Rapid urbanization and urban sprawl have a significant impact on conditions of urban land use changes due to the increasing demands for the construction of new residential, commercial, utility, and transport infrastructures [3,4].Urbanization is a process related to population increase in urban regions that is associated with settlement increases and automatically reduces vegetation areas due to limited land [5].Population growth and economic changes are two important factors influencing land use changes and the availability of agricultural land, especially in urban areas [6].One of the drivers of economic and population growth is tourism activities [7,8].Tourism has a strong impact on land use, both at particular sites and in adjoining zones and travel corridors [9].In Bali, Indonesia, an increase in the tourism industry and the population is creating huge social and environmental problems, and the main threats to irrigated agriculture and the related livelihoods are land conversion and the transfer of water from agriculture to non-agricultural uses [10].
Most land use changes in Indonesia at the local and regional scales are conversions from agricultural land into residential land and from forest to agriculture [11][12][13][14].Land use change can have a large impact on ecosystems, altering their composition and structure as well as their function, including matter and energy cycles [15].Several studies have explained that land use and land cover changes in urban areas produce global consequences, including effects on the climate, e.g., [16], hydrology, e.g., [17], and biodiversity, e.g., [18], and that these changes alter ecosystem services and affect the biological carbon uptake of urban vegetation, e.g., [19].Land use changes in urban areas play a significant role in determining local, regional, and global carbon emissions [20].Meanwhile, urban ecosystems can also account for a significant portion of terrestrial carbon storage at both local and regional scales [21].
Human activities are highly dependent on ecosystem goods and services, especially net ecosystem CO 2 exchange (NEE), the balance between gross primary production (GPP) and ecosystem respiration [5].GPP is a key measure of carbon mass flux in carbon cycle studies, which is the total carbon fixation by terrestrial ecosystems through vegetation photosynthesis [22].In addition, GPP can be used to quantify vegetation productivity and to understand temporal variability in productivity [23].The quantification of the magnitude of carbon uptake by terrestrial vegetation and how it varies due to land use changes in urban areas is an important question, as future potential sequestration is influenced by both increased atmospheric CO 2 and the changing climate [24], especially for tourism cities.Therefore, the estimation of GPP in an urban area is essential for the quantification of net terrestrial carbon and thus is of particular importance for global carbon cycle research.Previous studies on GPP impacts from land cover/land use changes associated with urban growth have been documented in Southeastern Michigan [19,25].However, changes in GPP due to land use conversion as a result of urbanization in the tourism area have not been reported in previous work.
The remote sensing of vegetation GPP is an important step in analysing terrestrial carbon cycles in response to changing climate [26].Repetitive and systematic satellite remote-sensing observations of vegetation dynamics and ecosystems allow us to characterize vegetation structure, including leaf area index (LAI) and the absorptivity of photosynthetically active radiation (APAR) [27,28].Remote-sensing techniques are advantageous because they allow monitoring of the GPP on a regional and local scale through its relationship with the fraction of absorbed photosynthetically active radiation (fAPAR) driven by vegetation indices [29,30].Moreover, remote sensing is the best choice for urban areas to provide a better understanding of the effects of land use change on vegetation dynamics [31][32][33].Moderate-resolution Imaging Spectroradiometer (MODIS) satellite remote sensing data are commonly used to calculate GPP, e.g., [34].However, for urban regions that have heterogeneous landscapes, Landsat data are more frequently used [19,35].
GPP is generally estimated using models because it is difficult to measure directly [36].Satellite-based, light-use efficiency (LUE) models have been widely used because they rely on simple algorithms to estimate GPP [37].The LUE model states that carbon exchange is a function of the amount of light energy absorbed by vegetation and the efficiency with which that light energy is used to fix carbon [38].One of the LUE models that is commonly used for estimated GPP is the vegetation production model (VPM) [39].The VPM is a light-use efficiency model that utilizes remote sensing imagery to estimate GPP based on the impacts of temperature, water stress, and phenology [40].GPP from Landsat data based on VPM estimates has been validated by several studies that indicate good agreement with CO 2 eddy flux towers [41,42].The potential of using the VPM model to scale up GPP estimation in Indonesia has also been evaluated for CO 2 flux tower sites that use MODIS data [43].
Based on these conditions, in this work, we attempt to use Landsat data and land use information from topographic maps to estimate vegetation carbon uptake by urban tourism areas and their changes over the past two decades due to land use conversion.Three types of Landsat data, including Landsat 5, Landsat 7 and Landsat 8, are used to calculate annual GPP for 1995, 2003 and 2014, respectively.
Annual changes in GPP values estimated using VPM in Denpasar, an urban tourism city on Bali Island, Indonesia, were analysed and extracted for each land use type.It is expected that the results can provide information related to the impact of land use change on ecosystem services and particularly to the biological carbon uptake capability of urban vegetation.A map of Denpasar city, including boundaries before and after reclamation, is presented in Figure 1.Denpasar is located between two neighbourhood regencies, namely, the Badung Regency to the north and west, the Gianyar Regency to the east, and the Badung Strait area to the south.The three districts, including Denpasar, have the highest tourism activity in Bali.In Denpasar, the centre of tourism activity is in the southeastern part, namely, the Sanur region.

Research Location
data, including Landsat 5, Landsat 7 and Landsat 8, are used to calculate annual GPP for 1995, 2003 and 2014, respectively.Annual changes in GPP values estimated using VPM in Denpasar, an urban tourism city on Bali Island, Indonesia, were analysed and extracted for each land use type.It is expected that the results can provide information related to the impact of land use change on ecosystem services and particularly to the biological carbon uptake capability of urban vegetation.

Research Location
Denpasar (Bali, Indonesia) is located at 08°35'31"-08°44'49" south latitude and 115°10′23″-115°16′27″ east longitude.Denpasar city has four districts: North Denpasar, South Denpasar, West Denpasar, and East Denpasar.In 1995, Denpasar had an area of 120.47 km 2 , which increased to 124.68 km 2 because of land reclamation at Serangan Island and Benoa port in 1996, which is located in South Denpasar.A map of Denpasar city, including boundaries before and after reclamation, is presented in Figure 1.Denpasar is located between two neighbourhood regencies, namely, the Badung Regency to the north and west, the Gianyar Regency to the east, and the Badung Strait area to the south.The three districts, including Denpasar, have the highest tourism activity in Bali.In Denpasar, the centre of tourism activity is in the southeastern part, namely, the Sanur region.The population of Denpasar reached 335,196 people in 1992 and increased to 880,600 people in 2015, an increase of more than 250%.Topographically, the study area has a flat relief with zero elevation in the south and southeast, a highest elevation of 89 m a.s.l. in the north, and an average elevation of 23 m a.s.l.In general, the rainfall patterns of the Denpasar area are influenced by monsoons, with the maximum amount of precipitation occurring during the peak of the wet season from December to February and decreasing to a minimum during the valley of the dry season from June to August [44,45].Moreover, the annual average temperature of Denpasar is 27.6 °C, the average annual rainfall is 1803.62 mm year −1 [44], and settlement was the dominant land use in Denpasar city in 2006 followed by rice field areas [46].The population of Denpasar reached 335,196 people in 1992 and increased to 880,600 people in 2015, an increase of more than 250%.Topographically, the study area has a flat relief with zero elevation in the south and southeast, a highest elevation of 89 m a.s.l. in the north, and an average elevation of 23 m a.s.l.In general, the rainfall patterns of the Denpasar area are influenced by monsoons, with the maximum amount of precipitation occurring during the peak of the wet season from December to February and decreasing to a minimum during the valley of the dry season from June to August [44,45].Moreover, the annual average temperature of Denpasar is 27.6 • C, the average annual rainfall is 1803.62 mm year −1 [44], and settlement was the dominant land use in Denpasar city in 2006 followed by rice field areas [46].

Data Use
The main remote sensing data used in this study were level-1 geometrically corrected (L1G) Landsat 5, Landsat 7, and Landsat 8 data acquired on 2 February 1995, 21 March 2003, and 27 March 2014, respectively, which were obtained from the United States Geological Survey (USGS) via the EarthExplorer database (https://earthexplorer.usgs.gov).However, clouds and cloud shadows obstructed the view of the land surface during dataset acquisition in the study area.Clouds and their shadows are one of the main problems when using Landsat data in the tropics [47,48].The brightening effect of the clouds and the darkening effect of cloud shadows create significant sources of noise in the Landsat data and cause problems in cloud detection, which is an initial step in most analyses, including error estimation for vegetation indexes [49].Clouds and their shadows covered 1.77%, 2.87% and 0.13% of the entire area of Denpasar for Landsat 5 on 2 February 1995, Landsat 7 on 21 March 2003, and Landsat 8 on 27 March 2014.Most clouds and cloud shadows cover the agriculture and grass, settlement, and forest land use types.The area of clouds and cloud shadows for the three main Landsat datasets was determined using the automated cloud-cover assessment (ACCA) algorithm [50].The application of the ACCA algorithm requires green, red, near-infrared (NIR), shortwave infrared (SWIR), and thermal infrared (TIR) spectral bands in surface reflectance form.To replace the data lost due to clouds and cloud shadows, Landsat data acquired adjacent to the main data are employed to eliminate the impact of pixels lost in GPP estimation.As a consideration before using these data to replace the lost data, first, the Landsat data must have the date of acquisition closest to the primary data to avoid land surface changes due to agricultural processes.Second, the land surface in the main data and the secondary data used as a patch for lost pixels must have similar spectral characteristics, especially for unchanged land surfaces such as built-up areas.Third, the Landsat data used as the missing pixel cover must be cloud-and cloud shadow-free in the area to be patched.Based on these three considerations, a Landsat 5 acquisition on 26 May 1995, a Landsat 7 acquisition on 28 March 2003, and a Landsat 8 acquisition on 11 March 2014 were used as patches for pixels lost in 1995, 2003 and 2014, respectively.However, only 1.11% of the 2.87% pixels lost in 2003 are covered by secondary data, and the Landsat 7 acquisition on 22 April 2003 was utilized to cover the nearly 1.76% missing pixels that remained.
The spectral similarity between the main data and secondary data is derived by correlating the remote sensing indices that are used to estimate GPP from VPM.The two remote sensing indices are the enhanced vegetation index (EVI) and the land surface water index (LSWI).EVI was proposed based on a feedback-based approach that incorporates both background adjustment and atmospheric resistance concepts into the normalized difference vegetation index (NDVI) [51].The EVI employs an additional blue band for atmospheric correction and variable soil and canopy background reflectance [52] and is calculated from blue, red, and NIR bands.The LSWI, which is proposed by Xiao et al. [53], is the normalized difference between NIR and SWIR spectral bands.Because the SWIR spectral band is sensitive to vegetation water content and soil moisture, a combination of the SWIR and NIR spectral bands has been utilized to derive water-sensitive vegetation indices.For instance, when the leaf water content or soil moisture increases, SWIR absorption increases and SWIR reflectance decreases, resulting in an increase in the LSWI value [40].The EVI and LSWI are calculated as: where ρ blue , ρ red , ρ nir and ρ SWIR are the Landsat (5, 7 and 8) surface reflectance for the blue, red, NIR and SWIR bands, respectively.The relationship between the two indices from the primary and secondary Landsat datasets for the study site represented very good agreement, with correlation (r) values ranging from 0.90 to 0.96 for EVI and from 0.90 to 0.95 for LSWI, indicating that secondary data can be used as a patch.The scatter plots between the indices from the primary and secondary Landsat datasets are presented in Figure 2.
ISPRS Int.J. Geo-Inf.2018, 7, x FOR PEER REVIEW 5 of 21 data can be used as a patch.The scatter plots between the indices from the primary and secondary Landsat datasets are presented in Figure 2. Topographic maps for 2003 and 2014 with a scale of 1:25,000, interpreted by the Japan International Cooperation Agency (JICA) [54] and Regional Planning and Development Board (Bappeda) of Denpasar city, were used to generate land use types for 2003 and 2014, respectively.However, due to the absence of land use data for 1995, the topographic maps in 1992 from the Indonesian Geospatial Information Agency (BIG, previously known as the Indonesian National Coordination for Survey and Mapping, BAKOSURTANAL) were used to produce land use types for 1995.Land use data from 1992 and 2014 were derived from a basic map that is used by stakeholders and policy makers published by the Indonesian government through central and local governments, respectively.Land use data for 2003 were interpreted by JICA based on Landsat data from 2003 that were utilized by the Public Works Department of Bali Province.Through the Indonesian Geospatial Information Agency, the Indonesian government has a regulation that when a map is used as a base map for the government, it should have an accuracy ≥85%.Therefore, all land use data used in this study have an accuracy level greater than or equal to that allowed.All land use data used in this study are in the form of vector data and are subsequently converted into raster data with the same spatial resolution as the Landsat data.To address the complexity of land use types, land use was labelled, and the clusters were combined into five land use types, similar to the classification scheme used by Zhao et al. [25], i.e., settlement, forest, agriculture and grass, others, and water (Table 1).Water and other land use types were not included for analysing changes in GPP values due to land use conversion in the current study as these cover types were assumed to be relatively constant with respect to GPP [25].Topographic maps for 2003 and 2014 with a scale of 1:25,000, interpreted by the Japan International Cooperation Agency (JICA) [54] and Regional Planning and Development Board (Bappeda) of Denpasar city, were used to generate land use types for 2003 and 2014, respectively.However, due to the absence of land use data for 1995, the topographic maps in 1992 from the Indonesian Geospatial Information Agency (BIG, previously known as the Indonesian National Coordination for Survey and Mapping, BAKOSURTANAL) were used to produce land use types for 1995.Land use data from 1992 and 2014 were derived from a basic map that is used by stakeholders and policy makers published by the Indonesian government through central and local governments, respectively.Land use data for 2003 were interpreted by JICA based on Landsat data from 2003 that were utilized by the Public Works Department of Bali Province.Through the Indonesian Geospatial Information Agency, the Indonesian government has a regulation that when a map is used as a base map for the government, it should have an accuracy ≥85%.Therefore, all land use data used in this study have an accuracy level greater than or equal to that allowed.All land use data used in this study are in the form of vector data and are subsequently converted into raster data with the same spatial resolution as the Landsat data.To address the complexity of land use types, land use was labelled, and the clusters were combined into five land use types, similar to the classification scheme used by Zhao et al. [25], i.e., settlement, forest, agriculture and grass, others, and water (Table 1).Water and other land use types were not included for analysing changes in GPP values due to land use conversion in the current study as these cover types were assumed to be relatively constant with respect to GPP [25].

Methods
GPP is estimated by the VPM.The VPM is an LUE model based on the conceptual partitioning of light absorption by chlorophyll pigments (FPAR chl ) and non-photosynthetic vegetation (NPV such as branches, trunks, stems or senescent leaves), where only FPAR chl is used for photosynthesis [39].A brief description of VPM is given below: where ε g is light-use efficiency (µmol m −2 s −1 , photosynthetic photon flux density, PPFD); FPAR chl is the fraction of photosynthetically active radiation (PAR) absorbed by chlorophyll; and PAR is the photosynthetically active radiation (µmol photosynthetic photon flux density, PPFD) calculated as 0.45 × R (R, incoming global solar radiation).In the current study, daily PAR is summed in a year to obtain an annual PAR value; thus, annual GPP can be obtained.Light-use efficiency (ε g ) is a function of the maximum light-use efficiency (ε 0 ), temperature (T scalar ), water condition (W scalar ), and leaf phenology (P scalar ) and is defined as: A maximum light-use efficiency, ε 0 , value of 0.40 gC mol −1 PPFD [55] was used in this study.The parameter T scalar represents the temperature sensitivity of photosynthesis and is calculated at each time step using the equation developed for the terrestrial ecosystem model [56]: where T is the monthly average air temperature (when PAR > 1 µmol m −2 s −1 ) obtained from the Indonesian Meteorology, Climatology and Geophysics Agency (BMKG) and the parameters for T min , T max , and T opt for photosynthesis were 0, 48, and 28 • C, respectively, based on Vetrita et al. [43].
The W scalar takes into account the complex impact of water stress on photosynthesis (i.e., changes in stomatal conductance, leaf water potential, etc.) caused by soil moisture and/or atmospheric water deficits [57].W scalar is calculated as follows [39]: where LSWI max is the maximum LSWI during the growing season for each pixel.
In the VPM, P scalar is included to account for the effect of leaf age on photosynthesis at the canopy level [58].The calculation of P scalar is dependent upon the longevity of leaves (deciduous, versus evergreen).For a canopy dominated by leaves with a life expectancy of one year, P scalar is calculated in two different phases as a linear function [39].Because urban areas have very complex and heterogeneous plant structures and always have new leaves emerging during the plant-growing season, P scalar is set to be 1.0 in this study.
The FPAR chl within the photosynthetically active period in the VPM is assumed to be a linear function of EVI, and the coefficient a is set to be 1.0 [39]: PAR irradiance measurements were not routinely carried out at radiometric sites, such as standard actinometric stations, where broadband solar irradiance components are routinely measured [59].PAR is the region of incoming global solar radiation reaching the earth surface ranging from wavelengths of 400 to 700 nanometers [60].Therefore, PAR is assumed to be approximately half of the global solar radiation [61].Based on the literature compiled by Escobedo et al. [62], the ratio between PAR and global solar radiation ranges from 0.42 to 0.52.For this study, a ratio of 0.45 is used to calculate PAR from global solar radiation because this factor is employed in most studies related to GPP estimation, e.g., [63].Global solar radiation was calculated using the modified Sayigh universal formula for the Indonesian region [64] by utilizing daily data on bright sunshine hours received during the day, maximum air temperature, and relative humidity obtained from the BMKG in a year to obtain annual PAR for 1995, 2003 and 2014.Details on the global solar radiation computation method using the modified Sayigh universal formula can be found in Halawa and Sugiyatno [64].
In this study, the remote sensing indices EVI and LSWI are assumed to be the same throughout the year.Consequently, the monthly and seasonal variability of vegetation cover and water stress cannot be known throughout the year.Only PAR and T scalar have monthly and seasonal variability, because they are calculated from daily meteorological data.However, the primary remote sensing Landsat data used in this study were taken during the end of March, which is the beginning of the transition season (March-April-May) between the northwest monsoon (rainy season) and the southeast monsoon (dry season) for the study area [45,65], except for the data for 1995, which were acquired in February.During the transition season, the water content is not too high because rainfall begins to decrease and is not too low because it is not the dry season.In Denpasar, the vegetation in settlement areas is dominated by annual and perennial plants, such as roadside trees, shade trees, and annual flowering plants.The forest land use type is also dominated by perennial plants and mangroves.Nuarsa et al. [66] showed that annual, perennial and mangrove plants in the tropics have relatively stable vegetation indices throughout the year and are unaffected by the monsoon seasons.The agriculture land use type exhibits a seasonal pattern that is influenced by monsoons due to the planting period in which the dominant plant are paddy crops [67].However, Jayanti et al. [68] found that paddy crops in Bali, Indonesia exhibit peaked vegetation index values in late June and early December (cropping period (reproductive stage)) with the lowest values in mid-September (fallow and flooding/transplanting period) and are at an average value in March and late October (cropping period (vegetative growth stage)).Moreover, As-syakur et al. [69] also showed that the vegetation health index (VHI), which can describe variations in vegetation indexes and water stress other than drought conditions, does not have high fluctuations throughout the year in the tropics except during El Niño events.Therefore, it can be assumed that the vegetation cover and water content in the transition season can represent the average annual condition for the study area.
To understand the impact of land use changes on GPP estimation in an urban tourism area of Denpasar, GPP values are extracted for each type of land use during the study year.The definition of land use types is based on Table 1.Previously, the area of each land use type and the speed of change were calculated.The total value of GPP is obtained from the sum of the GPP values in each pixel multiplied by the pixel area in all research areas or in each land use type, while the average GPP value is obtained by dividing the total values of GPP by the number of image pixels.In this study, the speed of GPP changes was also calculated in all research areas and in each land use type.
The analysis of land use change was performed for two periods (1992-2003 and 2003-2014) with the topographic data sources previously described.The land use change can be measured using Equation (8) in hectares and Equation (9) in percent, and the equations are defined as follows: where LUC a (ha) is changes in land use type a in hectares from time t (current year) to time t−1 (former year) and LUC a (%) is land use change in percent.LU a,t and LU a,t−1 are the total land area of land use type a in hectares at time t (current year) and time t−1 (former year).In the current study, land use expansion in each type was also calculated and can be measured using a simple equation of the land use expansion index (LUEI).The LUEI can reflect the temporal patterns of expansion or reduction in each type of land use and is defined as follows: where LUEI a (ha year −1 ) is the annual expansion index for land use type a in hectares and can be converted to percent per year (% year −1 ) by multiplying by 100 and can also be converted to hectares per day (ha day −1 ) after dividing by the number of days in the year (356).N is the total number of years from time t (current year) to time t−1 (former year).A positive LUEI value indicates expansion, and a negative value indicates reduction.

Land Use Changes
Table 2 lists the land use area and proportion of the total area in each land cover type in 1992, 2003 and 2014.In 1992, agriculture and grassland comprised the largest proportion of the total area, approximately 65.31% of the total area, while the water and other land use types covered less than 5% of the study area by 1992.Figure 3a shows that agriculture and grassland were located mainly in the areas surrounding settlements, except in the southeastern and southern parts of Denpasar.Interestingly, during 2003, the proportion of agriculture and grassland began to equalize with that of settlements, with settlement areas and agriculture and grass areas accounting for approximately 43.49% and 42.03% of the total area, respectively.During this period, the settlement land use type began to spread in all directions from the city centre and clustered in the southeast region.It is noteworthy that compared to the spatial pattern in 1992, agriculture and grassland became more fragmented and scattered between settlements in 2003 (Figure 3b), and a large area of forest began to establish in the southern part of Denpasar.The forest in the southern part of Denpasar is dominated by mangroves, which began to be replanted in the pond area in 1995, although the project itself began in 1992 [70]. Settlement became the most dominant land cover type in Denpasar city during 2014, accounting for 58.57% of the study area by 2014, followed by agriculture and grassland, accounting for 28.21% and 9.68% of the study area, respectively.The two other types of land use covered less than 5% of the total area.The area of settlement expanded broadly and resulted in clear visibility of the agriculture and grass land use type, especially in the southern part of the study area, which is adjacent to the centre of tourism activity (Figure 3c).Tourism centres are not only located in the southeastern part of Denpasar but also in another district neighbouring the southwestern side of Denpasar, namely, the Kuta tourism centre.    2 and 3).The fastest change in land use in Denpasar was found from 1992 to 2003, especially for settlement and agriculture and grassland.Based on the LUEI value, the settlement area increased by approximately 0.69 ha day −1 , and the agriculture and grass land cover type lost an area of approximately 0.60 ha day −1 .Meanwhile, forests decreased by approximately 0.04 ha day −1 during the first nine years.During the period 2003 to 2014, settlement still showed an expanding pattern with an expansion rate of 0.43 ha day −1 , and agriculture and grass declined with an LUEI rate of 0.39 ha day −1 .In general, the rate of change in the settlement areas in Denpasar over the past two decades increased by approximately 0.67 ha day −1 , and the forest and the agriculture and grass areas decreased by approximately 0.02 ha day −1 and 0.60 ha day −1 , respectively.Land use has changed in the tourism area of Denpasar from 1992 to 2014.Based on the area in 1992, settlement increased by approximately 204.75% in 2014, and the forest and agriculture and grass land cover types decreased to 12.68% and 55.29% ha in 2014, respectively (Tables 2 and 3).The fastest change in land use in Denpasar was found from 1992 to 2003, especially for settlement and agriculture and grassland.Based on the LUEI value, the settlement area increased by approximately 0.69 ha day −1 , and the agriculture and grass land cover type lost an area of approximately 0.60 ha day −1 .Meanwhile, forests decreased by approximately 0.04 ha day −1 during the first nine years.During the period 2003 to 2014, settlement still showed an expanding pattern with an expansion rate of 0.43 ha day −1 , and agriculture and grass declined with an LUEI rate of 0.39 ha day −1 .In general, the rate of change in the settlement areas in Denpasar over the past two decades increased by approximately 0.67 ha day −1 , and the forest and the agriculture and grass areas decreased by approximately 0.02 ha day −1 and 0.60 ha day −1 , respectively.The characteristics of the settlement land use type in Denpasar became increasingly complex from 1992 to 2014, with the geographical variations characterized mainly by the transformation from greenspaces into buildings for residence and government activity as well as into multi-functional buildings that provide general tourism services, including accommodation, shopping, restaurants, and entertainment.The changes in land use in the Denpasar region not only show expanding settlement areas but also expanding forestland due to the conversion of water and other land use types in the southern part of Denpasar and reclamation areas in Serangan Island.The clear land use change in this tourism city is consistent with the studies of Kytzia et al. [71], Xi et al. [72], and Mao et al. [73].However, the process of greenspaces changing into settlement areas in Denpasar is very fast compared to that in other cities in Indonesia and the world.In Jakarta, the capital city of Indonesia, the settlement area has grown to approximately 0.02 ha day −1 over the past four decades (1970 to 2008) [13], and in Beijing, China, the urban land cover from 1986 to 2001 expanded at a rate of 0.06 ha day −1 [74].
Denpasar has been the capital of Bali Province since 1958.Therefore, all activities of the provincial government and activities related to education and business are centred in Denpasar.Although Denpasar City is not a major tourist destination in Bali, being the centre of government and the presence of tourism destinations in the Sanur region have led to growth in tourism, which can have a direct impact on land use change.As shown in Figure 3a,b, most of the forest in southeastern Denpasar, which is the Sanur region, has was converted to settlement from 1992 to 2003, and most of the change was from forest to multi-functional buildings for tourism activities.Bali has been a tourism destination since before Indonesia became independent [75].Tourism in Bali has continued to grow, especially after the construction of the Bali Beach Hotel (now known as the Inna Grand Bali Beach Hotel), the first five-star accommodation for tourism located in the Sanur region of Denpasar.This hotel has been in operation since 1966, has a height of 10 floors, which was the tallest building at the time of its construction and is still the tallest in Bali today.For example, over 11,278 foreign tourists per year visited the island in 1969, and by 2014, the number had grown to 3,766,638.From 1990 to 2000, foreign tourist visits to Bali grew an average of 17.08% per year, and from 2000 to 2014, the average increase was only 11.11% per year.The rapid growth of foreign tourists visiting the island before the new millennium may have led to rapid changes in land use in Denpasar from 1992 to 2003 by increasing the number of multi-functional buildings for tourism services and the incomes of employees in the tourism industry.In addition, the change in land use resulting from the indirect impact of the development of the tourism industry on population growth through migration in the southern part of Bali, including Denpasar, has increased the need for residential buildings.Denpasar City is the most populous area in Bali Province, and in 1992, the population density reached 2704 people per km 2 and increased to 6622 people per km 2 in 2014.Bali has experienced the most rapid economic development of all the islands of Indonesia [76].One of the main reasons for this growth is the many income potential created directly by tourism and indirectly by industries that produce goods for the tourism industry.Consequently, tourism activities can ensure income increases and improvements to the standard of living.In total, 51 percent of people's income and 38 percent of Bali's employment opportunities are directly linked with expenditures by tourists and tourism investment [77].On average, the tourism sector has accounted for 30% of the gross domestic product (GDP) of Bali Province over the past two decades, where economic growth is largely contributed by the southern regions as centres of tourism activity, including Denpasar city [14].The growth in economic income, living standards, and population are common causes of land use change in urban areas [78,79].The relationship between population growth and regional GDP and land use change was also described by Wu et al. [80] in the Yangtze River Basin, China.Moreover, an increase in the number of tourists can affect the land needs associated with their activities [72], and similar results were presented by Allen et al. [81] in South Carolina.
As previously described, the Inna Grand Bali Beach Hotel is the tallest building in Bali today with a height of 10 floors.In Bali, the Decree of the Governor of Bali Province No. 13/Perbang.1614/II/a/1971, which was made in 1971, sets the maximum building height to be 15 metres, typically only four floors.This regulation stems from local wisdom of the Hindu religion and Balinese culture on the urban and tourism development processes.This rule ensures that the construction of buildings such as residences, government offices, businesses, multi-functional buildings for tourism, and others do not spread skyward.Consequently, sideways expansion is a logical choice for obtaining large buildings with the largest volumes but directly and indirectly sacrifices greenspaces for settlement areas.This regulation may be one of the causes of the rapid change in land use in Denpasar, although government regulation can control and block the acceleration of land use change [82][83][84].As mentioned earlier, Jakarta and Beijing did not experience any rapid land use change compared with Denpasar.Although both cities are capital cities, they allow the existence of taller buildings.Therefore, further studies related to the impacts of this regulation need to be undertaken to identify their linkage with declining greenspaces and the expansion of settlement.

Annual Changes in GPP
The VPM was used to estimate GPP in 1995, 2003 and 2014 in the tourism city of Denpasar.Table 4 reports the annual total GPP estimates for each of the land cover types and their percentage of the total.Meanwhile, the spatial distributions of annual GPP in Denpasar as analysed from the Landsat satellite data by the VPM method for 1995, 2003 and 2014 are presented in Figure 4.Because the current study only studied changes in GPP, the study only shows GPP values for land use types related to greenspaces and settlement and does not include GPP values for non-greenspaces and settlement land types.The annual GPP for the entire region in the years 1995, 2003 and 2014 was estimated to be 61,609.15,55,621.81and 53,815.19tC year −1 , respectively.In 1995, the estimated GPP in agriculture and grassland accounted for the largest portion, approximately 75.37% of the total estimated GPP, followed by forest with a proportion of 13.10%, and settlement accounted for 11.53%.It can be clearly seen in Figure 4a that the GPP was estimated to be higher than 550 gC m −2 year −1 in the agriculture and grass land use type.In 2003, the allocation of agriculture and grassland still accounted for the largest area, although it is lower than that in 1995, 55.33% of the total area, while the forest and settlement land use types exhibited increasing allocations of estimated GPP with proportions of 16.61% and 28.06% of the study area by the year 2003, respectively.There is an interesting spatial pattern of estimated GPP in 2003.Inside the area around the settlement, there is still a fairly high GPP value (>550 gC m −2 year −1 ), and the south part of Denpasar has been dominated by high GPP values as a result of the replanting of mangroves.However, the high estimated GPP regions in the southeastern part of Denpasar in 1995 became more fragmented and scattered with low values of GPP in 2003 (Figure 4b), while the same location in 2003 was dominated by multi-functional buildings that serve tourism activities.Interestingly, the GPP estimation for settled areas has overtaken that of agriculture and grass areas, with a value of 41.39% of the total GPP estimation in 2014.The GPP estimation from forest increased by 19.08%, while that of agriculture and grass areas also declined, although the proportion was still quite large with a share of 39.53% of the total area in 2014.It is noteworthy that compared to their spatial pattern in 2003, the medium to high value GPP values were more fragmented and scattered around the regions with low GPP estimates in 2014.High estimated GPP was clearly observed in the southern part of Denpasar and in small areas of the agriculture and grass land use type (Figure 4c).     .However, in general, the mean GPP estimated for the settlement and in agriculture and grassland land cover types was almost similar, and that for the forest increased over the past two decades.This study found that the estimated GPP for the settlement region is half of that for agriculture and grassland and that a change of 1 ha of agriculture and grassland into settlement will cause Denpasar city to lose 3 tC ha −1 year −1 of carbon uptake carbon by green vegetation (GPP).The current study found that the total value of estimated GPP for the settlement land use type increased dramatically.From the lowest proportion in 1995, this land use type achieved the highest ability to uptake total carbon through vegetation of all land use types in 2014.The average value of estimated GPP is similar in the three observation periods, although the average value is lower than the result obtained by Zhao et al. [19] (390.24gC m −2 year −1 ), which may be due to the different methods used in GPP calculation.The increase in the GPP estimation value is very high but is still less than the increase in the settlement area, where in general, there is a linear relationship between the increase in area and the total GPP in the settlement area (see Tables 3 and 5).The linear relationship between changes in settlement area and primary production is consistent with the study of Lu et al. [5].A linear relationship was also found in the agriculture and grass land use type, with decreases in area of approximately 55.29%, followed by a loss of ability to uptake total carbon of 54.19%.Although only approximately half of the agriculture and grass land use type has been converted, the area reached up to 4350.51 ha.The settlement areas increased to 4906.44 ha during the study period, which indicates that agriculture and grass have been converted to settlement areas.Consequently, the ability of cities to uptake carbon from the atmosphere is reduced, as displayed in Table 5. Settlement areas can uptake carbon due to the presence of vegetation.As-syakur et al. [85] have shown that not all settlement areas in Denpasar are covered by structures, but have various built-up land uses that include vegetated and non-vegetated areas such as backyards and front yards.However, this city is unique in that the majority of the population is Hindu, and the city is a "flat city" due to the absence of tall buildings.These characteristics are not found in other regions in Indonesia or in other regions in Southeast Asia.The The region-wide net decrease in GPP was approximately 7793.96 tC year −1 during the past two decades, a decrease of 12.65% from the GPP value in 1995 (Table 5).The estimated total GPP decreased during the 20-year period with an average decline of 0.63% per year.However, the fastest decline occurred in the first nine years compared to 2003 to 2014 (Table 6), which is similar to the decrease in the agriculture and grassland use type.Interestingly, based on land use data, GPP estimates increased for settlements and forests, although the forest area decreased, while a decline in the value of estimated GPP occurs only in the agriculture and grassland use type.The GPP estimate for settlements increased at a rate of 13.31% per year, while during 2003 to 2014, the increase in the GPP estimate was less than that for the first nine years, and over the past 20 years, the increase has been large, reaching 10.68% per year.The estimated GPP in forests over the past two decades has continued to grow, even though the forest area has declined, with the fastest growth occurring from 1995 to 2003.Decreasing agricultural and grass land use has also led to a decline in the value of GPP estimates, and the decline in the past 20 years reached more than 50%.The fastest change, a value of 3.75% per year, was found from 1995 to 2003, and from 2003 to 2014, the value of estimated GPP loss was approximately 2.57% per year.The current study found that the total value of estimated GPP for the settlement land use type increased dramatically.From the lowest proportion in 1995, this land use type achieved the highest ability to uptake total carbon through vegetation of all land use types in 2014.The average value of estimated GPP is similar in the three observation periods, although the average value is lower than the result obtained by Zhao et al. [19] (390.24gC m −2 year −1 ), which may be due to the different methods used in GPP calculation.The increase in the GPP estimation value is very high but is still less than the increase in the settlement area, where in general, there is a linear relationship between the increase in area and the total GPP in the settlement area (see Tables 3 and 5).The linear relationship between changes in settlement area and primary production is consistent with the study of Lu et al. [5].A linear relationship was also found in the agriculture and grass land use type, with decreases in area of approximately 55.29%, followed by a loss of ability to uptake total carbon of 54.19%.Although only approximately half of the agriculture and grass land use type has been converted, the area reached up to 4350.51 ha.The settlement areas increased to 4906.44 ha during the study period, which indicates that agriculture and grass have been converted to settlement areas.Consequently, the ability of cities to uptake carbon from the atmosphere is reduced, as displayed in Table 5. Settlement areas can uptake carbon due to the presence of vegetation.As-syakur et al. [85] have shown that not all settlement areas in Denpasar are covered by structures, but have various built-up land uses that include vegetated and non-vegetated areas such as backyards and front yards.However, this city is unique in that the majority of the population is Hindu, and the city is a "flat city" due to the absence of tall buildings.These characteristics are not found in other regions in Indonesia or in other regions in Southeast Asia.The city of Denpasar contains a "holy area" for Hindu religious activities at each residence that is characterized by the presence of vegetation with varying levels of cover and heights.However, the increase in settlement area has not impacted the average GPP estimates, unlike the results obtained by Zhao et al. [19], which show an increase in the mean annual GPP at the time of an increase in the built-up area in the eastern United States.Given that the development of settlement areas in the city is difficult to control due to high demand for tourism activities, an increased settlement area accompanied by good management of "holy areas", backyards, and front yards should be able to increase the ability of settlement areas to uptake atmosphere carbon through vegetation.
Interestingly, although the forest areas are degraded, forest restoration has been helpful in preventing a drastic decline in GPP in Denpasar city over the past two decades.A forest area of approximately 175.41 ha was lost during the study period, but the annual estimated GPP increased by approximately 27.23% from 1995 to 2014.This increase is due to the restoration of mangrove forests in the southern part of Denpasar from 1995 to the present.Reforestation is conducted by replanting mangroves in water (pond) land use areas and maintaining these areas by making them protection areas.The mangrove forest in southern Denpasar is part of a protected area established by the Indonesian government, the Ngurah Rai Grand Forest Park, through decision No. 544/Kpts-II/1993 in 1993.
The current study indicated the importance of forest (mangrove) restoration to sustain the carbon cycle in urban areas that have experienced rapid land use change due to increased tourism activity and urbanization.Ecological restoration is important to increasing the total carbon uptake by vegetation.Successful ecological restoration associated with primary production has been reported previously by Yang et al. [86], but this result was an increase in the amount of carbon sequestration due to an increased terrestrial vegetation area, not under the decline of forest area.However, Briber et al. [87] found an increase in tree productivity during forest conversion to urban areas.Forest productivity is a function of many factors such as species, age, resource availability, growing season length, and competition [88,89].Mangroves are among the most carbon-rich forests in the tropics [90].Duarte et al. [91] also suggest that mangrove carbon production is more rapid than other estuarine and marine primary producers.On the other hand, Landsat data are some of the primary sources for identifying areas of forests in space and time with good accuracy [92], including in urban areas [93] and mangroves [94].
However, this study still has some uncertainties.The land use data used in this study do not match the year, especially for GPP analysis in 1995, which may cause a mismatch in extracting the estimated GPP.This mismatch will lead to an underestimation of GPP values in settlements and an overestimation in agriculture and grass areas because most of the settlement areas come from agriculture and grass areas.From 1992 to 1995, Denpasar had an average population growth rate of 2.85% per year, and foreign tourist visits to Bali grew an average of 13.47% per year in the same period.One consequence is the increasing number of multi-functional buildings for tourism services and the number of residences needed to accommodate the population growth.However, the four-year difference may not have a significant impact on a study that had a long span period.In addition, the use of 30-m Landsat data can cause errors in detecting vegetation in heterogeneous settlement areas.As shown by As-syakur et al. [46], different spatial resolutions in the satellite data used for analyses of GPP in urban areas will lead to different averages for settlement areas, and more detailed spatial resolutions will give larger average values.Increased pixel sizes (or decreased spatial resolution) result in the loss of image detail [95].Satellite data with high spatial resolutions may be able to reduce the difficulties with remotely sensed data from coarse-resolution satellites [96].Finally, this is the first study that used the VPM method with medium-spatial-resolution satellite imagery to calculate GPP in urban areas; therefore, the results are questionable.Previously, Ciu et al. [97] tried to compare GPP-based VPM estimates (using MODIS data) and MODIS GPP data products (MOD17A2) with solar-induced chlorophyll fluorescence (SIF) over the most populous megacity area with better results compared to those produced with the MOD17A2 product.On the other hand, the total annual GPP estimated using the VPM in the current study was not very different from the GPP estimate produced by other LUE models or from the MODIS GPP product, as presented by As-syakur et al. [46] in the same location.This study can be used as an initial source of information related to land use changes and their impact on terrestrial carbon uptake by vegetation.However, urban areas are highly heterogeneous landscapes that have different or even opposing effects on overall urban vegetation productivity [98].Urban areas also have complex climate systems that are affected by complex socio-ecological systems [99] and the urban-rural proportion as related to urban heat islands [100], which directly and indirectly affect the calculation of light-use efficiency.Therefore, a validation of the estimated GPP results of the VPM method with eddy flux towers needs to be performed to advance our quantitative understanding of the capability of this method to contribute to analyses of vegetation carbon uptake in urban areas.

Conclusions
A recent study indicates that land use has changed in the tourism area of Denpasar over the past two decades and that this change has had overall negative effects on terrestrial GPP.The total GPP of the entire study area decreased by 7793.96tC year −1 during the 20-year period, with an average decline of 0.63% per year.The fastest decline in GPP values was found in the first nine years, especially for the agriculture and grass land use type, and a rapid increase in the settlement area was found at the same time.Although the growth in the settlement area can increase the total GPP value, the amount is not sufficient to cover the total amount of carbon that can be taken up by other types of greenspaces.Based on the mean GPP value for each land use type, a change of 1 ha from agriculture and grass area into a settlement will cause a loss of 3 tC ha −1 year −1 in the ability of the city to uptake carbon through vegetation.Moreover, despite losses in forest area of approximately 12.68%, forest areas will still be able to increase the ability to uptake carbon from the atmosphere as long as the forest is maintained and restored by making them protection areas.
Further studies related to the direct impacts of tourism activity on terrestrial carbon uptake by vegetation need to be undertaken in the future.Moreover, seasonal analysis and validation of the GPP calculated using VPM also requires future research using eddy flux towers to determine the capability of the method to estimate GPP in urban areas.Finally, this study provides information that is useful for carbon resource management, the tourism industry, policy makers and scholars.In this way, the adverse impacts of land use change due to tourism activities in urban areas on ecosystem services can be minimized.
Denpasar (Bali, Indonesia) is located at 08 • 35 31 -08 • 44 49 south latitude and 115 • 10 23 -115 • 16 27 east longitude.Denpasar city has four districts: North Denpasar, South Denpasar, West Denpasar, and East Denpasar.In 1995, Denpasar had an area of 120.47 km 2 , which increased to 124.68 km 2 because of land reclamation at Serangan Island and Benoa port in 1996, which is located in South Denpasar.

Figure 1 .
Figure 1.The research location: Denpasar, Bali Province, Indonesia.Dark grey regions in southern Denpasar are Serangan Island and Benoa port after reclamation, and these areas are used for Denpasar boundaries in 2003 and 2014.The dotted line inside Denpasar is the district boundaries.

Figure 1 .
Figure 1.The research location: Denpasar, Bali Province, Indonesia.Dark grey regions in southern Denpasar are Serangan Island and Benoa port after reclamation, and these areas are used for Denpasar boundaries in 2003 and 2014.The dotted line inside Denpasar is the district boundaries.

Figure 2 .
Figure 2. Scatter plots of the Landsat base index from the main data with the data used for patching the no-data regions caused by clouds and their shadows.Top panels present enhanced vegetation index (EVI) and bottom panels present land surface water index (LWSI).(a,b) for Landsat 5 on 3 February 1995 with 26 May 1995; (c,d) for Landsat 7 on 21 March 2003 with 28 March 2003 (black dots) and for 21 March 2003 with 22 April 2003 (red dots); and (e,f) for Landsat 8 on 27 March 2014 with 11 March 2014.

5 EVI 95 Figure 2 .
Figure 2. Scatter plots of the Landsat base index from the main data with the data used for patching the no-data regions caused by clouds and their shadows.Top panels present enhanced vegetation index (EVI) and bottom panels present land surface water index (LWSI).(a,b) for Landsat 5 on 3 February 1995 with 26 May 1995; (c,d) for Landsat 7 on 21 March 2003 with 28 March 2003 (black dots) and for 21 March 2003 with 22 April 2003 (red dots); and (e,f) for Landsat 8 on 27 March 2014 with 11 March 2014.

Figure 4 .
Figure 4. Annual gross primary production (GPP) distribution in Denpasar from Landsat satellite data analysis.The grey colour indicates that the area is not counted due to land use types such as water and others.(a) 1995; (b) 2003; and (c) 2014.

Figure 4 .
Figure 4. Annual gross primary production (GPP) distribution in Denpasar from Landsat satellite data analysis.The grey colour indicates that the area is not counted due to land use types such as water and others.(a) 1995; (b) 2003; and (c) 2014.

Figure 5
Figure 5 shows the mean estimates of annual GPP as analysed from Landsat satellite data for each land use type for 1995, 2003 and 2014.The mean GPP estimation for settlements decreased slightly (339.23 to 304.99 gC m −2 year −1 ) during the 20-year period and increased marginally for the agriculture and grass land use type (593.53 to 605.65 gC m −2 year −1 ).Meanwhile, the average GPP estimation in forests over the past two decades increased dramatically compared to that for the other two types of land use from 1995 to 2014 (544.97 to 846.40 gC m −2 year −1).However, in general, the mean GPP estimated for the settlement and in agriculture and grassland land cover types was almost similar, and that for the forest increased over the past two decades.This study found that the

21 Figure 5 .
Figure 5. Mean of annual GPP in Denpasar from Landsat satellite data analysis of each land use type.

Figure 5 .
Figure 5. Mean of annual GPP in Denpasar from Landsat satellite data analysis of each land use type.

Table 1 .
Definition of land use types for the current study.

Table 2 .
Area of each type of land use in Denpasar from 1992 to 2014 and their proportion of total area.

Table 2 .
Area of each type of land use in Denpasar from 1992 to 2014 and their proportion of total area.

Table 3 .
Changes in land use types in Denpasar from 1992 to 2014 in percent and the land use expansion index (LUEI) value (% year −1 ).

Table 3 .
Changes in land use types in Denpasar from 1992 to 2014 in percent and the land use expansion index (LUEI) value (% year −1 ).

Table 4 .
Annual gross primary production (GPP) in Denpasar from 1995 to 2014 for each different land use as a total (tC year −1 ) and a percentage of the total.

Table 5 .
GPP estimation changes in each type of land use in Denpasar from 1995 to 2014.

Table 6 .
Changes in GPP estimation per year (% year −1 ) for each type of land use from 1995 to 2014.

Table 5 .
GPP estimation changes in each type of land use in Denpasar from 1995 to 2014.

Table 6 .
Changes in GPP estimation per year (% year −1 ) for each type of land use from 1995 to 2014.