Land suitability assessments for yield prediction of cassava using geospatial fuzzy expert systems and remote sensing

Computers and Electronics in Agriculture 166 (2019) 105018

Contents lists available at ScienceDirect

Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag

Land suitability assessments for yield prediction of cassava using geospatial fuzzy expert systems and remote sensing Riska Ayu Purnamasaria, Ryozo Noguchib, Tofael Ahamedb, a b

T

⁎

Graduate School of Life and Environmental Sciences, University of Tsukuba, Japan Faculty of Life and Environmental Sciences, University of Tsukuba, Japan

ARTICLE INFO

ABSTRACT

Keywords: Land suitability Cassava Yield prediction Fuzzy expert systems Remote sensing

Cassava has the potential to be a promising crop that can adapt to changing climatic conditions in Indonesia due to its low water requirement and drought tolerance. However, inappropriate land selection decisions limit cassava yields and increase production-related costs to farmers. As a root crop, yield prediction using vegetation indices and biophysical properties is essential to maximize the yield of cassava before harvesting. Therefore, the purpose of this research was to develop a yield prediction model based on suitable areas that assess with land suitability analysis (LSA). For LSA, the priority indicators were identified using a fuzzy expert system combined with a multicriteria decision method including ecological categories. Furthermore, the yield prediction method was developed using satellite remote sensing datasets. In this analysis, Sentinel-2 datasets were collected and analyzed in SNAP® and ArcGIS® environments. The multisource database of ecological criteria for cassava production was built using the fuzzy membership function. The results showed that 42.17% of the land area was highly suitable for cassava production. Then, in the highly suitable area, the yield prediction model was developed using the vegetation indices based on Sentinel-2 datasets with 10 m resolution for the accuracy assessment. The vegetation indices were used to predict cassava growth, biophysical condition, and phenology over the growing seasons. The NDVI, SAVI, IRECI, LAI, and fAPAR were used to develop the model for predicting cassava growth. The generated models were validated using regression analysis between observed and predicted yield. As the vegetation indices, NDVI showed higher accuracy in the yield prediction model (R2 = 0.62) compared to SAVI and IRECI. Meanwhile, LAI had a higher prediction accuracy (R2 = 0.70) than other biophysical properties, fAPAR. The combined model using NDVI, SAVI, IRECI, LAI, and fAPAR reported the highest accuracy (R2 = 0.77). The ground truth data were used for the evaluation of satellite remote sensing data in the comparison between the observed and predicted yields. This developed integrated model could be implemented for the management of land allocation and yield assessment in cassava production to ensure regional food security in Indonesia.

1. Background The productivity of cassava is one of the most important issues related to national food policy in Indonesia. Inappropriate management decisions and unsuitable land selection could limit productivity. Land suitability assessment is essential for current and future land use planning. The land suitability assessment can be performed using the multicriteria decision method. Several approaches of the multicriteria decision method (MCDM) using expert systems have been used to conduct land suitability assessments (Malczewski, 2006; Ferretti and Pomarico, 2013; Elsheikh et al., 2013). The MCDM becomes more appropriate with geospatial references. In recent years, computing technologies combined with GIS have enabled geospatial references using ⁎

land suitability evaluation based on FAO classes (FAO, 1976). To accommodate the complexity in assigning suitability classes, the fuzzy expert system combined with MCDM has the potential to determine the land suitability classification for cassava production. In MCDM, fuzzy set membership can be used in to standardize the criteria (Aydi et al., 2016; Feizizadeh and Blaschke, 2013; Romano et al., 2015). The fuzzy set theory allows for continuous factors to be modeled for a suitability assessment within a GIS analysis. However, most of those studies either used the MCDM technique or fuzzy set alone, resulting in a poorly handled weight of each factor or an inappropriate calculation of the suitability index. Furthermore, the MCDM with fuzzy set theory has the potential to reduce the subjectivity in the assessment of results. The integrated approach of GIS, fuzzy set,

Corresponding author. E-mail address: [email protected] (T. Ahamed).

https://doi.org/10.1016/j.compag.2019.105018 Received 10 October 2018; Received in revised form 3 July 2019; Accepted 19 September 2019 0168-1699/ © 2019 Elsevier B.V. All rights reserved.

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

Indonesia

Banten Province

Fig. 1. Land suitability analysis in the Banten province of Indonesia.

and MCDM has great potential to increase the effectiveness and accuracy of land suitability assessment for crop production (Aguilar-Rivera et al., 2018). Moreover, the effectiveness and accuracy of land suitability assessment can give more advantages for the yield prediction model. Because, a yield prediction model based on reference points such as highly suitable area may have higher accuracy than statistically predicted yields used in previous studies (Feizizadeh and Blaschke, 2013; Romano et al., 2015). Furthermore, considering cassava is a root crop, the prediction of yield using vegetation indices (VIs) and biophysical properties is important to maximize yield before harvesting. The most straightforward approach to estimating crop yields is to establish empirical relationships between ground-based yield measures and VIs measured on a single date or integrated over the growing season (Tucker, 1979). In this regard, satellite remote sensing and GIS applications for monitoring crops have the potential for timely assessments of changes in the growth and development of crops on regional scales (Lobell et al., 2015; Campos et al., 2018). The rapid development of spectral reflectance and remote sensing technologies has been used in various models to estimate crop yield using remote sensing applications (Azzari et al., 2017; Lobell et al., 2015; Sakamoto et al., 2013; Zabihi et al., 2015). For example, healthy crops and stress conditions are distinguished by the absorption of red radiation and the reflectance of NIR radiation (Tucker, 1979). Moreover, in the further research, the combination of the reflectance’s is more complex, not only red and near-infrared reflectance, but also green, red-edge and others. This combination also defined as vegetation indices (VIs) (Frampton et al., 2013). However, there is a lack of research on cassava, which is a root crop

