Satellite Constellation Reveals Crop Growth Patterns and Improves Mapping Accuracy of Cropping Practices for Subtropical Small-Scale Fields in Japan

: Mapping of agricultural crop types and practices is important for setting up agricultural production plans and environmental conservation measures. Sugarcane is a major tropical and subtropical crop; in general, it is grown in small ﬁelds with large spatio-temporal variations due to various crop management practices, and satellite observations of sugarcane cultivation areas are often obscured by clouds. Surface information with high spatio-temporal resolution obtained through the use of emerging satellite constellation technology can be used to track crop growth patterns with high resolution. In this study, we used Planet Dove imagery to reveal crop growth patterns and to map crop types and practices on subtropical Kumejima Island, Japan (lat. 26 ◦ 21 (cid:48) 01.1 (cid:48)(cid:48) N, long. 126 ◦ 46 (cid:48) 16.0 (cid:48)(cid:48) E). We eliminated misregistration between the red-green-blue (RGB) and near-infrared band imagery, and generated a time series of seven vegetation indices to track crop growth patterns. Using the Random Forest algorithm, we classiﬁed eight crop types and practices in the sugarcane. All the vegetation indices tested showed high classiﬁcation accuracy, and the normalized di ﬀ erence vegetation index (NDVI) had an overall accuracy of 0.93 and Kappa of 0.92 range of accuracy for di ﬀ erent crop types and practices in the study area. The results for the user’s and producer’s accuracy of each class were good. Analysis of the importance of variables indicated that ﬁve image sets are most important for achieving high classiﬁcation accuracy: Two image sets of the spring and summer sugarcane plantings in each year of a two-year observation period, and one just before harvesting in the second year. We conclude that high-temporal-resolution time series images obtained by a satellite constellation are very e ﬀ ective in small-scale agricultural mapping with large spatio-temporal variations.


