A Multi-Decadal Change Analysis for Irrigation Ponds in Taoyuan, Taiwan Using Multi-source Data

A multi-decadal change analysis of the irrigation ponds in Taoyuan, Taiwan was conducted by using multi-source data including digitized ancient maps, declassified single-band CORONA satellite images, and multispectral SPOT images. Supervised LULC classifications were conducted using four textural features derived from the single-band CORONA images and spectral features derived from SPOT images. Post-classification analysis revealed that the number of irrigation ponds in the study area decreased during the post-World War II farmland consolidation period (1945 – 1965) and the subsequent industrialization period (1970 – 2000). However, efforts on restoration of irrigation ponds in recent years have resulted in gradual increases in the number (9%) and total area (12%) of irrigation ponds in the study area.


Introduction
Urbanization and anthropogenic disturbance can greatly alter the landscape pattern in a region. Understanding the process of landscape changes can lead to better urban and regional planning and management. For example, Chen et al. [1] investigated decadal land-cover changes in Taipei and found that urbanization resulted in a greater increase in urban heat absorption than in thermal inertia which in turn caused increase in urban diurnal land surface temperature. They suggested using measures that can reduce urban heat absorption (such as urban parks, community green spaces, green roofs and permeable pavements) to mitigate the adverse effects of urbanization. Sun et al. [2] studied the long-term trends and patterns of urban growth in Guangzhou, China and identified three urban growth types of infilling growth, outlying growth and edge-expansion growth.
Studies have shown that, comparing to urban area, paddy fields and their surrounding areas in East Asia region have lower ambient air temperatures [3][4][5][6]. Cheng et al. [6] studied land-cover specific surface and ambient air temperatures in the Taoyuan County of northern Taiwan. Comparing to urban built-up land-cover, ambient air temperatures of the irrigation ponds and paddy fields were roughly 3 -4 and 1 -2 degrees (C) lower, respectively. These studies demonstrate that irrigation ponds and paddy fields play an essential role in regional and local climate. Studies have also revealed the multifunctionality of irrigation ponds and paddy fields in East and Southeast Asian countries, including flood mitigation, temperature cooling, groundwater recharging, providing habitat for inhabitants, etc. [7 -12].
Many irrigation ponds scatter over the Taoyuan City in norther Taiwan nowadays. Construction of these irrigation ponds can be dated back to the early 18th Century when immigrants from the southeast region of Mainland China moved to northern Taiwan [13]. Most of these irrigation ponds were constructed in lower plain areas. At its peak, there were more than 3,200 irrigation ponds (only ponds with an area larger than 0.3 hectare were counted) in the Taoyuan area. These irrigation ponds have become a crucial element for paddy rice farming in the region for over two centuries. After the World War II, a major reservoir, the Shihmen Reservoir, was constructed in 1964 to supply irrigation water for downstream paddy fields and to boost the agricultural productivity. With the construction of the Shihmen Reservoir and subsequent industrialization and urbanization in the region, the number of irrigation ponds has greatly decreased. Industrial parks, residential blocks and transportation facilities replaced irrigation ponds and paddy fields, and Taoyuan City is now one of the most populated cities in Taiwan. In recent years, there have been growing concerns of the deteriorating environmental and ecological conditions in the Taoyuan area due to land-use/land-cover (LULC) changes. Quantitative assessment of the environmental and ecological effects of LULC changes requires an historical account of various LULC types and their spatial distributions. Such information and documents are largely unavailable or less accurate for the period prior to 1970s. Therefore, the objective of this study is to conduct a multi-decadal change analysis of the irrigation ponds in Taoyuan area by using multi-source map and remote sensing data.
The remaining parts of this paper is organized as follows. In Section 2, we describe the study area and multi-source data used in this study. In Section 3, we introduce the methods adopted for this study, including four textural features, details for calculation of textural features, and the maximum likelihood classifier. In Section 4, we present the textural features derived from single-band images and LULC classification results of different years. Detailed discussions on multi-decadal changes of irrigation ponds during three chronological periods are also given in Section 4. A summary of the findings and concluding remarks are presented in Section 5.