that requires yield prediction based on the canopy and biomass development over the growing seasons. To the best of our knowledge, there is no research reported on cassava yield prediction at the suitable locations for cassava production. The yield prediction also helps in regional food security policy and production inventories to understand the availability of cassava. Thus, the objective of this research is to develop a model from various vegetation indices and biophysical properties derived from satellite datasets to predict the yield of cassava within highly suitable area using MCDM and fuzzy expert system. 2. Methodology 2.1. Study area Banten Province is located between 5°7′50″–7°1′1″S and 105°1′11″–106°7′12″E having 9663 km2 of total area and is the most western point of Java, approximately 90 km from Jakarta. Banten Province is strategically positioned as the area of land connecting Java and Sumatra Islands. Since the productivity of cassava is decreasing in Indonesia, the productivity of cassava in Banten province should be increased to support national food security (Fig. 1). Cassava production in the Banten province has the strong historical aspect. In some regions of Banten, cassava is an important alternative source of food, especially for traditional cuisine that is prepared for traditional events. Regarding this condition, the research for cassava is very important for the Banten province. The research conducted initially for the Banten Province and can be expended in other parts of Java Island. 2

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

Fig. 2. Framework of the research.

references and a literature review. In ArcGIS 10, there are seven kinds of fuzzy functions. In this study, based on the ecological criteria in the Banten province, we selected four functions: Large, Small, Gaussian and Linear functions. These functions produce continuous fuzzy classifications of standardized criteria. The Fuzzy Large transformation function is used when the larger input values are more likely to be a member of the set. The Fuzzy Small transformation function is used when the smaller input values are more likely to be a member of the set. The Gaussian function is useful if the membership functions represent the uncertainty in measurement. Whereas, the Fuzzy Linear transformation function applies a linear function between the user-specified minimum and maximum values (ESRI, CA, USA). In the Fuzzy Large, Small and Gaussian membership functions, the control point must include a midpoint (f2) and spread (f1) (Table 1). A midpoint is a point that has a 0.5 value of membership in Large, Small, and Gaussian functions and is determined by the user based on references (ESRI, CA, USA). The spread is generally assigned a number between 1 and 10. A higher spread value causes the fuzzy membership curve to become steeper. Meanwhile, the fuzzy linear transformation function applies a linear function between the minimum and maximum values. Anything below the minimum will be assigned a 0 (not a member) and anything above the maximum a 1 (a member) (Figs. 4 and 5). The equations for the fuzzy Large (Eq. (1)), Small (Eq. (2)), Gaussian (Eq. (3)), and Linear functions (Eq. (4)) are,

2.2. Framework of land suitability analysis As depicted in Fig. 2, the model was built over several stages. First, in the land suitability analysis, all the criteria of ecological potential of Banten province were standardized using fuzzy membership functions. Second, the highly suitable area for cassava was obtained in the ArcGIS® environment using weighted overlay tools and influenced by weights from MCDM. Then, cassava yield data points within the highly suitable areas were collected. These points were correlated to the corresponding vegetation indices (VIs) and biophysical properties to generate a regression model for the early prediction of cassava yield for the given fields. The Sentinel-2 satellite data sets were obtained to analyze the VIs and biophysical properties. 2.3. Land suitability analysis (LSA) 2.3.1. Reclassification by fuzzy membership function The fuzzy set theory allows for the concept of these continuous factors to be modeled in a suitability assessment within a GIS or spatial domain (Fig. 3). In a standard approach, membership within a set or class, is clearly and crisply defined as either in the class or not in the class (Bellman and Zadeh, 1970). In the present study, fuzzy membership classification is used to accommodate the high uncertainty of scoring methods in assigning the suitability classes. In this study, fuzzy membership functions were used in ArcGIS 10® (ESRI, CA, USA) for standardization. Several fuzzy functions were determined based on 3

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

1.2

1

Fuzzy Membership

Highly Suitable (S1) 0.8

Moderately Suitable (S2)

0.6

0.4

Marginally Suitable (S3)

0.2

Not Suitable (N) 0

0

20

40

60

80

100

120

Input Value Fig. 3. Reclassification by fuzzy membership functions.

on slopes between 0% and 15%. These slope ranges were determined to be highly suitable membership values. Based on the case and the criteria, we selected the small membership function that explained the small slope values indicating high suitability. The control point had a midpoint value of 25 and a spread value of 5.

1

µ (x ) = 1+

()

f1

x f2

(1)

1

µ (x ) = 1+ µ (x ) = e (

() x f2

f1

(2) 2

µ (x ) = f (x ) =

x b

0

a a a

1

2.3.2.2. Rainfall. Banten is characterized by a tropical climate and has different variability in some areas. This factor could affect the spatial pattern of the cassava yield in the whole region. The average temperature and rainfall levels are 27.4 °C and 1500–2000 mm/year, respectively (BPS, 2017). The smallest amount of water, cassava can grow with is approximately 400 mm/year, but higher yields can be obtained with a greater amount of water supplied to the area (FAO, 2013). For this condition, the Gaussian fuzzy function was used. The control point was 1500 mm/year with a spread of 10-4.

(3)

f 1 (X f 2) )

