Geospatial Analysis of Horizontal and Vertical Urban Expansion Using Multi-Spatial Resolution Data: A Case Study of Surabaya, Indonesia

Abstract

Urbanization addresses urban expansion, and it leads conversion of the green space into the built-up area. However, previous studies mainly focused on two-dimensional (2D) urban expansion rather than three-dimensional (3D) growth. Here, the purpose of this study is to examine the urban expansion, including built-up and green space for both horizontal and vertical dimensions using geospatial analysis including remote sensing (RS) and Geographic Information System (GIS) in the sub-Central Business District (CBD) area of Surabaya, Indonesia. The medium resolution remote sensing data for both image and Digital Surface Model (DSM) acquired by Advanced Land-Observing Satellite (ALOS) were applied for time-1 (2010). The orthophoto and DSM derived by LiDAR were used for time-2 (2016). We quantified the built-up and green expansions in 2D (area), which were extracted from land use/land cover (LU/LC) by applying hybrid classification. The built-up and green expansions in 3D (volume) were estimated by generating a surface feature model. The spatial configuration of area expansion was investigated using patch metric, while the volume growth was examined using the volume expansion rate. We got three findings. (1) The built-up and green area had expanded about 11.54% and 95.61%, respectively, from 2010 to 2016. The expansion of green area presented in a notable portion, which was mainly contributed by the conversion of bareland to playground or park. However, the expansion of built-up area was less than the volume expansion of 20.6%. It revealed that built-up growth led to vertical rather than horizontal development. (2) The built-up area expansion tended to scatter configuration, whereas, the green area expansion tended to aggregate in a linear pattern. (3) The ratio of built-up volume expansion to green volume expansion showed a mean of 3.7, indicating that the development of built-up and green volume was imbalanced. The built-up growth presented higher than the green growth, mainly in the areas with more vertical building establishment. The pressing need for higher green volume in the study area was identified in several sites located at surrounding artery and toll roads. Overall, our approach can be applied as a reference in monitoring neighborhood environment through greening programs for sustainable urban development.

1. Introduction

It is predicted that urban population will gain 66% of the world’s total population by 2050, which means global urbanization, then continuing to increase at the current level, especially in developing countries [1]. Urbanization is always accompanied by spatial expansion of urban land [2]. In Southeast Asia, the urban expansion rate is approximately 2.8% higher when compared to other urbanized regions of the world [3,4]. Consequently, green space had increasing pressure during the urban expansion and caused many serious ecological problems, psychological well-being, and the health of urban dwellers [5]. Thus, the green space transformed into the built-up area has become one of the leading causes of urban environment destruction [6]. Therefore, it is essential to examine and to monitor the change of green space coverage in spatial complexity of the urban landscape when rapid expansion occurs.

Recently, urban expansion is measured by horizontal and vertical urban growth. A third dimension using urban volume is used as an indicator of vertical development [7,8]. Spatial information on the third dimension is a key indicator of urban expansion in dense urban areas where multiple land uses exist [9]. Attention has been addressed to urban modeling with the integration of horizontal and vertical dimensions, not only from the perspective of the scientific community but also from that of urban stakeholder [10]. The inclusion of the third dimension to investigate the thematic factors of urban systems is very important in monitoring urban development and detecting land use/land cover (LU/LC) changes [11].

Remote sensing (RS) data with various spatial resolution and image analysis techniques have been applied to monitor urban expansion with both horizontal and vertical dimensions, while Geographic Information System (GIS) can identify and investigate the configuration of urban expansion in the regard of land use pattern [12,13,14]. The spatial structures and configurations of land use patches have significant impressions on urban environmental attributes. Thus, numerous landscape metrics have been developed to measure the complexity of the urban landscape patterns [15]. Landscape metrics can be applied over an entire area to reflect the overall situation of the urban expansion pattern and can also be efficiently used to provide specific information of urban expansion pattern in patch or class levels [16]. The patch analysis has been extensively used to characterize, monitor, and assess landscape pattern and composition [17,18], and to study LU/LC changes [19,20].

Several studies have been taken to detect both horizontal (2D) and vertical (3D) urban expansions using RS data. Multi-temporal high resolution images were utilized to investigate changes in individual buildings [21]. 3D building information was extracted from Quickbird images to examine the vertical change of urban morphology [22]. Urban volume was estimated using elevation point data, topographic map and building block outline [9]. The study by [23] used the combination of historical 2D urban change and presented 3D city model acquired from LiDAR data. While the study by [24] used high spatial resolution RS data and a statistical method to estimate the population based on the building height. Those earlier studies employed the very high resolution RS data, however, the availability might be very limited, especially in developing countries. Therefore, our study attempted to use bi-temporal RS data, including image and Digital Surface Model (DSM), having very different spatial resolutions.

Advanced Visible and Near Infrared Radiometer type 2 (AVNIR-2) image and Panchromatic Remote Sensing Instrument for Stereo Mapping (PRISM) DSM acquired by the same platform, the Advanced Land-Observing Satellite (ALOS), were used for analyzing time 1 (2010). Orthophoto and LiDAR DSM by the aerial-LiDAR mapping survey were used for analyzing time 2 (2016). The same time and same platform for acquisition address an advantage in the accuracy enhancement of urban volume [10]. Moreover, our study investigated the green expansion both in the horizontal and vertical dimensions, which were not examined by previous studies. While considering the green space as an environmental indicator, the urban volume, including built-up and green volume, is crucial for the urban development. Green volume delivers a new perspective in quantifying the green quantity essential for urban greening study [10].

We selected sub-Central Business District (CBD) of Surabaya, Indonesia as a case study, since it has experienced rapid growth in population and capital investment. To counterbalance the urban expansion, the government has set up regulation to promote the harmonized development in providing the green environment. “The Regulation of Surabaya City Number 12 Year 2014 Concerning Spatial Plan of Surabaya for 2014–2034” states that at least 30% of city area have to be green [25]. Therefore, the purpose of this study is to examine the urban expansion, including built-up and green space for both horizontal and vertical dimensions using geospatial analysis. We have mainly three objectives: (1) the estimation of the built-up and green expansions in horizontal (area) and vertical (volume) dimensions, (2) the investigation of the structure and configuration of the area expansion using a patch analysis, and (3) the measurement of the volume expansion rate.

2. Methodology

2.1. Study Area and Data Used

Surabaya is located in the northeast of Java Island, and it has been growing as a center for business, commerce, and industry of eastern Indonesia (Figure 1). With a population of 2,765,487 and an area of 326.81 km², Surabaya is divided into 12 Units of Development (UD) [26]. These units were composed to promote balanced and harmonious development in the social and economic perspectives. The study area covers 22.08 km², including UD of A. Yani, which was proposed as a sub-CBD with the development spots sidelong the A. Yani artery road connecting Surabaya to other fringe cities [25]. A. Yani UD consists of three districts where the population increased about 13.4% from 2010 to 2016. The population density was 9438 persons per km² in 2010 and 10,701 persons per km² in 2016 [27].

