SUGARCANE PLANTATION MAPPING USING DYNAMIC TIME WARPING FROM MULTI-TEMPORAL SENTINEL-1A RADAR IMAGES

Updating of seasonal agricultural crop map is limited by the local knowledge of the mapper. Mapping of previously unaccounted agricultural plots involve massive field works aided by very high-resolution images. The phenological cycle of seasonal crops like sugarcane, with a range of ten (10) to twelve (12) months from planting to harvesting, exhibit a unique characteristic in terms of radar backscatter and time. In this paper, a pattern matching algorithm was tested to detect sugarcane plantations. Dynamic Time Warping (DTW), which was originally used for voice recognition, was used to detect sugarcane plantations from multitemporal Sentinel-1A images. Using known sugarcane plots, temporal signatures were gathered and used to detect other plantations in the area. The result helped the Sugar Regulatory Administration (SRA) in updating the inventory of sugarcane plantations faster with detection accuracy of more than 92 percent.


INTRODUCTION
Nationwide crop maps are integral in planning and decision making undertaken by mandated government agencies. These types of maps are usually inaccessible or outdated in the Philippines. One effort to do an updated agricultural map in a national scale was done by the Nationwide Detailed Resources Assessment Using LiDAR (Phil-LiDAR 2 Program) between 2014 and 2017 but the detection focused on areas with LiDAR data (Blanco et al., 2016) which is around 42.43% (Gatdula et al., 2017) of the land area of the country.
Large-scale crop mapping requires a large amount of data and manpower. When using remotely-sensed images, prior knowledge of crop location in an image helps. Mappers usually spend a considerable amount of time looking for crop areas to be digitized relying on visual classification keys and different image manipulation techniques to highlight the target crops. Although digitizing from Very High Resolution (VHR) image produces accurate plots, single scene image can only capture a snapshot of what is on the ground when the image was taken which poses a problem for areas with different planting or harvesting time. As an example, rice fields usually become bare soil or covered with grass for a few months in a season. Some farmers also practice multi-cropping adding complexity to the detection. This is why delineating seasonal crops need to be conducted in a multitemporal manner to properly determine what crops were grown in an area which may not be captured by a single image. Frequent cloud cover also adds to the problem which introduces gaps even in temporal datasets making temporal signature detection also difficult.

Statement of the Problem
The Sugar Regulatory Administration (SRA) of the Department of Agriculture (DA) embarked into a mapping project to delineate all sugarcane plantations in each sugarcane milling district to estimate potential yield to determine the country's remaining sugar reserves anytime throughout the year. The estimation of crop yield and growth monitoring is part of SRA's mandate to regulate sugar classification based on the estimated production of sugar in the country.
The current methodology employed by SRA involve the use of VHR optical images captured by Unmanned Aerial Vehicles (UAVs) and satellites in the known plantation sites. The plots were digitized in each milling district to produce a vector file which is one of the inputs in estimating probable yield.
Sugarcane plots appear like grasslands in some areas requiring cross-checking with other datasets which most of the time are unavailable. SRA needs reliable aide to digitizers to easily find and recognize sugarcane plantations with high certainty without needing to rely on outdated ancillary data and frequent site visits.

Objectives
The objective of this study is to develop a sugarcane mapping algorithm from multitemporal radar images. Radar satellite data were specifically chosen due to their cloud-penetrating characteristic and reliable spatial and temporal resolution. Temporal backscatter values from known sugarcane plantations will be used to find similar signature to create a map of sugarcane areas which can also potentially lead to discovery of previously unaccounted plantations.

Crop Mapping
Multi-temporal techniques in crop identification from satellite images has been implemented ever since the availability of temporal images from Landsat. Odenweller and Johnson (1984) were one of the first to develop a methodology in identifying crops from multi-temporal Landsat images based on their spectral characteristics. They used the Greenness Above Bare Soil technique to detect vegetation. Seasonal crops exhibit a pattern wherein low detectable vegetation cover can be observed at the start of the planting season that slowly increases in value and decreases again after harvesting. They even managed to identify different crops by matching a crop calendar to the observed temporal variation to determine which crops were planted.
A similar phenological sequence technique was implemented by Bargiel (2017), this time using multi-temporal radar images from Sentinel 1. He compared the detection of different crops based on Phenological Sequence Patterns (PSP), Random Forest and Maximum Likelihood Classifier. He demonstrated that the inclusion of PSP improves the detection of individual crops.