x a
(4)

2.3.2. Geographical extent and climatic factor 2.3.2.1. Slope. In Banten Province, most topography was classified as slopes between 0% and 45% in steepness. Cassava is easily cultivated Table 1 Suitability class by fuzzy membership for cassava production. Suitability classes

Fuzzy membership function

Criteria

S1

S2

S3

N

Mid-point

Spread

Slope (%)

0–8%

8–15%

15–25%

25–45%

25

5

Equation

1

µ (x ) = 1+

Rainfall (mm/year)

1500

1500–3000

< 1500

> 3000

1500

10

Elevation (m)

0–100

100–300

300–700

> 700

300

7

-4

µ (x ) =e (

0.0001 (X 1500)2)

1

µ (x ) = 1+

NDVI

1–0.7

0.7–0.5

0.5–0.3

0.3–0

0.53

5

µ (x ) = 1+

Distance from rivers (m)

1000

1000–1500

< 500

> 1000

1000

10

Distance from roads (m)

< 1000

1000–2000

2000–3000

> 3000

2000

5

-4

µ (x ) =e (

S1

S2

S3

N

Minimum

Maximum

LULC

Agri-culture dry field

Rice-field

Open field

Forest, Settle-ments, Waterbody

1

10

x 7 300 1

x 0.53

1

µ (x ) =

x 5 2000

0 (x ) = f (x ) =

Soil Type

Latosol

Rego-sol

Alu-vial

Podzol

1

5

0.0001 (X 1000)2)

1+

Criteria

x 5 25

x 1 1 10 1

1

25

0 µ (x ) = f (x ) =

x 1 1 25 1

1

4

x 1 < x < 10 x 10 x 1 < x < 25 x 25

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

1

0.8

Membership Values

Membership Values

1

0.6 0.4 0.2 0

0

10

20 30 Slope %

0.8 0.6 0.4 0.2 0

40

0

1000

1

1

0.8

0.8

0.6 0.4 0.2 0

0.6 0.4 0.2 0

0

200

400 600 Elevation (m)

800

0

0.2

1

1

0.8

0.8

0.6 0.4 0.2 0

500 1000 1500 Distance from river (m)

1

0.4 0.2 0

2000

0

1000 2000 Distance from road (m)

3000

(f) Distance from road membership function

1 Membershp Values

1

Membership Values

0.8

0.6

(e) Distance from river membership function

0.8 0.6 0.4 0.2 0

0.4 0.6 NDVI Value

(d) NDVI membership function

Membership Values

Membership Values

(c) Elevation membership function

0

5000

(b) Rainfall membership function

Membership Values

Membership Values

(a) Slope membership function

2000 3000 4000 Rainfall (mm/year)

0

2

4 6 LULC Values

8

0.8 0.6 0.4 0.2 0

10

0

5

10 15 Soil Values

20

(h) Soil membership function

(g) LULC membership function

Fig. 4. (a-h) Fuzzy membership function criteria.

5

25

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

(a) Soil

(b) Elevation

(c) Rainfall

(d) LULC

(f) Distance from road

(e) NDVI

(g) Distance from river

(h) Slope

Fig. 5. (a-h) Reclassification of the criteria.

6

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

determined as a linear function of fuzzy membership. The minimum membership value referred to a waterbody surface, and the maximum referred to a dry agricultural field surface. The LULC database was divided into ten classifications.

Table 2 Ground reference collection for cassava yield within highly suitable area. Field ID

Village

Latitude

Longitude

Yield (ton/ha)

31 33 57 58 61 63 65 66 69 70 71 76 79 80 82 83 84 85 86 91 92 93 94 98 101

Cililitan Pasireurih Bojongjuruh Sumberwaras Bejod Ciginggang Kramat Jaya Cijaku Lebak tipar Peucang pari Cirinten Panyaungan Pabubulan Lebak Cimayang Cimayang Bojong Cisimet Cisimet Raya Sukamarga Padasuka Bojongleles Bojong cae Sangiang Cileles

−6.51 −6.40 −6.53 −67.56 −6.76 −6.55 −6.62 −6.70 −6.96 −6.72 −6.67 −6.91 −6.92 −6.96 −6.60 −6.59 −6.58 −6.57 −6.56 −6.53 −6.31 −6.36 −6.31 −6.38 −6.48

105.98 105.94 105.99 105.99 105.94 106.02 106.07 106.06 106.34 106.13 106.19 106.16 106.29 106.34 106.17 106.17 106.20 106.24 106.31 106.32 106.20 106.21 106.25 106.41 106.11

12.00 13.00 23.00 17.30 16.00 21.00 9.25 25.00 15.00 15.00 16.00 20.00 22.00 6.00 5.00 22.00 10.26 8.00 17.00 23.00 3.33 5.80 16.00 15.00 4.00