Table 1. Data used in this study.

Data	Date	Source	Description	Purpose
ALOS image	17 July 2010	[37]	- The spatial resolution of the ALOS image is 10 m with 4 bands including blue (0.42–0.5 μm), green (0.52–0.6 μm), red (0.61–0.69 μm) and near-infrared (0.76–0.89 μm). - The spatial resolution of the ALOS DSM is 5 m with a wavelength of 0.52–0.77 μm [38].	- The ALOS image was used to determine the LU/LC 2010 from which the built-up and green area were extracted. - The ALOS DSM was used to generate the DTM 2010 and the Surface Feature Height (SFH) 2010.
ALOS DSM
LiDAR DSM	Acquired on May–August 2016	[33]	- The LiDAR DSM, LiDAR-derived DTM, and orthophoto which are acquired from aerial-LiDAR survey consist of 104 quads regarding a standard map number. - The orthophoto with RGB color has a spatial resolution of 0.08 m produced by orthorectification using a combination of DSM, DTM and several GCPs from a geodetic GPS [33].	- The LiDAR DSM and LiDAR-derived DTM were used to generate the SFH 2016 which was also used as a reference in the validation of SFH 2010. - The orthophoto was used to derive the LU/LC 2016 from which the built-up and green area were examined.
LiDAR-derived DTM
Orthophoto
Base map	2002	[39]	- A scale of base and building maps is 1:5000 with the horizontal accuracy (RMSe) of 0.5 m, processed by a photogrammetric method for frame photograph with sufficient numbers of GCPs to produce the orthophoto which from the maps were generated. - Geospatial Information Agency Indonesia has verified the accuracy of maps in order to establish NSDI.	- The base and building maps were used as a reference in the geometric correction of ALOS image and DSM also in the classification process. - The building map was applied to define the building footprint in the validation of SFH 2010.
Building map
Land use map	2008	[40]	Scale of 1:10,000	The previous land use map was used to conduct the accuracy assessment of LU/LC map 2010 and 2016.

shows the data used in this study. The accuracy of ALOS DSM has been assessed by several researchers. Using the ground control points (GCPs) as measured by high accurate global positioning system (GPS) as the validation reference, the relative accuracy of ALOS DSM showed a standard deviation of DSM height errors about 4.9 m in the flat terrains including urban area and paddy field, 5.8–7.1 m in the various mixed terrain, including mountains, farms, and cities [28]. Another study combined GCPs derived from a ground-based GPS tracks and the LiDAR DSM as the reference of the assessment. The Root Mean Square error (RMSe) was 4 m in the varied topography between 0 to 900 m, including flat terrain, e.g., paddy farm and manmade constructions at the mountainous terrain [29]. The accuracy was consistent and practically satisfying the target accuracy of 5 m (RMSe) [29]. For the 2016 data, the LiDAR point cloud with the density of 8–12 points/m² was acquired to produce the sub-meter resolution DSM of 0.25 m and Digital Terrain Model (DTM) of 0.4 m, which were acceptable ranging from 0.2–0.5 m for the dense urban area [30,31,32]. Vertical accuracy of the LiDAR data was approximately 0.050 m (RMSe) [33]. Both DSM and DTM of LiDAR data were used for orthorectifying aerial photo to produce the base map with scale 1:1000 with horizontal accuracy (RMSe) of 0.25 m [33], using a base map 2002 as a reference of the coordinate system. The accuracy of LiDAR products has been verified by the Geospatial Information Agency of Indonesia in order to establish National Spatial Data Infrastructure (NSDI) as a part of One Map Policy Program [34,35,36]. Thus, both ALOS DSM and LiDAR products were adequate to generate the urban volume.

In the pre-processing stage, the cloud covers masking and the geometric correction were conducted for ALOS DSM (Figure 2). The radiometric calibration was applied for ALOS image and DSM products [28,41]. Both ALOS image and DSM were projected in a Universal Transverse Mercator (UTM) map projection zone 49S with WGS1984 as a datum. However, the image and DSM position were not appropriately overlaid on the base map and had some displacements with the building map. The geometric correction using 55 GCPs from the base map was performed. The RMSe of 3.68 m, which is equivalent to less than a half pixel, was acceptable [42,43]. The same GCPs were applied for the geometric correction of the ALOS DSM. While for the LiDAR DSM, LiDAR-derived DTM, and orthophoto, the pre-processing was conducted for checking order and number of quads, and mosaicking. Similar to the ALOS data, these data had been in the UTM map projection with WGS 1984 as the datum. Indonesian Geospatial Reference System 2013 (SRGI2013) has been established for operationalizing NSDI. Hence, the Vertical Geospatial Reference System (VGRS) was applied while using the Indonesia geoid (Earth Gravitational Model 2008) which was determined by airborne gravity survey bounded to the national gravity control network [44]. The VGRS as a coordinate reference system was used in the elevation determination of LiDAR DSM and the corresponding LiDAR-derived DTM [33].

2.2. Land Use/Land Cover (LU/LC)

Regarding the limitation of single classification, a hybrid approach combining pixel-based and object-based methods was proposed. The development of the hybrid approach is more varied and robust [45,46]. Previous researches have attempted the integration of different classifiers to enhance the quality and accuracy of a land cover classification [47,48]. The advantageous technique has been applied in various spatial resolution data for the medium as well as fine resolution RS image. For the ALOS image, the previous study by [49] proposed the combination of several classification approaches that generated LU/LC of the complex urban landscape more accurately. Herein, the hybrid approach was applied to classify the ALOS image and the orthophoto in five classes consisting of built-up, road pavement, green space, waterbody, and bareland. The built-up contains all building types, including the commercial, institutional, and single house residential. Road pavement consists of all road type, pavement, and asphalt surfaces such as parking lot. Green space includes all type of green features, including all vegetation and plant, grassland, scrub, playground, and cemetery. Waterbody comprises river, pond, and reservoir pool. Bareland is defined by bare soil or vacant land. Those classes were determined in consideration with the interest of the study and dominant LU/LC; herein, the land use map of 2008 was used as a reference in classification.

2.2.1. LU/LC of ALOS Image 2010

For the ALOS image, the hybrid classification was performed in three levels. In level 1, the pixel-based classification was applied using maximum likelihood classification (MLC). 1012 ground truth points were involved, of which 512 points were used for the selection of training sets. The familiarity with the study area allowed in selecting those ground truth points carefully utilizing a combination of landscape signature, site knowledge, and visual interpretation keys, such as texture, tone, size, shadow, shape, and pattern [10,50]. The roadmap and the land use map of 2008 were used as reference data in the training sites determination. However, the mixed pixels produced “unclassified” objects. Therefore, those classes that could be distinguished clearly were then masked out, and the only spectral-mixed class was further examined using the object-based classification in level 2.