Dynamic Time Warping
In the absence of phenological sequence data and crop calendar, pure sequential backscatter is the only data that can be extracted purely from a temporal radar image stack with prior knowledge on the length of the cropping season. In this study, pure signature matching technique is employed.
Dynamic programming techniques in aligning sequences like this has been widely applied in biotechnological applications (Eddy, 2004). With a known temporal signature as template, it can be compared to a second series and determine if the same pattern appeared. Exact matching is nearly impossible to achieve. Sometimes, the length and magnitude of a sequence are different from the template. This type of problem is usually encountered on voice recognition problems which varies in timing, pitch, or pronunciation. The solution to this dynamic programming problem was proposed by Sakoe and Chiba (1978) by introducing symmetric and asymmetric pattern matching and with slope constraints. This evolved into the Dynamic Time Warping (DTW) technique which can align a time series to a specific template by minimizing a distance measure. The technique can be applied in detecting sugarcane plantations due to variabilities in planting and harvesting time producing shorter or longer temporal signature.
This study capitalizes on the warping capability of DTW to capture the early or late harvesting time of sugarcane whose season ranges between 10 to 12 months.
The DTW technique has been applied before in crop mapping. Guan et al. (2016) used temporal Normalized Difference Vegetation Index (NDVI) values from MODIS data in the Mekong River Delta to detect rice but found out that the resolution of the dataset is too coarse to properly identify the said crops. Belgiu and Csillik (2018) also implemented a variation of DTW which penalizes longer time series (Time-Weighted DTW) using Sentinel 2 images with better resolution at 10 meters and achieved promising results.

R Software
To implement the DTW algorithm, the R statistical software was chosen. R is an open source software that provides a comprehensive implementation of different DTW algorithms. It supports arbitrary local (e.g. symmetric, asymmetric, slopelimited) and global (windowing) constraints, fast native code, and several plot styles, etc. The DTW implementation used in this paper was based on the method of Sakoe and Chiba being the basic implementation of the technique. The title should appear centred in bold capital letters, at the top of the first page of the paper with a size of twelve (12) points and single-spacing. After one blank line, type the author(s) name(s), affiliation and mailing address (including e-mail) in upper and lower case letters, centred under the title. In the case of multi-authorship, group them by firm or organization as shown in the title of these Guidelines.

Study Area
The study area is a portion of the province of Tarlac as shown in Figure 1 which covers parts of Tarlac City, Municipalities of La Paz, Capas, and Concepcion; all under the Province of Tarlac. The area, which is part of the Tarlac Mill District found in Central Luzon, is one of the major producers of sugar in the region. The total sugarcane area in the district is around 10,695 hectares with an average cane yield of 51.42 tons of cane per hectare (Sra.gov.ph, 2016).   ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-3-2020, 2020 XXIV ISPRS Congress (2020 edition)

Sugarcane Plantation Vector Data:
A vector file of the sugarcane plots in the area was obtained from the SRA. The data was digitized from VHR images and validated in the field. The dataset was divided into three areas; Site A (green) for the extraction of sugarcane temporal signature, Site B (blue) for the determination of parameters, and Site C (red) which is the biggest, served as validation site as shown in Figure 2. (14) PlanetScope scenes from December 20, 2016 to March 16, 2018 were provided by the Philippine Earth Data Resource and Observation Center (PEDRO). The images were used to verify the vector file provided by SRA analyzing the temporal NDVI of all the plantations within the study area.

Process Workflow
The processing steps can be summarized into four (4) major steps; pre-processing of Sentinel 1A data, temporal signature extraction, DTW processing, and sugarcane classification as shown in Figure 3.

Pre-processing:
All Sentinel 1A images were preprocessed by applying orbit correction, radiometric calibration to beta nought, multilook, radiometric terrain flattening, coregistration, and range doppler terrain correction using ESA's Sentinel 1 Toolbox. Shuttle Radar Topography Mission's (SRTM) 1-second DEM product was used for geocoding resulting to a 42-band with 15-m spatial resolution data. Separate files for the co-polarized (VV) and cross-polarized (VH) backscatter values were generated.