2.3.2.8. Soil. The soil plays a significant role in crop growth and cassava yield (Howeler, 1991). The major soil types found in Banten are alluvial, red regosol, red-yellow podzolic, and latosol soils. According to Wargiono and Sudaryanto (2000), latosol areas are optimal for cultivating cassava. Latosol soils have good physical properties and are deep and tolerant to erosion. However, podzols include low levels of organic matter and tend to erode easily. Just as LULC was classified, the soil was classified with the linear function, with Latosol as the maximum value and podzols as the minimum value. 2.3.3. Weighted linear combination (WLC) The WLC method was used to calculate the suitability index (S) for the area, Wi is the relative importance weight of criterion, Ri is the standardized value of the area under criterion i and n is the total number of criteria. The suitability score for each of the land units was calculated using the following equation: n

S=

Wi × Ri i=1

(5)

In the FAO's framework for land evaluation, the land index classification was designated as suitable (S) or not suitable (N). These suitability classes can then be further subdivided as required. The suitability classification aimed to show the suitability of each land unit for cassava production. In this section, three classes (S1, S2, and S3) were used to distinguish land that is highly suitable, moderately suitable and marginally suitable for cassava production. For the verification of the F-MCDM model, we used the regression analysis using ground truth data collected from the cassava yield within all suitable areas in the Banten province. The correlation of the Land Suitability Index (LSI) with actual yield was developed (Azzari et al., 2017; Braimoh et al., 2004; Sys et al., 1991).

2.3.2.3. Elevation. In Indonesia, most cassava-growing areas are located in the lowland humid and subhumid tropics (Heumann et al., 2011). In some areas, cassava can be grown in hilly or mountainous terrain. In Banten Province, elevation ranges from 0 to 700 m. Most of the area is suitable for cassava production, although the optimal elevation for cassava production is approximately 62.5–137.5 m. In this case, the small function was fitted with a control point midpoint of 300 m and a spread of 7. 2.3.2.4. Vegetation index. The vegetation canopy is an important factor that could prevent soil erosion during cassava production. The NDVI is a vegetation index that is correlated with several important biophysical properties and is used to generate different crop indices (Elhag, 2014). Higher vegetation canopy indicated more highly suitable areas for cassava production. In this circumstance, the large function was used with 0.53 as the midpoint and a spread of 5 for the control point.

2.4. Yield prediction 2.4.1. Sentinel-2 imagery All 24-available satellite images from Sentinel-2A over Banten province taken between March 2017 and February 2018 were downloaded from the USGS earth explorer internet resource. The cloud masking and buffering were performed using idepix toolbox of Sentinel Application Platform (SNAP) 5.0. For calculating the vegetation indices, several bands at a 10-meter resolution were used. The bands included in this analysis were Band 4 (red), Band 8 (NIR) and three new bands in the red-edge region, Band 5, Band 6 and Band 7, which are centered at 705, 740 and 783 nm, respectively.

2.3.2.5. Distance from river. In Banten Province, there are four wide river regions (Cibaliung, Ciliman, Cidanau, and Cisadane) that include over 102 watershed areas. The physical factors associated with water supply were used to determine suitability levels for cassava production. Cassava can be planted on areas that are located not too far or too close to water resources. The Gaussian function was chosen when the midpoint referred to higher suitability. The control point was 1000 m for its midpoint and 10−4 for its spread.

2.4.2. Vegetation indices 2.4.2.1. NDVI. The NDVI was proposed by Rouse et al. (1973) and it has become the most widely used indicator for studying vegetation canopy, crop production, green biomass, chlorophyll content, and canopy water stress (Vrieling et al., 2011). The NDVI utilizes two essential wavelength bands: the red and near-infrared (NIR) bands. The combination of Band 4 and Band 8 were used, corresponding to red and NIR reflection, respectively. The relationship between NDVI and yield has been described in various experiments. The NDVI was calculated as follows:

2.3.2.6. Distance from road. In selecting areas suitable for cassava production, the distance from roads must be considered because such distances affect transportation costs for supply processes. Shorter distances between fields and roads facilitate access to transportation infrastructure and link farmers and farming activities to marketing channels. This means that the small function was fitted for this criterion. The control point was set with a midpoint of 2000 m and a spread of 5.

NDVI =

2.3.2.7. Land use and land cover. Land use and land cover (LULC) data files describe the vegetation, water, natural surfaces, and cultural features of a land surface (Akıncı et al., 2013). Rice fields covered most lands in the Banten province. The classifications of LULC was

RNIR RRed RNIR + RRed

(6)

2.4.2.2. SAVI. The calculation of vegetation indices is influenced by 7

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

0.7 NDVI

0.6

SAVI

Vegetation Indices

0.5 0.4 0.3 0.2 0.1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

Number of fields

(a) Variability of NDVI and SAVI in Banten province 30

Yield (Ton/Ha)

25 20 15 10 5 0

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Number of fields

(b) Variability of yield of cassava in Banten province Fig. 6. Variability of cassava fields in the same growing season.

soil background conditions in some cases (Huete, 1988; Gilabert et al., 2002). Huete (1988) used a soil adjustment factor, L, to reduce the soil background effect. The proposed constant adjustment factor (L = 0.5) is referred to as the soil-adjusted vegetation index (SAVI):

R RRed SAVI = (1 + L) x NIR RNIR + RRed

where RRE1 (Red Edge 1) is Band 5 (705 nm) and RRE2 (Red Edge 2) is Band 6 (740 nm). 2.4.3. Biophysical properties 2.4.3.1. LAI. The biophysical monitoring of essential plant properties, such as chlorophyll, nitrogen, LAI, leaf water content and crop health, is very important. The leaf area index (LAI) is an important biophysical variable for agricultural land monitoring and modeling studies. Since it plays a important role in ecological processes, LAI retrieval maps from satellites have been used in many land observation studies. LAI retrieval from Sentinel-2 was calculated using an algorithm provided by the SNAP® biophysical operations toolbox.