Multiresolution segmentation was applied using the scale of 2, the shape of 0.1, and the compactness of 0.7. Several experiments for investigating the behavior of scale, shape, and compactness parameters had been conducted. Using the scale parameter of 2, we found that the objects were very similar to the original pixel level that can be used to differentiate the textural properties of the LU/LC features. The compactness parameter of 0.7 was selected in the regard of the prominence of generated built-up class, which was more congested at a higher level. The experimental results of the segmentation parameters provide the eminence of generated objects for further classification, which affected the classification accuracy [51,52].

Level 2 in classification was applied using rule-based in accordance with Normalized Difference Water Index (NDWI) to delineate the waterbody, and Normalized Difference Vegetation Index (NDVI) to determine the green space. NDWI can effectively augment the water information in most cases [53] as a dimensionless quantity used as an indicator of the surface wetness [54]. NDVI is most widely used as the indicator to differentiate the vegetation in forests, agricultural land cover, and built environment [17,55]. A threshold of NDWI was applied to distinguish the waterbody and non-waterbody class by selecting several objects. Subsequently, the non-waterbody was classified using NDVI. Two thresholds were set to define the green space and built-up classes; therefore, the pixels without fulfilling the standard were categorized into “other class”, denoting mixed features. Herein, the green features, e.g., plants, trees, grasses, were selected to examine the green space class. To define the exact threshold values, several adjustment methods were employed to get information, including histogram, entropy, attribute similarity, and feature shape. Google earth, land use map 2008, and DSM were used to specify objects.

Level 3 was to classify the “other class” using the same training sample sets of MLC to keep the consistency of the spectral information in texture and pattern. The segmented areas for training set were presented by each ground truth point where the centroids of training samples were located close to the points. Nearest neighbor (NN) classifier was performed to produce three classes, including built-up, road pavement, and bareland. Finally, we merged the results of three levels to create the LU/LC map 2010. The accuracy assessment was then applied to the remaining of 500 ground truth points using post-classification to produce overall accuracy (OA) and kappa coefficient (K).

2.2.2. LU/LC of Orthophoto 2016

Two levels of object-based classification were adopted to classify the orthophoto, including rule-based for level 1 and NN classifier for level 2. Like the ALOS image classification, several experiments in selecting the appropriate parameters of segmentation were also conducted. Finally, the scale of 50, the shape of 0.1, and the compactness of 0.5 were selected for the clarity and similarity of the segmented object with original features. The scale parameter is the key factor of multi-resolution segmentation, which manages the maximum heterogeneity of image segments [56]. The scale of 50 was selected, as it showed the best visualization of the resulted LU/LC types in term of homogeneity and color, and gave the effectiveness of processing time. Texture information, including pattern variability, shape, and size, is considered as important indicator in examining the spatial heterogeneity [57]. Using the frequency of the gray-level changes or color scheme correlation, the texture analysis can examine the relationship between pixels and describe the image details [58]. According to a human texture discrimination study, the texture analysis can be naturally distinguished by second-order statistics, such as entropy, contrast, dissimilarity, homogeneity, and asymmetry [59]. Therefore, these indices were applied in the rule-based classification. Also, the appearance of the red band, which was mostly known sensitive to vegetation, was employed to differentiate the green features.

Threshold selection was a crucial step in extracting water pixels in the first stage. The use of an empirically selected threshold which was combined with the texture analysis of candidate water features could determine the waterbody. To be consistent, we preserved the same water objects that were used in the rule-based classification of ALOS image as the threshold definition. This idea was also applied in the determination of green space and built-up threshold. For the other two classes, road pavement and bareland, the shape histograms of contrast, asymmetry, and compactness were also employed in defining those thresholds. However, the spectral similarity between water, vegetation and other dark objects caused the pepper noise. Therefore, the homogeneous segment and its mean value were computed. Thus, we conducted the texture analysis of candidate features combined with histogram interpretation to eliminate the noise effect. Even though, it was difficult to eliminate the shadows from other dark objects. Therefore, the NN classifier was applied in level 2 for unclassified pixels using 512 ground truth points as the training sets. Three classes were generated, including built-up, road pavement, and bareland. Finally, we merged all of the classes determined by two levels of classification to produce the LU/LC map 2016. The integration of spectral, shape, texture and statistical characteristics as well as relationships between linked levels of the segmented objects reported can enhance the accuracy of classification in the aerial photo [60,61]. For the accuracy assessment, the same sets of checking points for LU/LC map 2010 were employed here as the reference.

2.3. Urban Volume

To maintain the same spatial resolution for LU/LC maps and urban volume, we resampled the DSMs and LU/LC maps in the spatial resolution of 5 m.

2.3.1. DTM and Surface Feature Height (SFH) Determination

For SFH 2016, the DTM was provided by ground point’s extraction of the 2016 complete LiDAR data. The GCPs from the high accuracy measurement by geodetic GPS were used as the reference [33]. Since the DTM for SFH 2010 was not provided, the DTM generation from ALOS DSM should be conducted with the consideration of preserving the single vertical reference system. The slope based method that was initially proposed by [62] was applied for DTM determination. Their study has been reported as an important reference for selecting the DTM generation method according to specific terrain condition, especially in urban areas where complexity of terrain relief exists with a variety of non-ground features [24]. The slope-based method is the most effective approach for removing low vegetation, and the approach is most sensitive to the slope parameters in the thresholds [63]. Thus, we applied this method to minimize error in producing bare-earth where the low vegetation (e.g., sparse grass) was frequently found in the bareland class.

The DTM determination from ALOS DSM was conducted using three parameters, including the search radius in kernel pixel, the approximate slope of the terrain, and the confidence interval. Derivation of the slope threshold is implemented as the adaptive slope-based filter which works well in either for flat terrain or oblique terrain [64]. Therefore, in the study area with flat topography, this approach was appropriate. The slope threshold of 3% was set based on the LiDAR-derived DTM. The 95% confidence interval was settled as fixed parameter to reject the outliers of bare-earth points. Search radius with the kernel of 3, 5, 7 until 45 pixels associated with a grid size of 15 m, 25 m, 35 m, and so on, was applied to assess the best accuracy of DTM. The identified bare-earth points then were employed to produce DTM by using the B-spline interpolation approach.

After the DTM determination, the Surface Feature Height (SFH) was generated using the equation, as follows:

$${\mathrm{SFH}}_{\mathrm{i}}={\mathrm{DSM}}_{\mathrm{i}}-{\mathrm{DTM}}_{\mathrm{i}}$$

(1)

where SFH, DSM, and DTM are surface feature height (m), digital surface model (m), and digital terrain model (m), respectively, for pixel i. The ALOS DSM and DTM generated by ALOS were used to derive the SFH 2010. The LiDAR DSM and LiDAR-derived DTM were used to derive the SFH 2016.