Temporal Signature Extraction:
Temporal signature for both VV and VH polarizations were extracted from a plot of sugarcane plantation with verified similar planting and harvest date. The average backscatter per band and polarization was related to the phenological cycle of the crop. Data points before and after planting and harvesting were removed to generate a complete temporal signature corresponding to one phenological cycle. Figure 4 shows the average backscatter values over the phenology of sugarcane.

DTW Processing:
Processing is performed through a Python script using GDAL, rpy2, and multiprocessing package to read/write raster files, call R commands, and do parallel processing, respectively. Both VV and VH image stacks were split into blocks and assigned to different compute nodes for processing. The DTW distances were calculated between the template signatures and the temporal series found in each pixel of the image stacks within the study area specifically for Sites B and C. The parameters used for the DTW distance computation are listed in Table 2. Additional data (Pearson correlation for both VV and VH) were computed to aid in the classification. The final raster file produced contains the information found in Table  3.

Name Description Value Distance Method
Pointwise (local) distance function to use.

Euclidean
Step Pattern Lists the transitions allowed while searching for the minimum distance path. Size of the windowing function.

Open End
When enabled, the end-point constraint is relaxed.

Open Begin
When enabled, the start-point constraint is relaxed.
False Table 2. DTW parameters and values used.

Classification/Determination of Threshold Values: DTW Distance, Pearson Correlation VV, and Pearson
Correlation VH of confirmed sugarcane plots from Site B were extracted from the resulting raster data. The Pearson correlation between the template temporal signature and the matched pattern was utilized as an additional measure of similarity between the two signals. This is due to some instances where the DTW distance is very low, but the pattern is not exactly followed. In principle, the DTW distance must be low and the correlation values must be close to 1. To remove outliers from the data, only the ninety fifth (95 th ) percentile of the DTW distance values and the fifth (5 th ) percentile of the correlation values were retained. The histograms of the DTW distances and correlation values for VV and VH are shown in Figure 5, 6, and 7, respectively. It was found that the combination of the values of less than 1.71643 for DTW distance, greater than 0.50606 for VV correlation, and greater than 0.57862 for VH correlation comprise the 95 th percentile of the DTW distances and the 5 th percentile of the correlation values found in Site B. The final threshold values are tabulated in Table 4.

RESULTS AND DISCUSSION
Upon extraction of temporal signature from Site A, the DTW algorithm was applied in Site B using the signature as template.
The DTW distance computed together with correlation values of VV and VH polarizations were also computed. The ordinary DTW algorithm can only align a pair of one-dimensional time series at once. Since we have two stacked images--(VV and VH image), multivariate DTW was applied. Applying traditional DTW not only makes each pairwise alignment independent but it can also cause each dimension to be compared separately. It can also utilize all dimensions to obtain the optimal path and align multiple signals simultaneously. After alignment, the regular similarity measurement can be used to the aligned signals. It was observed that some data (which belongs to the discarded 5 th percentile) lie on the edges of sugarcane plots (see Figure 9) which may contain pathways or irrigation facilities which are also sometimes covered with grass. These points are considered as noise and were discarded to preserve the integrity of the gathered temporal signature.
A 3x3 smoothing filter was applied to Site C before applying the thresholds from Table 4 to remove single-pixel noises. All pixels were classified into Sugar and Non-Sugar classes where random points were extracted to calculate the accuracy of the method.
The technique attained an overall accuracy of 92.75% with User's accuracy of 88.80%, Producer's Accuracy of 89.47% (Table 5), and Kappa Index of 0.84 which indicates high confidence in its detection rate. Upon visual inspection of misclassified points for sugar, it was observed that some of them were cogon grass.

CONCLUSION
The sugarcane plantation detection was implemented using basic DTW algorithm which locates similar temporal signatures from a known plantation within a temporal image stack. The detection technique uses Radar images from Sentinel 1A satellites to take advantage of its temporal reliability and cloud-penetrating characteristic avoiding data gaps which is very important in signature-matching.
The overall accuracy of 92.75% attained by the proposed methodology is very high in the standards of image classification in remote sensing, making it very promising in detecting seasonal crops which exhibit a similar pattern, wherein crops have low backscatter values at the start of the planting season while gradually increasing at different growth stages and drops again upon harvest. The technique is also being implemented to detect rice and corn plantation in cooperation with mandated agencies in the Philippines.