Study Area and Data
Twelve districts in the Taoyuan City and two townships in the neighboring Hsin-Chu County were selected as the study area ( Figure 1) which encompasses approximately 900 km 2 . In order to investigate multi-decadal changes of irrigation ponds in the study area, maps and remote sensing data of different periods (1904,1969,2004, and 2014) were collected.
Digitized ancient land maps (originally, in scale of 1:20,000) of the study area are available from Taiwan Ancient Land Map Digital Archive (TALM) [14]. The ancient land maps, originally created in 1904, reveal land-cover types and their spatial distributions within the study area in the beginning of the 20th Century. A total of 22 TALM maps covering the study area was collected. An example of the TALM maps is shown in Figure 2.
Four frames of the declassified military satellite CORONA KH-4B imagery of the study area acquired in February, 1969 were collected [15]. The CORONA satellite acquired photographs with a telescopic camera system. High performance photogrammetric film scanners were used to create digital black-and-white images at 7-micron (3,600 dpi) or 14-micron (1,800 dpi) resolution. The CORONA images collected for this study have a spatial resolution of 6 feet (approximately 1.8 meters). Figure 3 shows the CORONA black-and-white images of the study area in 1969.
SPOT-5 HRG multispectral images of the study area acquired in January 2004 and November 2014 were also collected. Table 1 shows the characteristics of SPOT-5 HRG multispectral images. These images were preprocessed for radiometric and geometric corrections by the Center for Space and Remote Sensing of the National Central University in Taiwan. We purposely chose images acquired in January and November, during which paddy fields were left fallow (without standing water), to avoid misclassification between paddy fields and irrigation ponds. A pseudo-color SPOT-5 image of the study area is shown in Figure 4.

Methods
Multi-decadal change analysis of irrigation ponds in the study area consists of two steps: (1) identification of irrigation ponds in the study area in 1904,1969,2004, and 2014 using TALM maps and CORONA and SPOT-5 remote sensing images; and (2) post-classification change analysis.
Irrigation ponds in the study area at the beginning of the 20th Century were extracted from TALM maps which had been digitized through a Map and Remote Sensing Images Digital Archive Project (MIDAP) [16]. Multivariate supervised classification was applied to the CORONA and SPOT-5 multispectral images for land-cover classification.
The original CORONA black-and-white images have only one broadband 8-bit (0-255) grey levels and the multivariate supervised classification cannot be directly applied. Therefore, instead of using the spectral information of CORONA images, several textural features were calculated and used as the classification features in subsequent multivariate supervised LULC classification. For SPOT-5 images, reflectances of the green, red and near infrared bands were used as classification features. In this study, we adopted the maximum likelihood classification method for LULC classification. Details of the maximum likelihood classification method have been widely documented [17,18].