2.3.2. Validation of SFH

The validation was conducted for SFH 2010 using SFH 2016 as the reference. The SFH was produced by the subtraction of DTM from DSM, as presented in Equation (1). The SFH 2016 used the DSM and DTM derived by LiDAR data with the specification that was described in Section 2.1. Herein, both DSM and DTM LiDAR had the high accuracy processed by rigorous methods. Therefore, the validation of SFH 2016 was not necessary. 30 buildings were selected as the height reference which was visible in all data, including the ALOS DSM, the ALOS image, the LiDAR DSM, and the orthophoto. Also, their heights and shapes had been confirmed without changes since 2010. The knowledge on the study area, the field survey to verify the building height and shape, as presented in the previous study [10], combined with the building map as the reference allowed in selecting those buildings carefully. Total 1946 pixels of 30 buildings were then used in validation with the height varied from 10.7 m to 41.8 m.

2.3.3. Built-Up and Green Volume

The built-up footprint and green footprint were extracted from built-up and green space classes, respectively. The equation to define the built-up volume (BV) is shown, as follows:

$${\mathrm{BV}}_{\mathrm{j}}=\mathrm{PA}\times {\mathrm{SFH}}_{\mathrm{j}}$$

(2)

where BV, PA, and SFH are the built-up volume (m³/pixel), pixel area (25 m²), and the surface feature height (m), respectively, of pixel j, a member of built-up class.

The equation for green volume (GV) is shown as follows:

$${\mathrm{GV}}_{\mathrm{k}}=\mathrm{PA}\times {\mathrm{SFH}}_{\mathrm{k}}$$

(3)

where GV, PA, and SFH are the green volume (m³/pixel), pixel area (25 m²), and the surface feature height (m), respectively, of pixel k, a member of green space class.

2.4. Built-Up and Green Expansions

Two different approaches were conducted to evaluate the built-up and green expansions. First, a patch analysis was adopted to assign the structure and configuration of horizontal (area) expansion. Secondly, the volume expansion rate was calculated to examine the vertical (volume) expansion in the consideration with the neighborhood growth. In this study, the term of area and volume expansions denoted to the increase of area and volume, respectively. Grid size with the best accuracy by the validation of SFH 2010 was applied to visualize the horizontal and vertical growth. The grid analysis is capable of figuring out the urban growth pattern and characterizing the spatial phenomena [65].

2.4.1. Patch Metrics of Built-Up and Green Area

In two metrics, we measured the configuration and spatial pattern of built-up and green patches in certain grids using the 8-neighbor rule (

Table 2. The selected patch metrics.

Name (Unit)	Abbreviation (Range)	Description	Equation
Area of built-up (hectares)	AREA (AREA > 0, without limit)	To examine the area expansion in term of size.	$\mathrm{AREA}={\mathrm{a}}_{\mathrm{ij}}\left(\begin{array}{c}\frac{1}{\mathrm{10,000}}\end{array}\right)$ aij = area (m2) of patch ij.	(4)
Contiguity index (none)	CONTIG (0 ≤ CONTIG ≤ 1)	To examine the pattern of area expansion regarding the spatial contiguity of patches.	$\mathrm{CONTIG}=\frac{\left[\frac{{\displaystyle \sum _{\mathrm{r}=1}^{\mathrm{z}}{\mathrm{c}}_{\mathrm{ijr}}}}{{\mathrm{a}}_{\mathrm{ij}}}\right]-1}{\mathrm{v}-1}$ cijr = contiguity value for pixel r in patch ij. v = sum of the values in a 3-by-3 cell template (13 in this case). aij = area of patch ij in terms of number of cells.	(5)

). These metrics were investigated based on area and shape features which are commonly applied to analyze landscape patterns [18,66]. The area features represent collection metrics assigning with the size of patches. The shape features assess the complexity of patch geometry, as well as the influence of the interaction of patch within each of the grid. The AREA metric is reported as the most important and useful component of the landscape information [66]. CONTIG metric assesses the spatial connectedness, or contiguity, of the pixel within a grid patch to provide an index of patch boundary pattern and thus patch shape. The larger contiguous patches produce a larger contiguity index [67]. The CONTIG metric shows the connectedness of patches based on the distance between neighboring patches of the same type. The CONTIG metric equals to 0 for single-pixel patch, and gains to a limit of 1 as the patch contiguity increases.

Then, we used Pearson correlation analysis to observe the relationships between the AREA and CONTIG metrics for both built-up and green area expansions.

2.4.2. Built-Up and Green Volume Expansion Rates

We assessed the volume expansion rate with the consideration of the mean expansion as representative of neighborhood expansion. Thus, the definition of volume expansion was addressed for built-up and green pixels having the volume increase more than a mean of volume expansion. For instance, in Figure 3, the mean of BV expansion was 591 m³/pixel. Thus, the BV expansion rate would be calculated for the pixel member of grid that had the BV expansion more than that mean value. The definition was based on the assumption that the volume increase for each grid would be different from one another. For instance, the sub-CBD area had more rapid growth than others, and the industrial area might have moderate expansion. Hence, it was unfair if the expansion rate was not considered on the typical growth for each grid.

Therefore, the expansion rate of BV (BR) was calculated using the equation, as follows:

$${\mathrm{BR}}_{\mathrm{i}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1\mathrm{to}\mathrm{n}}^{\mathrm{j}=1\mathrm{to}\mathrm{m}}{\mathrm{If}\mathrm{BV}}_{{\mathrm{j}}_{\mathrm{expansion}}}>{\mathrm{BV}}_{{\mathrm{i}}_{\mathrm{expansion}},\mathrm{mean}}}{{\mathrm{BV}}_{{\mathrm{i}}_{\mathrm{expansion}},\mathrm{mean}}}.\%$$

(6)

where BR, ${\mathrm{BV}}_{{}_{\mathrm{expansion}}}$ and ${\mathrm{BV}}_{{}_{\mathrm{expansion}},\mathrm{mean}}$ are the built-up volume expansion rate, the built-up volume expansion, and the mean of built-up volume expansion, respectively, of pixel j, a member of the built-up class for each grid i.

The expansion rate of GV (GR) was calculated using the equation, as follows:

$${\mathrm{GR}}_{\mathrm{i}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1\mathrm{to}\mathrm{n}}^{\mathrm{k}=1\mathrm{to}\mathrm{m}}{\mathrm{If}\mathrm{GV}}_{{\mathrm{k}}_{\mathrm{expansion}}}>{\mathrm{GV}}_{{\mathrm{i}}_{\mathrm{expansion}},\mathrm{mean}}}{{\mathrm{GV}}_{{\mathrm{i}}_{\mathrm{expansion}},\mathrm{mean}}}.\%$$

(7)