(7)

2.4.2.3. IRECI. IRECI (Inverted Red Edge Chlorophyll Index) utilized all red-edge bands provided by Sentinel-2. The combined bands have substantial effects on chlorophyll absorption and internal leaf scattering. The IRECI also has potential as the optimal red edge index for evaluating the biomass status using Sentinel-2A imagery (Castillo et al., 2017).

IRECI =

(RNIR

2.4.3.2. fAPAR. Some research found that total biomass production is closely related to the fraction of photosynthetically active radiation (fAPAR). fAPAR is absorbed in the canopy over the growing season. Estimations of fAPAR are often derived from Vis. In the SNAP®

RRed)

RRE 2 RRE1

(8) 8

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

0.7 0.6

Vegetation Indices

0.5 0.4

NDVI SAVI

0.3

IRECI

0.2 0.1 0

06/03/2017

14/06/2017

24/07/2017

01/11/2017

26/12/2017

19/02/2018

(a) Vegetation Index values for cassava during growing season 2.5

Biophysical Properties

2

1.5 LAI fAPAR

1

0.5

0

06/03/2017

14/06/2017

24/07/2017

01/11/2017

26/12/2017

19/02/2018

(b) Biophysical properties for cassava during growing season Fig. 7. The optimum growing season based on VIs (a) and biophysical properties (b).

software, the calculation of fAPAR is also carried out using the biophysical operations toolbox. 2.4.4. Ground reference data The location points of the cassava fields were collected by GPS Garmin eTrex (Garmin Inc. KA, USA) around Banten province in 2016–2017. Yield estimates were evaluated using the ground-truth yield data based on the farm surveys at the cassava field locations. The harvested time was recorded for each cassava field. The field location and yield data collected on the highly suitable area listed in Table 2. In the Banten province, the farmers have different times to start to plant cassava. The variability of the fields within the various growing seasons is shown in Figs. 6 and 7.

Fig. 8. Suitable area for cassava production.

9

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

biophysical properties had higher prediction accuracy (R2 = 0.70) than fAPAR. The combined model using NDVI, SAVI, IRECI, LAI, and fAPAR reported the highest accuracy (R2 = 0.77). The combination model was used to generate the yield prediction map. The ground truth data were used to evaluate the satellite remote sensing data between the observed and predicted yields. Finally, the yield map was generated to show the variability of cassava yield in the various suitable area in Banten province (Fig. 12). The higher yield prediction was observed in the southern part of the Banten province. This result was consistent regarding the highly suitable area that also located in the same region (Fig. 8). The prediction formula was applied for estimating the yield before harvesting around December 2017. For comparison, the prediction was performed in the cassava fields within highly suitable and moderately suitable areas. In a highly suitable area, the average yield for cassava production was 12.21 ton/ha. This prediction was much higher than the prediction of cassava yield in the moderately suitable area (8.56 ton/ha) (Tables 5 and 6).

Table 3 Suitable areas for cassava production. Suitability classes

Fuzzy-MCDM

Highly suitable Moderately suitable Marginally suitable Not suitable

Percentage area (%)

Area (ha)

42.17 43.10 6.25 8.47

407,488 416,475 60,393 81,845

3. Results 3.1. Land suitability analysis The results showed that 42.17% of land area was highly suitable for the F-MCDM model (Fig. 8 and Table 3). In this regard, F-MCDM land suitability maps reflected interaction between the fuzzy membership function values and their weights. To evaluate the accuracy of the methods, the correlation between land suitability index (LSI) and actual yield of cassava in various suitable area was obtained (Fig. 9). This regression shows that 55% of cassava fields in the study area were located in suitable lands for cassava production.

4. Discussion In the land suitability analysis, F-MCDM could be good approach when the criteria were more complex and required a continuous value. The F-MCDM method was validated by the fact that approximately 55% of cassava fields were located in a proper suitable location. In the other way we can assume the higher yield that located in higher suitability area and lower yield that located in less suitable area is around 55%. The suitable areas were mostly located in the southern part of the Banten Province. Therefore, farmers in the Banten province required highly suitable land to cultivate cassava. In further assessments, the highly suitable areas were determined to predict the yield of cassava production. Moreover, to help the farmer obtain more yield in the optimum window for harvesting, the prediction of cassava yield was needed. The best time to predict the cassava yield was observed between yield and vegetation indices in the green-up season around July 2017. This result indicated that cassava yield can be predicted in July, five months before schedule harvest in December or February (Fig. 11). The vegetation phenology analysis and biophysical properties of cassava plants has the potential for estimating yield prediction as root crops with good accuracy in highly suitable areas. In general, all indicators explain the vegetation canopy coverage. A good prediction model was obtained from the combination of vegetation indices (NDVI, SAVI, and IRECI) and biophysical properties (LAI, fAPAR). NDVI, LAI,