Textural Features
Textural features are measures that characterize the structural or spatial variation of the variable under consideration within a pre-specified window. Four textural features -homogeneity (HOM), entropy (ENT), contrast (CON), and angular second moment (ASM), described below were considered in this study.
Given an image with M rows and N columns, we define ) , ( n m I as the grey level of the pixel located at the m-th row and n-th column of the image. The grey level difference of the pixel located at ) , ( n m , i.e. the central pixel, and its neighboring pixels can be represented by where 1 d and 2 d represent the distance (in unit of pixels) between the central pixel and its neighboring pixels in the y (row) and x (column) directions, respectively. Generally, the values of 1 d and 2 d are set to be either 0 or predetermined d  . Therefore, a pixel-pair can only align in one of the four angular directions - 0 ,  45 ,  90 , and  135 , as depicted in Figure 5 for the case of 3  d . Without considering the exact location of the central pixel, the grey level difference g is dependent on both the angular direction  and the distance d , and can be expressed as Hereinafter, textural features calculated using grey level difference of distance d is termed d distance  textural features. The d distance  textural features represent the average spatial relationship of grey levels in a moving block or kernel window of size × (in unit of pixels) . Values of textural features computed from a moving block are assigned to the central pixel. The moving block is moved, one pixel per move, around the whole image, except the peripheral pixels, to generate a new set of images of textural features. The size of moving block can influence the results of land-cover classification using textural features [19 -23]. The optimal size of the moving block is dependent on the pixel resolution and the typical size of land-cover objects. Small block sizes cannot capture the particular pattern of most land-cover classes, while the large ones can include pixels from more than one class [20]. In this study the following textural features were calculated, using distance-1 grey level differences of 8-bit CORONA images, and used for land-cover classification.
and N represent the count (frequency) of grey level difference of value g and the total number of pixel pairs defined by Equation (1) [19] and is adopted in this study. We tested various moving block sizes for calculation of textural features and found that textural features derived from moving blocks of size

13
pixels yielded better results of irrigation pond identification. Details of the test are described in Section 4.

Multivariate Supervised Classification Using Spectral and Textural Features
The maximum likelihood classification method was adopted for LULC classification in this study. Although not necessarily required, many LULC applications assume multivariate Gaussian distributions for the classification features of different LULC types. For LULC classification of CORONA images, textural features derived from moving blocks of size

13
pixels and distance-1 grey level differences were used. Land-cover types for LULC classification of CORONA images include ponds, paddy fields, buildings, and line features. Line features may include roads and man-made irrigation channels owing to their similar textural patterns. Unlike LULC classification of CORONA images, spectral reflectances of the green, red, and near infrared bands were used as classification features for LULC classification of SPOT-5 images. Four land-cover types including ponds, paddy fields, woods, and built-up were considered in LULC classification.
For Let = ( , … , ) be a k-dimensional feature vector of a particular pixel. The joint Gaussian density of the ith class A pixel with feature vector X is assigned to the ith class if ) ( X d i is the highest of all class-dependent discriminant functions; that is, Figure 6 shows the CORONA image (after mosaic and georeferencing) of the study area. Subarea A (4340 x 4340 pixels) was used for investigating the spectral and textural features of different land-cover types and their differences. Subarea B (990 x 998 pixels) was used for determining the optimal size of moving block (kernel window) for textural features calculation. From Subarea A, parcels of four land-cover types, namely ponds, paddy fields, buildings, and line features (for examples, roads and man-made channels) were selected for evaluating their spectral and textural differences. Figure 7 demonstrates that using only the spectral feature (grey levels) these land-cover types cannot be differentiated. By contrast, Figure 8 shows that these land-cover types can be better separated by textural features of ASM, CON, ENT, and HOM. Generally speaking, ponds are well separated from buildings and line features for all textural features, whereas ponds and paddy fields are less distinguishable.

Textural Features of the CORONA Images
(a) (b) Figure 7.  moving block sizes. The ground truthing/training data and producer's accuracies of different land-cover types with respect to various moving block sizes are shown in Figure 9. Figure 10 illustrates  4 m), it achieved the highest producer's accuracies for ponds and paddy fields. Since ponds are less likely to be misclassified to buildings and line features and vice versa, the optimal moving block size was then set to be 13 × 13 pixels under which the producer's accuracies for ponds and paddy fields were highest. Using this optimal moving block size, four CORONA textural feature images of the whole study area were generated and used for multivariate supervised LULC classification. Figure 11 illustrates CORONA textural feature images of the Subarea A.   Figure 11. CORONA textural feature images of the Subarea A using a 13 × 13 pixels moving block. (a) ASM; (b) CON; (c) ENT; (d) HOM.

LULC Classification Results
Classification results of the CORONA and SPOT images were evaluated by using the confusion matrices presented in Table 2. For LULC classification of CORONA images using textural features, most of misclassified ponds (omission error) were assigned to paddy fields and most of mis-commissioned ponds (commission error) were from paddy fields. Producer's and user's accuracies of ponds are 96% and 56%, respectively. These results indicate that although textural features enable conducting multivariate LULC classification using single-channel (black-and-white) CORONA images, noticeable commission error of ponds from paddy fields exists. For LULC classification of SPOT images using spectral features, misclassification between ponds and other land-cover types are rare, with the producer's and user's accuracies of ponds close to or higher than 98%.
It is worthy to note that, by using remote sensing images acquired in periods during which paddy fields were left fallow (without standing water), ponds and paddy fields are texturally similar but spectrally different, whereas the paddy fields and built-up are texturally different but spectrally similar. As a result, for textural-features-based classification using CORONA images most classification confusions occur between ponds and paddy fields, whereas for spectral-features-based classification using SPOT images most confusions occur between paddy fields and built-up.

Post-classification Processing of CORONA Images LULC Classification
The fact that LULC classification using textural features of 1969 CORONA images yielded lower user's accuracy (56.2%) for ponds indicates that the number of ponds were significantly overestimated. To circumvent such results, locations of the identified irrigation ponds in the 1969 CORONA mages were examined through the following post-classification processing to exclude unreliable ponds: (1) Identify individual ponds using a GIS based raster-to-vector conversion, (2) Ponds which only appeared in the 1969 CORONA images, but not in the TALM maps and SPOT images of 2004 and 2014 were excluded. (3) Woods are mostly located in mountainous regions where irrigation ponds do not exist. Also, LULC classifications using SPOT images achieved very high producer's and user's accuracies for the woods land-cover type. Therefore, all irrigation ponds located in mountainous regions were excluded. (4) Ponds with areas smaller than 0.3 hectares were excluded. This threshold is based on the definition of irrigation ponds by the Taoyuan Irrigation Association.

Multi-decadal Changes of Irrigation Ponds
The numbers and total areas of irrigation ponds estimated using TALM, CORONA and SPOT images in different years are shown in Figure 12. Generally speaking, the number of irrigation ponds have decreased significantly from the beginning of the 20 th century. The earliest irrigation ponds were established in the early 18 th century by immigrants from southern China. By the middle of the 19 th century, population in the study area grew fast and widespread irrigation ponds and irrigation canals were established. For example, an irrigation pond (the Hsiao-Li pond) and its connecting irrigation canal (the Hsiao-Li canal) were built in 1741 and are still in use. Three chronological periods, namely the post-World War II period (1945 -1965), the industrialization period (1970 -2000), and the restoration period (after 2000), and several policy decisions which contributed to decadal variation of irrigation ponds are also shown in Figure 12. Figure 13 shows the decadal variation and spatial distribution of irrigation ponds. During the post-World War II period (1945 -1965), Taiwan instituted the Arable Rent Reduction Act (ARRA) and land reforms which consolidated farmlands, improved irrigation systems and helped to boost the economic growth. Farmland consolidation resulted in larger paddy plots which were more suitable for mechanized farming and many ponds were converted to paddy fields or irrigation canals. Construction of a multi-purposes (irrigation water supply, power generation, flood prevention, domestic and industrial water supply, and tourism, in order of priorities) reservoir, the Shihmen Reservoir, in upstream of the study area started in 1956 and completed in 1964. After completion of the Shihmen Reservoir, the number of irrigation ponds largely decreased since the need for irrigation water supply by irrigation ponds sharply reduced. In around 1970's, many manufacturing factories and industrial parks were established and the study area gradually emerged to become a center of manufacturing industry in Taiwan. State highways and other transportation networks were also constructed after mid 1970's, including the Taiwan Taoyuan International Airport, Taiwan's largest international airport. Industrialization resulted in disappearance of many irrigation ponds (approximately 75% of the number of ponds in 1969), particularly in the Luju, Taoyuan, Badeh, Chungli, and Dayuan Districts (in blue color in Figure 13) where most manufacturing factories and industrial parks are located. After more than three decades of industrialization, local residents and government agencies became aware of the landscape aesthetics and ecological and flood prevention functions of the irrigation ponds. Two laws, the Taoyuan Ponds and Canals Restoration Program (PCRP) and the Taoyuan Pond and Canal Conservation and Restoration Incentive Act (PCRIA), aiming for conservation and restoration of irrigation ponds came into effect in 2003 and 2009, respectively. Such restoration efforts have resulted in gradual increases in the number (9%) and total area (12%) of irrigation ponds in recent years.

Summary and Conclusions
A multi-decadal change analysis of the irrigation ponds in Taoyuan City, Taiwan were conducted by using multi-source data, including digitized ancient land maps, declassified military satellite CORONA KH-4B single band imagery, and SPOT-5 HRG multispectral images. Supervised LULC classifications using textural features derived from the CORONA images and spectral features of the SPOT images were applied for identification of irrigation ponds in different years and for multi-decadal change analysis. Several concluding remarks are drawn as follows: 1. Textural features may be able to differentiate land-cover types which have similar spectral characteristics. In our study, ponds are well separated from buildings using textural features but they are nearly indistinguishable using single-band grey levels.
2. The optimal moving block size for textural features is dependent on the pixel resolution and the typical size of land-cover objects. For usage of the CORONA images in Taoyuan, the optimal moving block size is 13 × 13 pixels (or equivalently, 23.4 m  23.4 m).
3. Ponds and paddy fields are texturally similar but spectrally different, whereas the paddy fields and built-up are texturally different but spectrally similar. For textural-features-based classification using CORONA images most classification confusions occur between ponds and paddy fields, whereas for spectral-features-based classification using SPOT images most confusions occur between paddy fields and built-up. However, efforts on restoration of irrigation ponds in recent years have resulted in increases in the number and total area of irrigation ponds.