where GR, ${\mathrm{GV}}_{{}_{\mathrm{expansion}}}$ and ${\mathrm{GV}}_{{}_{\mathrm{expansion}},\mathrm{mean}}$ are the green volume expansion rate, the green volume expansion, and the mean of green volume expansion, respectively, of pixel k, a member of the green class for each grid i.

The ratio BR to GR (without unit) was also calculated to examine the comparison between the built-up and the green expansion rates. The balanced development was represented by the ratio BR to GR of value 1. The ratio BR to GR of less than 1 indicates that BR is lower than GR within a grid. While the ratio BR to GR of greater than 1 indicates that BR is higher than GR within a grid. Herein, we also denote “no BV expansion” (BR = 0), which means that a particular grid has no pixel with BV expansion more than the mean of BV expansion. Conversely, “no GV expansion” (GR = 0) indicates that a grid has no pixel with the GV expansion more than the mean of GV expansion. The linear regression analysis was applied to examine the expansion stability between built-up volume and green volume.

3. Results

3.1. Land Use/Land Cover (LU/LC)

For LU/LC 2010, the largest coverage is built-up class with an area of 9.49 km², and the smallest coverage is waterbody class with an area of 0.47 km². The built-up area is spread through the study area, especially in the north and east part (Figure 4a). The single residential houses are widely distributed evenly in the whole study area. While for the commercial buildings functioned as offices for banks and private companies are mostly found in the surrounding artery road of Ngagel and Darmo (see Figure 1b). For the surrounding artery road of Mayjend. Sungkono, many private companies and single shops are established. The institutional and company offices are frequently found in the sidelong the artery road of A. Yani. The big city park, namely KBS Park, has an area of 15.95 ha in a side of Darmo road. Other green spaces, which mostly include grasslands and scrublands, are distributed in the southern part. Cemeteries as a part of green space are also identified in the middle area. The backyard gardens of residential house with trees and grass are shown as the part of green space class. Bareland areas in large spots are mostly found in the southern part mixed with the grassland. The small spots of bareland are mostly spread in the western area. The historical Google Earth images coupling with the land use map 2008 are used to interpret the LU/LC map 2010.

The noticeable change is a replacement of the bareland area to the built-up area and green space. This alteration is in line with the increase of built-up area about 11.54% from 2010 to 2016 (

Table 3. Proportion of LU/LC for 2010 and 2016.

Type of LU/LC	Area (km2)		Change (%)
	2010	2016
Built-up	9.49	10.58	11.54
Road pavement	4.82	3.80	−21.30
Green space	3.09	6.06	95.61
Waterbody	0.47	0.69	46.93
Bareland	4.64	1.39	−70.00

). The significant increment also occurs for the green space portion as much as 95.61%. For LU/LC 2016, as presented in Figure 4b, the green space is evenly distributed throughout the study area. Some large spots of natural green space have shrunk. However, continuous patches of green space are visible alongside the roads, especially along the artery and primary roads. Green space is also seen in regular shape identified as new parks, which are a total of seven parks constructed. The small green patches are also recognized as the part of backyard residential house or commercial buildings, such as offices and shopping malls. The increase of the built-up area is mostly contributed by the decrease of road pavement portion by 21.30%.

The commercial buildings such as hotels and offices dispersedly exist sidelong of A. Yani road (see Figure 1b). While for the residential buildings, some mansions and apartments are found in the north eastern of Jagir Wonokromo road. The new buildings are mostly seen in the industrial area surrounding the Rungkut Industri road where some pavement areas have been transformed into fabrics and warehouses. Those areas were formerly used as a part of manufactory areas with asphalt or pavement surfaces which are functioned as a lot for placing huge container boxes. The increase of green space portion is obtained because of the removal of the bareland. However, in LU/LC 2016, bareland still appears in some large areas, which are frequently mixed with grassland. While for the waterbody, the increase is mostly contributed from several new ponds that are deliberately made for rainwater storage. The history function of Google Earth images combined with DSM, DTM, and orthophoto, have assisted in interpreting the LU/LC change.

The results of the accuracy assessment show OA about 83.2% and 86.43%, K about 0.79 and 0.83 for LU/LC map 2010 and 2016, respectively (

Table 4. Accuracy assessment of LU/LC.

Classes	LU/LC in 2010		LU/LC in 2016
	PA (%)	UA (%)	PA (%)	UA (%)
Built-up	85.84	89.82	87.61	95.52
Green space	81.9	80.00	82.86	84.00
Road pavement	83.00	83.80	87.00	86.10
Bareland	81.32	77.10	84.78	80.40
Waterbody	83.52	84.40	90.11	88.2
OA (%)	83.20		86.43
Κ	0.79		0.83

). The overall accuracy shows a good level within the standard range of 80–90% [19].

3.2. Urban Volume

3.2.1. DTM and SFH Determination

We used 11 numbers of kernel pixel associated with the grid size in meter unit to produce 11 DTMs (

Table 5. The lowest and highest values of DTM and SFH across grid size.

Kernel Size (pixel)	Grid Size (m)	DTM (m)		RMSe of Interpolation (m)	The Highest Values SFH (m)
		Lowest	Highest
3	15	−31.23	76.17	4.79	87.41
5	25	−30.65	61.77	4.87	89.43
7	35	−30.53	58.12	4.99	89.78
9	45	−30.47	54.27	5.04	89.01
11	55	−30.48	46.55	5.29	89.38
15	75	−31	39.23	5.6	91.03
21	105	−31.04	36.97	6.09	93.07
27	135	−31.03	30.24	6.42	93.79
33	165	−31.03	28.27	6.63	95.54
39	195	−31.05	25.68	7.29	95.85
45	225	−31.03	25.69	7.75	97.75

). For the lowest value of DTM in all mess sizes, no significant difference is recognized (Table 5). Nevertheless, for the highest value of DTM, the trend shows a substantial reduction from 76.17 m to 25.69 m. In terms of the LiDAR-derived DTM, the highest value presents at 26.7 m from Indonesia geoid. With the undulation height about 28.5–28.6 m for the study area [44], the highest DTM is approximately on 55.2–55.3 m from ellipsoid of WGS 1984. However, selecting the most accurate DTM is not straightforward, regarding the RMSe of interpolation and the highest value of SFH, which show the increasing tendency. Thus, validation of DTM in corresponding to the SFH maps is necessary.

3.2.2. Validation of SFH

30 buildings were selected as the reference height for validating the SFH 2010, as presented in Figure 5a. Total of 3 buildings have the height of 7–14 m with 3–4 stories, which are frequently functioned as single shops selling books, foods, or electronic items. 26 buildings have the height of 14–34 m with 5–10 stories which are mostly used for offices or shopping malls. While only one building with the height more than 34 m is used as a hotel in more than 10 stories.