3.2. Yield prediction using vegetation indices and biophysical properties The yield data from the field were compared with data retrieved from Sentinel-2. The analysis was conducted to find the relationship between yield, VIs, and biophysical values at each growing stage. This study investigated the use of VIs and biophysical properties in a more detailed field area located in a highly suitable location based on a ground-truth survey. The results of the study show that the VIs and biophysical properties during the 24 July 2017 growing season had slightly stronger correlations than the other (Figs. 10, 11, and Table 4). The yield prediction model was developed using the vegetation index retrieved from the Sentinel-2 datasets (10 m resolution). The vegetation indices were used to predict cassava growth, biophysical condition, and phenology over the growing seasons. The NDVI, SAVI, IRECI, LAI, and fAPAR were used to develop the predictive model for cassava growth. The generated models were validated using regression analysis to estimate the variations among the observed and predicted yields. As vegetation indices NDVI showed higher accuracy in the yield prediction model (R2 = 0.62) compared to SAVI and IRECI. The LAI as

30

Observed Yield (ton/ha)

25 y = 0.1348x + 3.5722 R² = 0.55

20 15 10 5 0

0

20

40

60 Land Suitability Index

80

Fig. 9. Verification of the F-MCDM land suitability model. 10

100

120

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

Fig. 10. Relationship between VIs, biophysical properties and the actual yield of cassava.

0.6 0.5 NDVI

0.4

R2

SAVI 0.3

IRECI LAI

0.2

fAPAR

0.1 0

06/03/2017

14/06/2017

24/07/2017

01/11/2017

26/12/2017

Fig. 11. Growing season for cassava in Indonesia.

11

19/02/2018

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

5. Conclusions

Table 4 The yield prediction models. Vegetation indices

R2

Yield expression

NDVI SAVI IRECI

0.62 0.48 0.37

Yield = 39.234*NDVI − 3.6631 Yield = 49.895*SAVI + 0.4108 Yield = 23.831*IRECI + 6.1878

Biophysical properties LAI fAPAR All Combination

R2

Formula

0.70 0.27 0.77

Yield = 9.4269*LAI + 0.6379 Yield = 21.636*fAPAR + 2.19 Yield = 33*NDVI – 54*SAVI + 8*IRECI + 2.2*LAI + 4.8*fAPAR

Inappropriate decisions relating to land selection limit the productivity of cassava cultivation and increase the associated costs to farmers. This research was conducted to develop a land suitability model on the provincial scale to find the best suitable areas for cassava production. Fuzzy set theory is advantageous for the standardization of criteria using fuzzy membership functions. In ArcGIS10®, there were 7 fuzzy functions classified into Large, Small, Gaussian and Linear functions. In other words, this research tried to demonstrate the sufficiency of fuzzy functions for cassava production area selection, which was based on the midpoint and spread defined control points. After this stage, the final suitability maps were determined using the weighted overlay method. A suitable land map based on fuzzy-based MCDM was generated. In further assessments, the highly suitable areas were determined for predicting yield of cassava production. The spectral bands of Sentinel-2 satellite imagery, vegetation indices, and biophysical properties were used as input parameters to generate cassava yield prediction models. The study was carried out in cassava fields that were located in highly suitable locations for cassava production in the Banten province in Indonesia. The regression models were generated using the vegetative indices and biophysical properties. The yield prediction model was developed using the vegetation indices and biophysical properties from the Sentinel-2 datasets with 10 m resolution. The vegetation indices (NDVI, SAVI, IRECI) and biophysical properties (LAI and fAPAR) were used to predict cassava growth, biophysical condition, and phenology over the growing seasons. The ground truth data were used to evaluate of satellite remote sensing data between the estimated yield within highly suitable and moderately suitable areas. From this result, the estimation of cassava yield was better when the field was located in a highly suitable area. The developed decision support system was integrated with the expert system, GIS, and Sentinel-2 satellite datasets to evaluate land suitability and predict cassava yield from canopy biomass and soil adjusted indices. The developed model can be used at the regional and country levels in land suitability assessment and yield estimation of cassava to maximize production. In further research, machine learning will be added to validate the yield prediction model. The developed integrated decision support system model can be used to make recommendations to policy

Fig. 12. Cassava yield prediction map (ton/ha).

fAPAR related with chlorophyll content. Some related indices like SAVI and IRECI can be used to extract information on protein content and crop yield prediction before the harvest time.

Table 5 Yield prediction for cassava in highly suitable area. Field ID

Village

Latitude

Longitude

Observed yield (ton/ha)

NDVI

SAVI

LAI

fAPAR

IRECI

Predicted yield (ton/ha)

31 33 57 58 61 63 65 66 69 70 71 76 79 80 82 83 84 85 86 91 92 93 94 98 101 Average

Cililitan Pasireurih Bojongjuruh Sumberwaras Bejod Ciginggang Kramat Jaya Cijaku Lebak tipar Peucang pari Cirinten Panyaungan Pabubulan Lebak Cimayang Cimayang Bojong Cisimet Cisimet Raya Sukamarga Padasuka Bojongleles Bojong cae Sangiang Cileles

−6.51 −6.40 −6.53 −68.00 −6.76 −6.55 −6.62 −6.70 −6.96 −6.72 −6.67 −6.91 −6.92 −6.96 −6.60 −6.59 −6.58 −6.57 −6.56 −6.53 −6.31 −6.36 −6.31 −6.38 −6.48

105.98 105.94 105.99 105.99 105.94 106.02 106.07 106.06 106.34 106.13 106.19 106.16 106.29 106.34 106.17 106.17 106.20 106.24 106.31 106.32 106.20 106.21 106.25 106.41 106.11