Introduction
Agricultural land is a key component of land use and land cover. Nearly 40% of the Earth's land surface is now being used for agriculture [1]. The expansion of agricultural land to increase food production and intensification of its use [2] has altered the structure and function of ecosystems and in some cases has diminished their ability to continue providing valuable resources [1,3,4]. Many land management applications require reliable and timely land cover mapping over large heterogeneous landscapes [5]. Therefore, mapping and monitoring of agricultural land use is an important topic of high public interest.
The use of satellite remote sensing to monitor the use and management of agricultural land over a large area is an efficient approach (e.g., [6][7][8][9][10][11][12][13][14][15][16][17]). In agricultural mapping, vegetation indices are used to The study area was an agricultural area of 17.1 km 2 (Ministry of Agriculture, Forestry and Fisheries, Japan, MAFF; https://www.maff.go.jp/j/tokei/kouhyou/sakumotu/menseki/#c), which is covered by sugarcane, pasture, purple yam, pineapple, agricultural facilities such as greenhouses, and other crops such as squash, okra, and rice. Sugarcane is the predominant crop [36]. On Kumejima Island, there are three key sugarcane management practices (Table 1): (1) Ratooning, in which sugarcane grows from sprouts of underground stubble left in the field after harvest of the main crop (i.e., ratoon), and is harvested about every 12 months, usually for at least 4 years, and then the crop is renewed due to decreasing yield; (2) spring planting of seedlings (between February and May, just after harvesting); and (3) summer planting of seedlings (between July and October). The sugar industry equipment usually operates from January to March. In ratooning, sugarcane is harvested 12 months after harvesting in previous year, in spring planting it is harvested 12 months after planting, and in summer planting it is harvested in the sugar-production period 18 months after planting. About 63% of farmers manage a sugarcane area of less than 1 ha [48]. The study area was an agricultural area of 17.1 km 2 (Ministry of Agriculture, Forestry and Fisheries, Japan, MAFF; https://www.maff.go.jp/j/tokei/kouhyou/sakumotu/menseki/#c), which is covered by sugarcane, pasture, purple yam, pineapple, agricultural facilities such as greenhouses, and other crops such as squash, okra, and rice. Sugarcane is the predominant crop [36]. On Kumejima Island, there are three key sugarcane management practices (Table 1): (1) Ratooning, in which sugarcane grows from sprouts of underground stubble left in the field after harvest of the main crop (i.e., ratoon), and is harvested about every 12 months, usually for at least 4 years, and then the crop is renewed due to decreasing yield; (2) spring planting of seedlings (between February and May, just after harvesting); and (3) summer planting of seedlings (between July and October). The sugar industry equipment usually operates from January to March. In ratooning, sugarcane is harvested 12 months after harvesting in previous year, in spring planting it is harvested 12 months after planting, and in summer planting it is harvested in the sugar-production period 18 months after planting. About 63% of farmers manage a sugarcane area of less than 1 ha [48]. Cultivation of purple yam for confectionery products and as feed for meat calf husbandry is also active, and new pastures are being developed in some parts of the island. Pineapples were cultivated extensively on the island in the 1960s and 1970s, but are now grown only by some farmers and in private gardens. Flowers, bitter gourd, mango, and banana are cultivated in agricultural facilities such as greenhouses. Precise classification is required for various crop type fields on the island; this need also provides us an opportunity to evaluate the performance of classification of crop types and practices using high-frequency and high-spatial imagery from a satellite constellation.

Dove Satellite Image Pre-Processing
Planet Dove imagery was obtained from Planet Explorer (https://www.planet.com/explorer/) and processed as shown in Figure 2. The image product is delivered as a continuous, split-frame strip, with half-frames containing red-green-blue (RGB) and near-infrared (NIR) imagery, with an observation interval of approximately 0.5 s [49]. Cultivation of purple yam for confectionery products and as feed for meat calf husbandry is also active, and new pastures are being developed in some parts of the island. Pineapples were cultivated extensively on the island in the 1960s and 1970s, but are now grown only by some farmers and in private gardens. Flowers, bitter gourd, mango, and banana are cultivated in agricultural facilities such as greenhouses. Precise classification is required for various crop type fields on the island; this need also provides us an opportunity to evaluate the performance of classification of crop types and practices using high-frequency and high-spatial imagery from a satellite constellation.

Dove Satellite Image Pre-Processing
Planet Dove imagery was obtained from Planet Explorer (https://www.planet.com/explorer/) and processed as shown in Figure 2. The image product is delivered as a continuous, split-frame strip, with half-frames containing red-green-blue (RGB) and near-infrared (NIR) imagery, with an observation interval of approximately 0.5 s [49]. We used Dove ortho scenes (Planet Labs image product Level 3B: 4-Band PlanetScope Scene; spatial resolution = 3 m) from 1 July 2017 to 31 December 2019, which are orthorectified and radiometrically and sensor corrected. The geometric correction was based on digital elevation models with post-spacing between 30 and 90 m, with a root mean square error specification of <10 m [41]. We removed bad pixels using the unusable data mask provided with the products. We converted We used Dove ortho scenes (Planet Labs image product Level 3B: 4-Band PlanetScope Scene; spatial resolution = 3 m) from 1 July 2017 to 31 December 2019, which are orthorectified and radiometrically and sensor corrected. The geometric correction was based on digital elevation models with post-spacing between 30 and 90 m, with a root mean square error specification of <10 m [41]. We removed bad pixels using the unusable data mask provided with the products. We converted radiance to top-of-atmospheric reflectance using the scaling factor provided with the products. Dove satellites have been observed by two sensor types: Conventional type (PS2) and next-generation type (PS2.SD). Since these sensors have different relative spectral response curves, we transformed the next-generation-type values to match those of the conventional type using the coefficient described in the metadata.
Mismatches in the band-to-band registration were confirmed by visual inspection, and all images were re-registered by the zero-mean normalized cross-correlation method [50] with edge enhancing by a Laplacian filter; images with a shift of more than 10 pixels were excluded from the analysis. To maintain the continuity of pixel-based time series change, the registration between images was processed by the same method as the band-to-band re-registration, and was based on the reference image with an accurate position; images with a shift of more than 10 m were excluded from the analysis.
We performed simple dark pixel subtraction to obtain surface reflectance [51]. Assuming the presence of a pixel with zero reflectance (i.e., the radiance recorded by the sensor is solely from atmospheric scattering), the minimum pixel value was subtracted from those of all other pixels to remove the radiance derived from atmospheric scattering. Pixels from deep ocean areas were used for this correction. In addition, we normalized the value of each band according to Ono et al. (2002) [52] by dividing it by the summed value of all bands. Even on the same ground surface, the spectra observed by the satellite vary considerably because of the influence of the surface slope and atmospheric effects, but the shape of the spectra can be regarded as similar [52]. This characteristic allows suppression of the topographic and atmospheric effects by normalizing with the total sum of reflectance [53].
These pre-processed Planet Dove images were also used to evaluate the temporal variation of the pixels available for analysis. Specifically, the acquisition rates of analyzable pixels (i.e., not cloud pollution pixels) in the five days composite image were calculated.

Time Series of Vegetation Indices
The vegetation indices are highly correlated with the biomass of vegetation, and their temporal changes represent the growth pattern of crops [54]. NIR is the spectral region most responsive to changes in vegetation density [42]. Rouse et al. (1974) [55] and Tucker (1979) [56] developed the normalized difference vegetation index (NDVI), which is based on NIR and red bands and is widely used. The soil-adjusted vegetation index (SAVI) [57], modified soil-adjusted vegetation index (MSAVI) [58], and enhanced vegetation index (EVI) [59] were all developed on the basis of NDVI. However, the usefulness of evaluating the vegetation temporal change pattern by using indices based on the visible bands has also been demonstrated [60]. In Planet Dove satellite images with a split-frame strip, observation times differ between RGB and NIR bands, which may result in misregistration (see above). Thus, it is also necessary to consider vegetation indices that use visible wavelengths only. These include the visible atmospherically resistant index (VARI) [17], green-red ratio vegetation index (GRVI) [61,62], and normalized green ratio (GR) [60]. We calculated the following seven vegetation indices from pre-processed Dove satellite images to create a vegetation index time series: where ρ Blue , ρ Green , ρ Red and ρ NIR are Planet Dove blue, green, red, and near-infrared band reflectance, respectively. To suppress atmospheric noise, we applied the best index slope extraction (BISE) algorithm [63] and refined the time series data of the vegetation indices. The search range was set to 40 days for 10-day composite data (the missing periods were linearly interpolated), considering the observation frequency and agricultural calendar in the target area. To deal with unequal observation intervals in the 10-day composite data, we used an integrated algorithm based on BISE and the maximum value interpolated (MVI) algorithm [64] to refine the vegetation index time series data on a daily basis.

Ground Truthing
Agricultural field polygon data were obtained from the Ministry of Agriculture, Forestry and Fisheries, Japan (available at: https://www.maff.go.jp/j/tokei/porigon/) and used to draw the shape of each parcel of agricultural land on satellite images ( Figure 1). To gather reference data (i.e., training and validation data) for image classification, we conducted field surveys in the study area. Crop types and land cover conditions were recorded in each field and period (November 2018 and February, May, July, and October 2019). Main crop types and practices were recorded to the agricultural field polygon data. The following classes were used: (1) Sugarcane (ratoon), (2) sugarcane (spring planting), (3) sugarcane (summer planting in 2018), (4) sugarcane (summer planting in 2019), (5) pasture, (6) purple yam, (7) pineapple, (8) agricultural facility (e.g., vinyl greenhouse).
The use of available agricultural field boundary data improved the accuracy of classification because it allowed us to eliminate the effect of spectral variability within the agricultural field and the mixed-pixel effect [65]. Therefore, we used the agricultural field polygon data and extracted pure pixels corresponding only to the agricultural land. Sample sizes of reference data sets for each class were: (1) 36,902, (2) 12,671, (3) 17,743, (4) 2897, (5) 55,335, (6) 6628, (7) 874, and (8) 4281 pixels.

Classification and Accuracy Assessment
Considering the growing periods of pineapple and summer planting of sugarcane, the time series should span multiple years [37]. We used the time series data between 1 August 2017 and 30 November 2019 to generate a dataset for classification and validation. In the analysis of high-dimensional datasets, machine learning classification yields better results than the parametric classification algorithms do [66]. In particular, Random Forest (RF) [67] classification performs well without the need for complicated initial parameters, and is frequently used in land cover classification based on remote sensing data [68][69][70][71]. In general, machine learning algorithms, including RF, can be considered as black-box-type classifiers: The split rules for classification by RF are unknown [5], but RF can indicate the degrees of importance of variables for the general classification by the mean decrease in the Gini index [5]. The greater this decrease, the more important the variable is in the performance of the classification model. The relative importance of each variable was evaluated by normalizing to the sum of the mean decreases in the Gini index of all variables obtained during the classification model construction based on RF; this sum was considered to be 1. Such estimation is important for choosing appropriate dates for images that should be used to classify fields with high spatio-temporal variability. Therefore, we adopted RF to classify crop types and practices. For model development, tuning, and accuracy assessment, we used the add-on packages within version 3.6.1 of the 64-bit version of R [72]. The randomForest package [73] was used to build classifications with RF.
The classification parameters were set as follows. (1) The size of the training data was set to 500; this choice was based on our preliminary test with different sizes of the training data (n = 100, 200, 300, 400, 500). (2) The number of trees (ntree) was set to 500 after preliminary verification, in line with common practice [74]. (3) The other adjustable RF tuning parameter, mtry, controls the number of variables randomly considered at each split during tree building, and performance of RF is believed to be "somewhat sensitive" to this parameter [75]. The default value of mtry is the square root of p, where p equals the number of predictor variables within a dataset [73]. Generally, the default value is used [74], but in this study, mtry was tuned by the tuneRF function.
Classifications and validations of each vegetation index time series were conducted 10 times, each creating a confusion matrix. Training data were randomly sampled from the dataset (500 pixels in each class). Validation data were randomly sampled from a dataset that was not used as training data according to the area ratio of the class in the classification map (class of minimum area ratio was 30). Then the confusion matrix was used to calculate the evaluation indices: Overall accuracy (OA), producer's accuracy of each class, user's accuracy of each class, and Kappa [76]. Kappa is defined as: Remote Sens. 2020, 12, 2419 where P o is the proportion of cases correctly classified (i.e., overall accuracy) and P e is the expected proportion of cases correctly classified by chance [77]. The Kappa ranges from −1 to +1, where +1 indicates perfect agreement, and values of zero or less indicate a performance no better than random [78]. The Kappa has been traditionally used in accuracy assessment for classification map based on remote sensing. One such interpretation scale that has been widely used in remote sensing applications is that proposed by Landis and Koch (1977) [79] criteria [77]: Poor for negative Kappa values, slight for 0.00 < Kappa ≤ 0.20, fair for 0.20 < Kappa ≤ 0.40, moderate for 0.40 < Kappa ≤ 0.60, substantial for 0.60 < Kappa 0.80, and almost perfect for >0.80. Figure 3 shows the results of validation of the original band-to-band registration. Misregistration between RGB bands was observed in images with heavy cloud cover, images of the sea, and images taken by the next-generation sensor types (PS2.SD); this sensor took about 1% of all images analyzed. Misregistration was not observed in clear land images taken by the conventional sensor types (PS2). Therefore, no re-registration between RGB bands was performed for images taken with the PS2 sensors, and images taken by PS2.SD sensors were removed from the analysis. In many images, misregistration was seen between the RGB and NIR bands; these images were re-registered by matching the NIR band image to the RGB band image. Figure 4 shows the average NDVI time series of each crop type and practice in the field survey area created from the Planet Dove images. These series represent the crop growth patterns of each crop type and practice. Ratooning and spring planting of sugarcane had a similar crop growth pattern. In the summer planting of sugarcane, changes in the percentage of vegetation cover with time differed greatly from those in ratooning and spring planting. Pastures had high NDVI values throughout the year. Pineapple had a lower NDVI value throughout the year than did pastures. Agricultural facilities had lower NDVI values than any crop types and practices. The NDVI value of purple yam peaked in July, then decreased gradually, and was lowest in November. While the acquisition rates of the original NDVI values based on Planet Dove imagery were low from January to February and in June, spectra related to land cover information were usually obtained within 5 to 10 days for the rest of the periods.

Crop Growth Patterns Estimated by Planet Dove Imagery
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 21 Classifications and validations of each vegetation index time series were conducted 10 times, each creating a confusion matrix. Training data were randomly sampled from the dataset (500 pixels in each class). Validation data were randomly sampled from a dataset that was not used as training data according to the area ratio of the class in the classification map (class of minimum area ratio was 30). Then the confusion matrix was used to calculate the evaluation indices: Overall accuracy (OA), producer's accuracy of each class, user's accuracy of each class, and Κappa [76]. Kappa is defined as: where is the proportion of cases correctly classified (i.e., overall accuracy) and is the expected proportion of cases correctly classified by chance [77]. The Kappa ranges from −1 to +1, where +1 indicates perfect agreement, and values of zero or less indicate a performance no better than random [78]. The Kappa has been traditionally used in accuracy assessment for classification map based on remote sensing. One such interpretation scale that has been widely used in remote sensing applications is that proposed by Landis and Koch (1977) [79] criteria [77]: Poor for negative Kappa values, slight for 0.00 < Kappa ≤ 0.20, fair for 0.20 < Kappa ≤ 0.40, moderate for 0.40 < Kappa ≤ 0.60, substantial for 0.60 < Kappa ≤ 0.80, and almost perfect for >0.80. Figure 3 shows the results of validation of the original band-to-band registration. Misregistration between RGB bands was observed in images with heavy cloud cover, images of the sea, and images taken by the next-generation sensor types (PS2.SD); this sensor took about 1% of all images analyzed. Misregistration was not observed in clear land images taken by the conventional sensor types (PS2). Therefore, no re-registration between RGB bands was performed for images taken with the PS2 sensors, and images taken by PS2.SD sensors were removed from the analysis. In many images, misregistration was seen between the RGB and NIR bands; these images were re-registered by matching the NIR band image to the RGB band image.   Figure 4 shows the average NDVI time series of each crop type and practice in the field survey area created from the Planet Dove images. These series represent the crop growth patterns of each crop type and practice. Ratooning and spring planting of sugarcane had a similar crop growth pattern. In the summer planting of sugarcane, changes in the percentage of vegetation cover with time differed greatly from those in ratooning and spring planting. Pastures had high NDVI values throughout the year. Pineapple had a lower NDVI value throughout the year than did pastures. Agricultural facilities had lower NDVI values than any crop types and practices. The NDVI value of purple yam peaked in July, then decreased gradually, and was lowest in November. While the acquisition rates of the original NDVI values based on Planet Dove imagery were low from January to February and in June, spectra related to land cover information were usually obtained within 5 to 10 days for the rest of the periods. Positional deviation between Planet Dove bands by the zero-mean normalized cross-correlation method: Maximum displacement between red-green-blue (RGB) bands in the transverse direction (a) and vertical direction (b), and maximum displacement between RGB bands and the near-infrared (NIR) band in the transverse direction (c) and vertical direction (d). Since the validation was performed in the range of ±10 pixels in the XY directions of the image, a shift of 10 pixels or more was rounded to 10.     Figure 5 shows the classification accuracy of each vegetation index time series evaluated from 10 random samplings of training data. High classification accuracy was obtained with all vegetation indices ( Figure 5). All vegetation indices have an average OA exceeding 0.91, and the average Kappa exceeding 0.89. Vegetation indices based on both visible and NIR bands (NDVI, SAVI, MSAVI, and EVI) provided slightly better classification accuracy than those based on visible bands only (VARI, GRVI, and GR). This can be confirmed by the threshold based on the average value of classification by all vegetation indices (red line in Figure 5). No significant difference was found between NDVI and other vegetation indices based on visible and NIR bands. This suggested that the NIR band, which reflects vegetation better than visible wavelengths [80], could be effectively used after band-to-band re-registering. Because NDVI did not differ from other NIR-based vegetation indices, we recommend the traditional and versatile NDVI classification. We created a crop type and practice map using the NDVI time series (Figure 6). In this map, the class with the maximum average value of the class probabilities when the classification was conducted 10 times by RF. Class probability indicates the probability that a pixel belongs to each class when classified by RF. Even though classification in this map was on a pixel basis, the segments of crop types and practices were adequately represented in most of the agricultural parcels. We created a crop type and practice map using the NDVI time series (Figure 6). In this map, the class with the maximum average value of the class probabilities when the classification was conducted 10 times by RF. Class probability indicates the probability that a pixel belongs to each class when classified by RF. Even though classification in this map was on a pixel basis, the segments of crop types and practices were adequately represented in most of the agricultural parcels.

Mapping and Classification Accuracy of Crop Types and Practices Using NDVI Time Series
We compared the percentage of pure pixels that can be classified with the Planet Dove imagery with other representative medium-resolution satellite images. The ratio of pure pixel area to agricultural land area (agricultural field polygon data) was 85.2% in Planet Dove imagery (spatial resolution, 3 m), 55.5% in Sentinel-2 (10 m), and 11.3% in Landsat images (30 m). The presence or absence of target pixels in the polygon data was 99.8% in Planet Dove imagery, 90.3% in Sentinel-2 imagery, and 16.4% in Landsat imagery. Table 2 shows a confusion matrix that combines the results of 10 trials of classification using the NDVI time series. The classification accuracy of crop type and practice using NDVI time series was 0.93 for OA and 0.92 for Kappa ( Table 2). The accuracy of classification into sugarcane vs. others was 0.98 for OA and 0.97 for Kappa. The producer's accuracy and user's accuracy of sugarcane were 0.98 and 0.99 respectively, and the other's were 0.99 and 0.98 respectively. This result is superior to the published values of OA (0.84-0.96) and Kappa (0.59-0.87) obtained using remotely sensed satellite imagery (combination of Formosat/MS, ALOS/AVNIR2, SPOT/HRV and HRG [37], combination of Landsat sensors [81,82], MODIS [83]). However, it should be noted that the previous studies were not necessarily limited to crop type classification as in this study (for example, forests were included in other class [37]). We created a crop type and practice map using the NDVI time series (Figure 6). In this map, the class with the maximum average value of the class probabilities when the classification was conducted 10 times by RF. Class probability indicates the probability that a pixel belongs to each class when classified by RF. Even though classification in this map was on a pixel basis, the segments of crop types and practices were adequately represented in most of the agricultural parcels. We compared the percentage of pure pixels that can be classified with the Planet Dove imagery with other representative medium-resolution satellite images. The ratio of pure pixel area to agricultural land area (agricultural field polygon data) was 85.2% in Planet Dove imagery (spatial resolution, 3 m), 55.5% in Sentinel-2 (10 m), and 11.3% in Landsat images (30 m). The presence or  Figure 7 shows the time series changes in the importance of variables based on normalized reduction of Gini; the importance seemed to have a seasonal variation. The peak days of time series with 10-day intervals coincided with the timing of crop planting and harvesting ( Figure 4). Figure 8 shows the normalized mean decrease of Gini reassessed by classification using only NDVI images of 14 peak days. We found that classification accuracy varied with the number of images used for crop type and practice classification (Figure 9), where the variables used were added from the date with the most variable importance (i.e., normalized mean decrease Gini) shown in Figure 8. Table 2. Confusion matrix of the map shown in Figure 6. From the reference dataset based on the field survey, the reference data not used for training were used for validation. The number of samples in each class was determined by the area ratio, with the minimum set to 100.

Band-To-Band Registration of Planet Dove Satellite Images
The deliberately inexpensive sensor designs and commercial off-the-shelf components of small satellites worsen the consistency of acquired data quality in comparison with that of traditional large satellites [42]. Measurements over slightly mismatched areas from different spectral bands lead to much less accurate science data products when they are used together in heterogeneous areas [84]. Because the imager of PS2 had a split-frame visible (VIS) + NIR filter [41] and the RGB bands were recorded in the same frame based on the Bayer array [83], there was no misregistration between RGB bands. On the other hand, misregistration between RGB and NIR bands recorded in different frames was seen in many images (Figure 3c,d) of wide areas with contamination by small clouds. This bandto-band misregistration suggests that the registration was matched to the drifting cloud that could not be removed because of the difference in observation timing of approximately 0.5 s between the RGB bands and the NIR band. Images with misregistration between some RGB bands were taken by the PS2.SD sensors. Because PS2.SD has a four-band frame imager with a butcher-block filter providing blue, green, red, and NIR stripes [43], the observation timing differs slightly among all four bands. Therefore, misregistration may occur in all bands, even between RGB bands.

Crop Growth Patterns Estimated by Planet Dove Imagery
Ratooning and spring planting of sugarcane had similar crop growth patterns (Figure 4). The only major difference in land cover condition between ratooning and spring planting was whether or not there was tillage after harvesting and if it was completely exposed to bare land, and the NDVIbased crop growth pattern is considered to be similar because the growing periods were the same. The NDVI time series clearly represented the summer planting of sugarcane. NDVI time series of summer planting did not always maintain the lowest NDVI value from harvest to re-planting. This process is due to the growth of weeds or residual sugarcane stocks in the agricultural field left after harvest, and these being embedded in the soil by tillage just before planting. Pineapple had a lower

Band-To-Band Registration of Planet Dove Satellite Images
The deliberately inexpensive sensor designs and commercial off-the-shelf components of small satellites worsen the consistency of acquired data quality in comparison with that of traditional large satellites [42]. Measurements over slightly mismatched areas from different spectral bands lead to much less accurate science data products when they are used together in heterogeneous areas [84]. Because the imager of PS2 had a split-frame visible (VIS) + NIR filter [41] and the RGB bands were recorded in the same frame based on the Bayer array [83], there was no misregistration between RGB bands. On the other hand, misregistration between RGB and NIR bands recorded in different frames was seen in many images (Figure 3c,d) of wide areas with contamination by small clouds. This band-to-band misregistration suggests that the registration was matched to the drifting cloud that could not be removed because of the difference in observation timing of approximately 0.5 s between the RGB bands and the NIR band. Images with misregistration between some RGB bands were taken by the PS2.SD sensors. Because PS2.SD has a four-band frame imager with a butcher-block filter providing blue, green, red, and NIR stripes [43], the observation timing differs slightly among all four bands. Therefore, misregistration may occur in all bands, even between RGB bands.

Crop Growth Patterns Estimated by Planet Dove Imagery
Ratooning and spring planting of sugarcane had similar crop growth patterns (Figure 4). The only major difference in land cover condition between ratooning and spring planting was whether or not there was tillage after harvesting and if it was completely exposed to bare land, and the NDVI-based crop growth pattern is considered to be similar because the growing periods were the same. The NDVI time series clearly represented the summer planting of sugarcane. NDVI time series of summer planting did not always maintain the lowest NDVI value from harvest to re-planting. This process is due to the growth of weeds or residual sugarcane stocks in the agricultural field left after harvest, and these being embedded in the soil by tillage just before planting. Pineapple had a lower NDVI value than that of pastures, possibly because of the difference in vegetation cover ratio. On Kumejima Island, 50% of purple yam is planted in spring (March to May), 30% in summer (June to August), and 20% in autumn (September to November); it is harvested 5 or 6 months after planting. Since the purple yam fields are usually transposed in one or two years, the crop growth patterns of 2017 and 2018 may have been unclear due to a mixture of other crop type and practice patterns. These results adequately represent the agricultural calendar of the crop types and practices according to the interviews with farmers, and indicate that Planet Dove images can reveal crop growth patterns.
The NDVI time series generated from Planet Dove images contains noise components due to individual differences in sensors [85]. The Planet Dove's broadband sensor has a smaller NDVI range and a negative bias compared with the Landsat-8/OLI sensor. However, the relationship between NDVIs is linear, and the coefficient of determination is very high (top-of-atmosphere reflectance: 0.96; surface reflectance (6SV): 0.99). The variations in NDVI calculated from Planet Dove were within the range of approximately ±0.1 [85]. This means the variations had no considerable effect on our crop growth pattern observations because they were smaller than the range of temporal NDVI variations observed in this study ( Figure 4).

Mapping and Classification Accuracy of Crop Type and Practice Using NDVI Time Series
The classification accuracy of SAVI was slightly higher than that of NDVI ( Figure 5). SAVI was developed to reduce soil effects, which affect NDVI [57]. This may be the reason for the slightly higher classification accuracy of SAVI, because the amount of vegetation could be estimated more accurately than with NDVI on agricultural land with low vegetation cover. MSAVI and EVI had slightly lower classification accuracy than NDVI, possibly because the parameters in their equations were not properly estimated owing to the dispersion of reflectance observed by Planet Dove sensors. EVI uses the blue band in addition to the red and near-infrared bands [59], so the increased number of bands may have increased the variability of the index. Moreover, EVI used the parameters developed for the MODIS sensor, which may not be suitable for broadband Planet Dove sensors. The effects of these parameters on EVI should be clarified by more detailed research. Therefore, we used NDVI to create a crop type and practice map.
To the best of our knowledge, this study was the first to attempt to classify typical sugarcane crop practices in Japan using satellite remote sensing. The NDVI time series generated by Planet Dove satellite imagery classified sugarcane crop practices with good accuracy, exhibiting 0.93 for OA and 0.88 for Kappa, calculated based on Table 2. The slight confusion was mainly between ratoon and spring planting of sugarcane. There is no large difference in the growing season between these crop practices ( Figure 4). Variations in the NDVI values due to individual differences in sensors [85] could have obscured the slight differences in the NDVI time series of the ratoon and spring planting, which have similar growth periods, resulting in a slight confusion of the two classes shown in Table 2. The only difference is that for spring planting season, vegetation is completely removed in spring planting fields, whereas the soil of ratoon fields is covered with dry sugarcane leaves remaining after harvest (mulching); therefore, in the NDVI time series, spring planting showed a slightly lower value than that of the ratoon (Figure 4). Hence, the presence or absence of mulching in the early stages of growth is the key to discriminating between ratoon and spring planting. It is necessary to examine the sub-classification algorithm for ratoon and spring planting with a focus on the spring planting season. At the same time, the classification between crop practices that show similar crop growth patterns requires further efforts to suppress noise in the NDVI time series based on Planet Dove imagery caused by individual differences in sensors [85].
The results for the user's and producer's accuracy of each class were good ( Table 2). The producer's accuracy of ratoon was the lowest and ratoon was the most underestimated class ( Table 2), indicating that ratoon sugarcane fields had large spatio-temporal variability. The pineapple was the most overestimated class, and pastures and sugarcane (summer planting in 2018 and ratoon) were often misclassified as pineapple (Table 2). Pasture and ratoon misclassification may have been caused by the presence of vegetation throughout the year in unharvested areas of agricultural fields. The confusion may have also been caused by the overlap between summer planting of sugarcane and pineapple growth season; pineapples are harvested in July and August, plowed between August and October, and then replanted. Classification may be improved by extending the time series period from 3 to 5 years.
Only pure pixels, not mixed pixels with non-agricultural areas, were classified in this study. Hence, mixed pixels can be complemented by the neighboring pure pixels, and the coverage ratio of pure pixels within an agricultural field at the spatial resolution of Planet Dove imagery (99.8%) is very important in the classification of small-scale agricultural areas.

Assessment of the Optimal Dates and Required Observation Periods for Classification
The classification accuracy using 14 peak days (OA = 0.90, Kappa = 0.87) (Figure 9) was lower than that using all variables (OA = 0.93, Kappa = 0.92); yet good classification accuracy continues to be obtained. The number of peak days (14) almost coincided with the most frequent value (mtry = 15) of the number of variables determined by tuning during classification.
On a single date, Kappa showed fair (0.20-0.40), making good classification impossible ( Figure 9). Therefore, multi-date images are required for mapping subtropical agricultural areas. The most efficient approach for classification is to use five images in a two-year observation period ( Figure 9): Two images of the spring and summer sugarcane plantings in each year, and one image just before the harvesting in the second year. The first four images may contribute to the identification of sugarcane and other classes, and of sugarcane cropping practices. We infer that the last image would contribute mainly to distinguish sugarcane (in particular, ratoon and spring planting) from purple yam, which are all harvested mainly at this time of the year and have little vegetation cover. Using these images, excellent classification performance (Kappa >0.80; Figure 9) was obtained.
The images of the harvest season did not contribute to achieve high classification accuracy (Figures 7 and 8). This is supported by the poor classification accuracy for the similar crops based on images taken in the harvest season [37]. The classification accuracy in the study of Ishihara et al. [37], who used combined images from various medium spatial-resolution satellites (8-20 m) for an area similar to ours, was lower than in our study (OA = 0.64 < 0.93, Kappa = 0.47 < 0.92). This difference may be because they used only images taken in harvesting (January to March) and summer (July to September), but not in the spring planting season. Furthermore, these authors intentionally used images with low cloud cover, and may have missed the appropriate timing for classification because of the observation frequency. The greatest advantage of creating a time series by high-frequency earth observation using a satellite constellation is the improved probability of obtaining NDVI values at a time suitable for classification.

Conclusions
Using a time series of Planet Dove images with high spatio-temporal resolution, we investigated the performance of crop type and practice classification in subtropical small-scale agricultural areas with a focus on sugarcane crop practices, which have large spatio-temporal variations. These images, taken by nanosatellites, showed a mismatch in band-to-band registration, and it was necessary to re-register the NIR band. Cloud contamination remained in the vegetation index time series created from Planet Dove images, so it was necessary to reduce its effect by using a noise removal algorithm such as the integrated BISE and MVI algorithm. The number of classes was set to eight, including sugarcane subclassifications, and classification was implemented by a supervised RF algorithm. Although high classification performance was achieved for all vegetation index time series, we recommend using the commonly used NDVI, which includes the more effective NIR band. The NDVI showed high classification accuracy (OA = 0.93, Kappa = 0.92) in discriminating between crop types and practices across small fields in Japan, thus confirming the utility of high spatio-temporal resolution satellites like Planet Dove. Analysis of the importance of variables obtained by the RF algorithm revealed that the best dates for classification were the dates of planting and cropping. The most efficient way to classify crop types and practices is to use five image sets taken during a two-year observation period: Two image sets of the spring and summer sugarcane plantings in each year of a two-year observation period and one just before harvesting in the second year. These findings are applicable to other regions and land use and land cover classes, and would allow a high level of classification performance and effective use of images according to the subject.
Author Contributions: All authors contributed to the conceptualization and research design. A.S. conducted all analyses, including software preparation, satellite data collection, classification, validation, and visualization. A.S. and H.Y. conducted the field survey, collected reference data, and contributed to project administration. The original draft was written by A.S. and was reviewed and edited by H.Y. All authors have read and agreed to the published version of the manuscript.