Based on the trend of the average difference between SFH 2010 and SFH 2016, the inclination is presented in negative value from grid size 15 m to 35 m (Figure 5b). It reveals that the building heights of SFH 2010 are mostly lower than those of SFH 2016. The average difference reaches 0.3 m in grid size 45 m then, the graph line tends to increase again from grid size 55 m until grid size 255 m. The trend of RMSe shows the declination from grid size 15 m to grid size 45 m, afterward, the trend becomes the increase from grid size 55 m to grid size 255 m. The RMSe graph presents that RMSe of grid size 45 m has the lowest number of 5.03 m. This trend is consistent with the average difference, which shows that the grid size 45 m is the most accurate DTM.

3.2.3. Built-Up Volume (BV)

For 2010, the mean BV was 141.43 m³/pixel, the standard deviation was 103.12 m³/pixel, and the highest BV was 2227.42 m³/pixel (Figure 6a). Those values and SFH combined with the knowledge on the building characteristics and the previous study by [10], were allowed to classify the building height and number of storey. The mean BV was presented by the building with 1–2 stories which were frequently functioned as the single residential house [10]. The highest BV was shown by building functioned as a condominium in Mayjend. Sungkono road (see Figure 1b) with height about 89 m for 26 stories. Other buildings with more than 31 m were mostly established in the sub-CBD area. Those buildings were typically used for apartments or hotels. Buildings with the height about 14–23 m for 5–7 stories generally functioned as schools or offices were found in the southern area. Others were found in sidelong the Jagir Wonokromo road.

For 2016, the mean and the highest BV increased to 153.98 m³/pixel and 2993.48 m³/pixel, respectively (Figure 6b). The highest BV was recorded for the hotel located in A. Yani road with the height of 120 m and 36 stories. The buildings with more than 31 m were also identified in some spots sidelong other artery roads. The BV of 162.4–368.64 m³/pixel was frequently identified throughout the study area. These buildings generally had the height about 7–14 m with 3–4 stories, which were established for single shops. Interestingly, the BV of 575–781.13 m³/pixel was also identified in several locations along collector road, meaning that the development was not only in the sub-CBD area and the surrounded artery road, but also in the region with a high demand of residential, such as near the industrial area.

The most significant reduction was experienced by built-up with the height of 6–7 m about 16.41% and the increment was reached by built-up with the height of 7–14 m approximately on 15.96% (

Table 6. The proportion of built-up volume across the building height and number of storey for 2010 and 2016.

Built-Up Volume (m3/pixel)	Height (m)	Number of Storey	Portion
			2010	2016	Change (%)
<162.4	<7	1–2	82.07%	65.66%	−16.41
162.4–368.64	7–14	3–4	14.71%	30.67%	15.96
368.64–575	14–23	5–7	1.91%	2.13%	0.22
575–781.13	24–31	8–9	0.57%	0.51%	−0.06
>781.13	>31	>9	0.74%	1.03%	0.29

). Building reconstruction happened when the residents needed more spaces for their family member addressing the establishment of more stories in their houses. The new construction of the buildings with 3–4 stories was generally designed for single shops. The increase of the built-up with more than 9 stories was significant, showing a portion of 0.29%. The height associated with number stories was defined by the statistical analysis combined with an observation during the field work. The total BV also increased about 20.6% presenting the rapid development in the study area.

3.2.4. Green Volume (GV)

For 2010, the mean GV was 8.32 m³/pixel, the standard deviation was 69.46 m³/pixel, and the highest GV was 1125.35 m³/pixel (Figure 7a). Those values and SFH were used to classify the green space height. The mean GV was typically presented for grass with the height lower than 0.3 m, as a part of grassland and playground. The largest portion GV of 58.48% showed the green space with a height of 3–14 m, which was frequently found as tree planted in parks or cemeteries. The green space with the height more than 14 m was recorded as a tree which was naturally planted in several spots mixed with bareland.

The mean and the highest GV of 2016 gained at 16.79 m³/pixel and 1181.52 m³/pixel, respectively (Figure 7b). The largest portion GV of 49% remained on range 77.78–355.61 m³/pixel showing as trees distributed evenly in the study area, generally planted in the park or garden of residential single house or office. Several big patches of green space in this GV range were recognized as new parks. While another green space was recorded in a sufficient portion about 35.48%, generally as the trees planted along the road both artery or collector road.

The green space showed the largest increment of 8.99% with the height of less than 0.3 m.

Table 7. The proportion of green volume across the height for 2010 and 2016.

Green Volume (m3/pixel)	Height (m)	Portion
		2010	2016	Change (%)
<8.32	<0.3	3.73%	12.73%	8.99
8.32–77.78	0.3–3	30.68%	35.48%	4.80
77.78–355.61	3–14	58.48%	49.00%	−9.48
>355.61	>14	7.10%	2.79%	−4.31

shows that the portion of green space with the height of more than 3 m reduces significantly because of the removal of naturally planted trees in the grassland and scrubland for establishing built-up. Herein, the development of green space was addressed in another form, such as the street trees, which were planted in the side or the median of the road and sidelong the river.

3.3. Urban Expansion

Considering the most accurate of DTM derived from DSM ALOS and addressing the consistency, herein, the grid size of 45 m × 45 m (9 × 9 pixels) was applied in the patch metric analysis and the estimation of the volume expansion rate. There were 10,906 grids counting 133 rows and 82 columns.

3.3.1. Patch Metric of the Built-Up and Green Area Expansions

The built-up area (BA) expansion of 100.48 ha was mostly located sidelong the artery roads and the riverbank (Figure 8a). Applying the AREA metric, the mean value (μ) was 0.038 ha, meaning that the average area of built-up as much as 0.038 ha of 0.2025 ha (total area per grid) has expanded. The standard deviation (σ) showed 0.035 ha close to the mean value. The built-up area expansion was occurred in the vacant land, such as bareland, while a dense built-up area implied without expansion. Thus large variance was produced. The highest value was 0.2025 ha, indicating that all area within the grid was transformed to built-up (Figure 8b). Those values were applied to define five classes to present the level of BA expansion from the low to high with the AREA scheme: class 1 (value of 0), class 2 (0 − μ), class 3 (μ − μ + 1σ), class 4 (μ + 1σ − μ + 2σ), and class 5 (>μ + 2σ). The largest portion AREA of 41.97% was class 2 spread in the whole study area. The class 1 presenting no BA expansion had the portion of 17.39%, identified as road, bareland, or cemetery. For CONTIG metric of BA expansion, the mean, standard deviation, and the highest value were calculated at 0.38, 0.221, and 0.856, respectively (Figure 8c). Those values were then used to classify the CONTIG metric using the same scheme of AREA. The largest portion CONTIG of 41.34% was reached by class 3 indicating the scatter configuration, which was mostly presented by the residential area located far from the artery roads. In line with the AREA metric result, a high CONTIG metric (classes 4 and 5) was mainly found in a side of A. Yani, Mayjend. Sungkono and toll roads.