12.00 13.00 23.00 17.30 16.00 21.00 9.25 25.00 15.00 15.00 16.00 20.00 22.00 6.00 5.00 22.00 10.26 8.00 17.00 23.00 3.33 5.80 16.00 15.00 4.00

0.58 0.62 0.60 0.78 0.80 0.76 0.57 0.34 0.45 0.64 0.37 0.73 0.75 0.40 0.31 0.68 0.73 0.70 0.71 0.79 0.82 0.76 0.61 0.66 0.77

0.37 0.38 0.41 0.54 0.54 0.54 0.32 0.23 0.14 0.42 0.21 0.47 0.49 0.26 0.19 0.44 0.47 0.42 0.40 0.46 0.57 0.48 0.37 0.41 0.47

2.08 2.77 1.81 3.76 3.93 3.97 1.96 1.11 1.18 1.94 1.08 1.89 2.14 0.85 0.63 2.19 2.64 2.44 0.91 2.75 3.70 2.57 1.53 2.00 2.74

0.65 0.62 0.61 0.83 0.83 0.84 0.62 0.52 0.35 0.64 0.44 0.71 0.71 0.39 0.28 0.67 0.75 0.70 0.60 0.72 0.83 0.72 0.54 0.64 0.74

0.49 0.39 0.51 1.12 1.00 1.08 0.40 0.26 0.10 0.59 0.22 0.45 0.86 0.24 0.16 0.65 0.94 0.66 0.56 0.79 1.24 0.77 0.47 0.53 0.81

10.72 12.51 8.90 18.18 18.18 17.36 12.05 5.87 12.01 10.60 6.71 9.85 13.42 5.13 3.95 11.88 15.40 14.31 10.84 17.09 18.34 14.58 9.69 11.71 15.86 12.21

12

Computers and Electronics in Agriculture 166 (2019) 105018

R. Ayu Purnamasari, et al.

Table 6 Yield prediction for cassava in moderately suitable area. Field ID

Village

Latitude

Longitude

Observed yield (ton/ha)

NDVI

SAVI

IRECI

LAI

fAPAR

Predicted yield (ton/ha)

37 40 87 88 100 Average

Menes Citalahab Ciminyak Sobang Cikulur

−6.39 −6.38 −6.54 −6.62 −6.40

105.93 106.11 106.31 106.30 106.18

13.00 12.00 21.00 10.62 12.00

0.60 0.35 0.62 0.48 0.52

0.23 0.33 0.46 0.57 0.39

0.25 0.40 0.79 1.24 0.62

0.88 1.08 2.75 3.70 1.94

0.36 0.61 0.72 0.83 0.66

13.02 2.53 11.80 7.03 8.42 8.56

makers and planners in Indonesia to increase the practices of producing cassava in highly suitable areas. Furthermore, the developed model can be employed for yield prediction to develop an inventory planning system to ensure regional food security.

Elsheikh, R., Shariff, A.R.B.M., Amiri, F., Ahmad, N.B., Balasundram, S.K., Soom, M.A.M., 2013. Agriculture Land Suitability Evaluator (ALSE): A decision and planning support tool for tropical and subtropical crops. Comput. Electron. Agric. 93, 98–110. https:// doi.org/10.1016/j.compag.2013.02.003. FAO, 1976. A Framework for Land Evaluation, first ed. FAO, Rome, Italy. FAO, 2013. Save and Grow: Cassava A Guide to Sustainable Production Intensification. FAO, Rome, Italy. Frampton, W.J., Dash, J., Watmough, G., Milton, E.J., 2013. Evaluating the capabilities of Sentinel-2 for quantitative estimation of biophysical variables in vegetation. ISPRS J. Photogramm. Remote Sens. 82, 83–92. https://doi.org/10.1016/j.isprsjprs.2013.04. 007. Feizizadeh, B., Blaschke, T., 2013. GIS-multicriteria decision analysis for landslide susceptibility mapping: comparing three methods for the Urmia lake basin, Iran. Nat. Hazards 65 (3), 2105–2128. https://doi.org/10.1007/s11069-012-0463-3. Ferretti, V., Pomarico, S., 2013. An integrated approach for studying the land suitability for ecological corridors through spatial multicriteria evaluations. Environ. Dev. Sustain. 15 (3), 859–885. https://doi.org/10.1007/s10668-012-9400-6. Gilabert, M.A., González-Piqueras, J., Garcıa-Haro, F.J., Meliá, J., 2002. A generalized soil-adjusted vegetation index. Remote Sens. Environ. 82 (2–3), 303–310. Huete, A.R, 1988. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 25 (3), 295–309. https://doi.org/10.1016/0034-4257(88)90106-X. Heumann, B.W., Walsh, S.J., McDaniel, P.M., 2011. Assessing the application of a geographic presence-only model for land suitability mapping. Ecol. Inf. 6 (5), 257–269. https://doi.org/10.1016/j.ecoinf.2011.04.004. Howeler, R.H., 1991. Long-term effect of cassava cultivation on soil productivity. Field Crops Res. 26 (1), 1–18. https://doi.org/10.1016/0378-4290(91)90053-X. Lobell, D.B., Thau, D., Seifert, C., Engle, E., Little, B., 2015. A scalable satellite-based crop yield mapper. Remote Sens. Environ. 164, 324–333. https://doi.org/10.1016/j.rse. 2015.04.021. Malczewski, Jacek, 2006. GIS-based multicriteria decision analysis: a survey of the literature. Int. J. Geographical Inf. Sci. 20 (7), 703–726. https://doi.org/10.1080/ 13658810600661508. Romano, G., Dal Sasso, P., Liuzzi, G.T., Gentile, F., 2015. Multi-criteria decision analysis for land suitability mapping in a rural area of Southern Italy. Land Use Policy 48, 131–143. https://doi.org/10.1016/j.landusepol.2015.05.013. Rouse, J.W., Haas, R.H., Schell, J.A., Deering, D.W., 1973. Monitoring vegetation systems in the Great Plains with ERTS, Third ERTS Symposium, NASA SP-351 I, 309–317. Sakamoto, T., Gitelson, A.A., Arkebauer, T.J., 2013. MODIS-based corn grain yield estimation model incorporating crop phenology information. Remote Sens. Environ. 131, 215–231. Sys, C., Van Ranst, E., Debaveye, I.J., 1991. Land evaluation. Part I: Principles In Land Evaluation and Crop Production Calculations. General Administration for Development Cooperation, Agricultural publication-No. 7, Brussels-Belgium, 274. Tucker, C.J., 1979. Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ. 8 (2), 127–150. Vrieling, A., de Beurs, K.M., Brown, M.E., 2011. Variability of African farming systems from phenological analysis of NDVI time series. Clim. Change 109 (3–4), 455–477. https://doi.org/10.1007/s10584-011-0049-1. Wargiono, J., Sudaryanto, B., 2000. Cassava leaves and forage crops for ruminant feed in the establishment of sustainable cassava farming system in Indonesia. In National Workshop-Seminar on Sustainable Livestock Production on Local Feed Resources, pp. 496–503. http://ciat.library.ciat.cgiar.org/Articulos_Ciat/proceedings_workshop_02/ 496.pdf. Zabihi, H., Ahmad, A., Vogeler, I., Said, M.N., Golmohammadi, M., Golein, B., Nilashi, M., 2015. Land suitability procedure for sustainable citrus planning using the application of the analytical network process approach and GIS. Comput. Electron. Agric. 117, 114–126. https://doi.org/10.1016/j.compag.2015.07.014.