The green area (GA) expansion of 287.98 ha was mostly found in a side or in median road of both artery and collector roads (Figure 9a). The result of AREA metric showed the mean of 0.044 ha with the standard deviation of 0.032 ha (Figure 9b). Those statistical values were used to determine five classes using the same scheme, as presented in BA expansion metrics. Interestingly, the highest AREA of GA expansion was recorded on 0.2025 ha, indicating the whole area within the grid was fully expanded for the green class that was generally detected for the new park’s construction. Class 2 presenting disjoined patches had the largest portion of 45.2%. Class 1 had the portion of 9.88%, presenting no GA expansion, detected mainly as built-up and waterbody. The result of CONTIG metric showed the mean of 0.235 with the standard deviation of 0.17 and the highest value was 0.856 (Figure 9c). The largest coverage was shown by class 2 with a portion of 43.26% accumulated in several residential sites. Similarly to the AREA metric, the high class of CONTIG was also detected alongside A. Yani road.

The “Pearson correlation coefficient” showed a significant positive correlation between the AREA of BA expansion and CONTIG of BA expansion, also between the AREA of GA expansion and the CONTIG of GA expansion (

Table 8. Pearson correlation coefficients * between AREA metric of BA expansion, CONTIG metric of BA expansion, AREA metric of GA expansion, and CONTIG metric of GA expansion.

	CONTIG BA	CONTIG BA	AREA GA	CONTIG GA
AREA BA	1	0.762	−0.111	−0.12
CONTIG BA		1	−0.21	−0.05
AREA GA			1	0.735
CONTIG GA				1

Note: * 95% confidence interval.

). For BA and GA, both metrics, AREA and CONTIG, they had a high correlation, as being expected. The establishment of new built-up created the large AREA expansion and the high connectedness between built-up patches. While, the increasing GA by the construction of new parks or planting trees along the road and river affected the green configuration, which was also increasingly contiguous. However, both AREA and CONTIG of BA expansions had a negative correlation with the AREA and the CONTIG of GA expansions.

3.3.2. Built-Up and Green Volume Expansion Rate

The BV expansion was 10.92 million m³ with the mean of 176.71 m³/pixel presenting for the built-up establishment with 7 m in height. The BR generated the mean of 75.39%, the standard deviation of 6.27%, and the highest of 81.07% (Figure 10a). The low BR (<69.12%) was frequently identified in the residential area in the north eastern study area and surrounding sub-CBD. The high BR (>75.39%) had the largest portion of 43.27%, which was distributed sidelong the artery roads and the river.

The GV expansion was 7.65 million m³ with the mean of 94.55 m³/pixel denoting for the tree with the height around 3.78 m. The GR showed the mean of 78.38%, the standard deviation of 3.17%, and the highest value of 81% (Figure 10b). The low GR (<75.21%) was accumulated in KBS Park due to the reduction of the green space with high trees. The high GR (>79.96%) was generally identified in the residential sites.

We highlighted four sites, as presented in Figure 10, using circle with a diameter of 450 m (10 grids). The circle size was selected in the regard of volume expansion configuration, herein, the building size was also considered according to the study area characteristics: site 1 having high BR (>75.39%) and high GR (>78.38%), site 2 having high BR (>75.39%) and low GR (<75.21%), site 3 having low BR (<69.12%) and high GR (>78.38%), and site 4 having low BR (<69.12%) and low GR (<75.21%). Site 1 was identified in the sidelong river where the high BR was mostly generated by the establishment of new single shops or offices. The trees producing the high GR were planted in sidelong rives to prevent water runoff and reduce the threat of flooding. Site 2 was located in the side of A. Yani road where new commercial buildings functioning as offices and hotels were constructed. However, the low GR was identified due to the limited space for green space because of the dense buildings. Site 3 was detected in the dense residential area with the lack of space for new building construction. The trees in the backyard of the residential house had been growing up to produce the significant GR. While site 4 was located in the high density of commercial buildings where sufficient space for new construction of building and green space was not available.

The ratio BR to GR with a mean of 3.7 means that the development of BV and GV was imbalanced (Figure 11a). It indicates that the BV expansion was almost four times of GV expansion. “No BV expansion” was mostly accumulated along A. Yani and Jagir Wonokromo roads, also surrounding toll road. Conversely, “no GV expansion” was mainly located in several spots sidelong of A. Yani road. This class was also accumulated in several sites surrounding toll road. While for the ratio value of more than 1, it was mostly found sidelong of both artery and collector roads. The regression revealed a moderate R of −0.356 with R² of 0.127 (Figure 11b). This moderate level in negative correlation showed the increasing BR but the decreasing GR. It is demonstrated that the BR and GR ratio was not proportional as the scatter plots were mainly located on the upper right side.

4. Discussion

The BA and GA expansions were 100.48 ha and 287.98 ha, respectively. The low BA expansion was generally presented in the dense residential area in the north eastern part of Jagir Wonokromo road. The high BA expansion for offices or hotels was frequently identified in sidelong the artery roads. The GA expansion presented larger than BA expansion since the transformation of bareland to playground or park. The tree planting sidelong road was seen in sufficient numbers. The program by the government for increasing the green space accommodates the policy that new commercial building construction has to provide green patches, such as garden or backyard. However, the presence of green space portion in LU/LC 2016 revealed at 26.89%. Another strategy to accelerate the greening coverage should be considered especially in the dense residential area. The green roof or green garage may be adopted and applied in the local policy to create green neighborhood.

Applying the patch analysis, the result indicated that the mean of BA expansion was smaller than that of GA expansion. The connectedness of BA expansion was greater than the one of GA expansion. The BA expansion was smaller due to the limited land availability which caused the higher contiguity in the dense area. The low expansion for both BA and GA in discrete configuration was identified in the dense residential area. The sufficient BA expansion in the sprawl pattern was detected in industrial areas, which were dominated by the new fabric construction among the existed buildings. Interestingly, the sufficient GA expansion was also detected, which was mostly contributed by the construction of the new park in surrounding fabric to provide the green environment for their workers. Surprisingly, the high GA expansion was recorded in the sub-CBD area where the several buildings in a commercial block were established. In conjunction with the commercial block establishment, the garden designed as aesthetical façade was also constructed, then the high GA with an aggregate pattern was generated.

The result of the correlation analysis revealed that the BA expansion affected the connectedness, likewise for the GA expansion. Herein, new parks and green space into the new form of linear arrangement with high connectedness provided the high accessibility for residents. The residents with higher proximity to green space were satisfied with their accessibility to green space [68], and they paid significant amounts of money to get a house for living close to the park [69]. The spatial planning should be considered in the development of connective green space in neighborhood landscape at various levels [70]. This study found that the built-up establishment filled merely the bareland or vacant land, which was located among the existed built-up to facilitate the market-oriented facilities. Thus, the Surabaya government wished to create the policy of urban land transformation in the sustainable context.