Acknowledgments We would like to thank the University of Tsukuba to support this research to develop the multi-criteria modeling for land suitability analysis for cassava production in Indonesia. We also express our sincere thanks to the, Indonesian Geospatial Agency, the United States Geological Survey (USGS) and European Space Agency (ESA) for geographical and satellite data information. We sincerely thank the Indonesia Endowment Fund for Education (LPDP) for providing scholarship to continue this research in Japan. We also expressed our gratitude’s to Indonesian experts and field surveyors to participate in this research. References Aguilar-Rivera, N., Algara-Siller, M., Olvera-Vargas, L.A., Michel-Cuello, C., 2018. Land management in Mexican sugarcane crop fields. Land Use Policy 78, 763–780. https:// doi.org/10.1016/j.landusepol.2018.07.034. Akıncı, H., Özalp, A.Y., Turgut, B., 2013. Agricultural land use suitability analysis using GIS and AHP technique. Comput. Electron. Agric. 97, 71–82. https://doi.org/10. 1016/j.compag.2013.07.006. Aydi, A., Abichou, T., Nasr, I.H., Louati, M., Zairi, M., 2016. Assessment of land suitability for olive mill wastewater disposal site selection by integrating fuzzy logic, AHP, and WLC in a GIS. Environ. Monit. Assess. 188 (1), 59. https://doi.org/10.1007/s10661015-5076-3. Azzari, G., Jain, M., Lobell, D.B., 2017. Towards fine resolution global maps of crop yields: Testing multiple methods and satellites in three countries. Remote Sens. Environ. 202, 129–141. https://doi.org/10.1016/j.rse.2017.04.014. Bellman, R.E., Zadeh, L.A., 1970. Decision-making in a fuzzy environment. Manage. Sci. 17 (4), B-141. https://doi.org/10.1287/mnsc.17.4.B141. BPS, 2017. Regional Statistics Banten Province 2017. Badan Pusat Statistik (The Central Statistics Agency of Indonesia). BPS Banten Province. Braimoh, A.K., Vlek, P.L., Stein, A., 2004. Land evaluation for maize based on fuzzy set and interpolation. Environ. Manage. 33 (2), 226–238. https://doi.org/10.1007/ s00267-003-0171-6. Campos, I., González-Gómez, L., Villodre, J., Calera, M., Campoy, J., Jiménez, N., Calera, A., 2018. Mapping within-field variability in wheat yield and biomass using remote sensing vegetation indices. Precis. Agric. 1–23. https://doi.org/10.1007/s11119-0189596-z. Castillo, J.A.A., Apan, A.A., Maraseni, T.N., Salmo III, S.G., 2017. Estimation and mapping of above-ground biomass of mangrove forests and their replacement land uses in the Philippines using Sentinel imagery. ISPRS J. Photogramm. Remote Sens. 134, 70–85. https://doi.org/10.1016/j.isprsjprs.2017.10.016. Elhag, M., 2014. Sensitivity analysis assessment of remotely based vegetation indices to improve water resources management. Environ. Dev. Sustain. 16 (6), 1209–1222. https://doi.org/10.1007/s10668-014-9522-0.

13