The LU/LC change showed the significant BA and GA expansion of 11.54% and 95.61%, respectively. The rapid development had addressed the increase of built-up area for commercial use along artery roads. The economic indicator became the driving force for the urban development [53]. Thus, the industrialization and urbanization had stimulated economic growth demand on manufactory and housing [54].

The result confirmed to BV and GV expansion approximately on 20.6% and 54.9%, respectively. The larger portion BV expansion than BA expansion showed that built-up growth led to vertical development rather than the horizontal. The new development was mostly designed for buildings with medium height about 7–14 m functioned as single shops or offices. The vertical dimension of urban growth reported by the previous study can enhance the land use efficiency and the economic benefit, including employment opportunity and population density [55].

The prominent green area expansion was determined by the removal of the bareland. The new form of green had figured out the transformation from the big patches of natural grassland to a long structure of green shape, such as the street tree. The street trees play an essential role in the urban environment by improving air quality, providing shade, decreasing carbon footprint, and reducing the urban heat island effect [56]. Moreover, the long configuration of trees can link among green spaces to provide a cohesive node of Urban Green Infrastructure (UGI) [57].

The mean GR was higher than the mean BR for 75.394% and 78.375%, respectively. The built-up expansion was caused by the establishment of buildings with spontaneous development in the dense neighborhood. The green expansion was contributed by trees that were planted in the garden or as the street trees. The urban development integrated with small scale of green infrastructure provides an important mechanism for gaining the benefits in multi-functional level [70]. However, the ratio of BR to GR exposed that the imbalanced expansion between BV and GV was mainly detected as the new building construction in the surrounding artery and toll roads.

The further policy and strategy of sustainable planning should be a prime concern to lead the harmonious development of built-up and green environment. The big sites without GV expansion were mainly detected in the bareland, which could be reformed to efficient land use by integrating economic and environmental development. The commercial and residential building can be established in conjunction with the green space construction. Therefore, the workers or residents will be provided by the green neighborhood environment. Positively, the proposed land conversion can improve the imbalanced development (Figure 11b). The green volume with high trees might also be a concern for providing the high green quality environment. Regarding the landscape ecology, the dense vegetation can support a habitat diversity [71]. Moreover, the complex configuration and structure of vegetation are considered by residents to promote aesthetical aspect [72].

Our approach is in accordance with the framework by [73] in the regard of six fundamental aspects of urban form which are organized into three overarching components. We focused on the first component of material including built-up as aspect 1 and green space as aspect 2. The built-up feature was used to represent the built infrastructure of city, and the green signature was used to represent land surface features with biological activity. Both materials reflect the visible urban form, as well as they change over time [73]. The configuration as the second component was involved by employing the 2D and 3D space in the spatial pattern analysis. The building height shows the positive correlation with land use, and the height information reflects the plant growth over time. Without 3D information, high-rise building and single storey are represented in the same mode; moreover, the building volume can address the building occupancy or population density [74]. Timeframe as the third component was applied in bi-temporal (2010 and 2016) to investigate the built-up and green expansion. The time selection was considered for evaluating the green program, which was formally launched by The President of The Republic of Indonesia in Law of The Republic of Indonesia Number 26 Year 2007 Concerning Spatial Management [75].

Using the hybrid classification approach, the high overall accuracy can be achieved for both medium and fine resolution image. However, some limitations raised in this study regarding the accuracy of the classified ALOS image, the difference of the LU/LC maps, and the vertical reference system of DSM. The bareland and green space classes were speciously over-classified in ALOS image. More land area was assigned to those classes in the classified image compared to the ground truth references. Those over-estimations were indicated by their user accuracies being smaller than the producer’s accuracies (Table 4). The green space containing grass with drying appearance was classified to bareland, generally occurred in the residential area. On the other hand, the bareland was occasionally classified into green space. It normally occurred when the weed grass was mixed with the bareland. This situation was frequently happened in Asian countries, as reported by [57], where cropland is converted tentatively to bareland when the land is cleared for preparing the construction site, and this bareland would frequently be covered with weed grass during the transition period.

The LU/LC map differences were mostly caused by the difference in the spatial resolution of dataset used in this study. The classified ALOS image produced classes in more fragmented structure with slight boundary. Contrary, the classified orthophoto generated classes in more aggregated configuration with a clear boundary. The differences were mainly identified in the urban expansion area, especially in the residential areas and the commercial mixed residential areas. The overestimated bareland in LU/LC 2010 by ALOS image accounted for the reduction of the built-up and green area existence. The green space along the road was also failed to be detected in LU/LC 2010. When we checked through google earth image history, it appeared that small green patch was classified into bareland. Thus, the improvement of classification by applying an advanced fusion method will be proposed for a future study. The difference of the vertical reference system for DSM might also produce the gap between the generated SFHs. Even the highest DTM generated by grid size 45 m is close to the highest LiDAR-derived DTM. However, the systematic error may address because of undulation of the geoid. The future study is necessary to verify the detailed assessment.

5. Conclusions

This study empirically investigated the urban expansion of the horizontal and vertical development in the sub-CBD area of Surabaya. The geospatial approach was conducted for examining the urban form for bi-temporal (2010–2016) in conjunction with two material aspects (built-up and green types) and two configuration aspects (2D-3D space and spatial pattern).

The result revealed that the notable increment of built-up and green area was shown in the particular area. The green area expansion presented the adequate portion which was mainly contributed by the conversion of bareland to playground or parks. However, the built-up area expansions were less than the volume expansions, about 20.6%. It revealed that built-up growth led to the vertical development rather than the horizontal. The expansion of built-up area tended to scatter in fragment configuration. The expansion of green area tended to aggregate in the linear pattern. The high built-up volume expansion was frequently detected in sidelong artery roads where the high buildings for hotels or offices were established. While for green space, the high green volume expansion was partially detected in the low built-up volume expansion when the trees planted in sidelong roads or the backyard gardens had grown significantly. However, imbalanced development presented that the built-up growth was higher than the green growth, because more vertical buildings were established. Thus, the green improvement was needed to promote the sustainable neighborhood environment.

The proposed approach of volume expansion rate can assess the vertical development in consideration with the neighborhood growth. This approach was convenient to investigate the vertical development of built-up. However, the approach seems tough to apply for the green class, due to the irregular shape and the similar structure with road.

In this study, we used the multi-spatial resolution datasets due to the lack of the fine resolution data for 2010. A future study could refine the differences in the result generated by applying some advanced methods to enhance the accuracy of LU/LC. Also, the detailed assessment in DTM generation would improve the accuracy of SFH. Presenting the overall accuracy of LU/LC maps in the adequate threshold by applying the hybrid classification and the successful validation of SFH 2010, our result can be accepted for understanding the urban area extension and urban volume expansion. This study can be used as a reference to improve the neighborhood environment through greening programs in the context of sustainable urban development.