Integrating Remote Sensing and a Markov-FLUS Model to Simulate Future Land Use Changes in Hokkaido, Japan

: As the second largest island in Japan, Hokkaido provides precious land resources for the Japanese people. Meanwhile, as the food base of Japan, the gradual decrease of the agricultural population and more intensive agricultural practices on Hokkaido have led its arable land use to change year by year, which has also caused changes to the whole land use pattern of the entire island of Hokkaido. To realize the sustainable use of land resources in Hokkaido, past and future changes in land use patterns must be investigated, and target-based land use planning suggestions should be given on this basis. This study uses remote sensing and GIS technology to analyze the temporal and spatial changes of land use in Hokkaido during the past two decades. The types of land use include cultivated land, forest, waterbody, construction, grassland, and others, by using the satellite images of the Landsat images in 2000, 2010, and 2019 to achieve this goal to make classiﬁcation. In addition, this study used the coupled Markov-FLUS model to simulate and analyze the land use changes in three different scenarios in Hokkaido in the next 20 years. Scenario-based situational analysis shows that the cultivated land in Hokkaido will drop by about 25% in 2040 under the natural development scenario (ND), while the cultivated land area in Hokkaido will remain basically unchanged in cultivated land protection scenario (CP). In forest protection scenario (FP), the area of forest in Hokkaido will increase by 1580.8 km 2 . It is believed that the ﬁndings reveal that the forest land in Hokkaido has been well protected in the past and will be protected well in the next 20 years. However, in land use planning for future, Hokkaido government and enterprises should pay more attention to the protection of cultivated land.


Introduction
The changing land use/land cover (LULC) exerts impacts on terrestrial environment, climate conditions and biosphere, and more dramatically, human social environment [1]. Land use change is an essential factor for different natural and social elements such as ecological service [2], hydrological characteristic [3], and urban expansion [4,5]. Detecting and analyzing accurate, historical, and continual land use change are essential for land resource management, which is recognized as a measure to better understand and solve socioeconomic and environmental problems [6,7].
Applications of detecting historical LULC change have developed with the advances in Earth observation technologies. Preceding data were interpreted artificially from images taken by satellite and complementary field-based observations [8]. In general, to save time and material cost, remote sensing technologies are commonly applied in monitoring changes in land use type area over time, such as regression analysis models, Markov models, and system dynamics models [31][32][33]. Another is to simulate the spatial distribution of land use, including the land change modeler (LCM) model, The conversion of land use and its effects (CLUE) model, cellular automata (CA) model and future land use simulation (FLUS) model, etc. [34][35][36][37]. Based on the causal relationship between past land use change and related driving force factors, they compute the land use demand and distribution probability, for the aim of simulating future spatiotemporal land use under targeted scenarios [38][39][40].
The FLUS is a model for simulating land use changes and future land use scenarios under the influence of human activities and nature. It has strong temporal and spatial dynamic simulation capabilities [34]. Compared with the traditional CA-based models which need to sample from the changes between the two periods of land use data, this model samples from only the most recent period, which can handle non-linear relationships more effectively by avoiding error transmission. In addition, the FLUS model proposes an adaptive inertial competition mechanism based on roulette selection, which can advantage the demonstration of uncertainty and complexity of the mutual transformation of multiple land use types under the joint influence of natural and human activities.
Hokkaido is the 21st biggest island in the world and second largest in Japan, accounting for 22.4% of the national territory. Unlike other relatively "central" places in Japan which are famous for their manufacturing industry, Hokkaido shows an agriculture industry predominance. The arable land area of Hokkaido accounts for 26.0% of the total arable land in Japan, and its per capita area is 2.17 ha, still ranking at the top in Japan. Besides, the agricultural economic output of Hokkaido was 13.8% of the national total and its calorie-based food self-sufficiency rate is 206%, which means Hokkaido is the veritable food base of Japan [41]. The status of agriculture in Hokkaido will probably affect land use dynamics, and the improvement of agricultural techniques and less population would make this change more complicated.
On the other hand, Hokkaido is significant for its protected areas. There are 11 national parks, 12 prefectural parks, two UNESCO global geo-parks, five national geoparks, and 13 wetlands listed in Ramsar Convention and one World Natural Heritage site in Hokkaido (uu-hokkaido.com (accessed on: 12 March 2021)). These protected areas have considerable ecological and sightseeing value, but simultaneously, might have an impact on future Hokkaido land use patterns. Therefore, it is necessary to figure out the past-to-future land use dynamics of Hokkaido for better understand the land use situation and make optimal governmental decisions about future land use patterns.
Recently less research has been conducted on the entire island of Hokkaido's land use analysis and simulation, possibly due to the large extent of study area and the shortage of data. By using remote sensing tools and a land use simulation model, this study can give fundamental insights into the land use dynamics of Hokkaido. Therefore, the objectives of this study are: (1) to understand the recent spatiotemporal dynamics of land use in Hokkaido by mapping the land use classification in 2000, 2010, 2019, respectively, and (2) to simulate the scenario-based future land use maps in 2040 using an auto-improved FLUS model.

Study Area and Data
Hokkaido is located at the north of Japan (Figure 1), covering the area of about 8.3×10 4 km 2 (excluding controversial territory). Geographically it lies in 40 • 33 N to 45 • 33 N and 139 • 20 E to 148 • 53 E. The elevation ranges from −20 m to 2272 m. Hokkaido is the only one administrative unit called "do" among 47 prefectures in Japan. In subunits, there are 14 subprefectural bureaus, which consist of 64 counties, 35 cities, 129 towns, and 15 villages. The terrain of Hokkaido is high in altitude in its center and relatively low in altitude at all sides near the coast. Mountainous areas account for 60% of the total land, and 40% of this area is occupied by volcanos. On January 1st in 2019, the population of Hokkaido was 5,268,352 and nearly half of is was living in Sapporo, the capital city of Hokkaido. The rest are mostly distributed in other central cities of subprefectural bureaus such as Hakodate, Obihiro and Asahikawa. The age structure of Hokkaido's population is characterized by a serious old-aging trend. People over 65 years old account for nearly 40% of Hokkaido's population.
Meanwhile, the negative natural growth rate of Hokkaido has led to a decreasing total population trend in recent years and future predictions (Hokkaido Population Vision, 2020). The National Institute of Population and Social Security Research (NIPSSR) conducted a population prediction showing the Hokkaido population is expected to drop by nearly 1/5 by 2040 ( Figure 2).
Hokkaido is also one of Japan's main forest lands, and its forests account for about 22% of Japan's total forested area. The national forests of Hokkaido account for more than 50% of the total forest area of Hokkaido, including Mt. Daisetsu, the Hidaka Mountain and other major mountain ranges, creating one of the richest ecosystems in Hokkaido. The national forest of Hokkaido consists of conifers (Sakhalin fir and Sakhalin spruce) and broadleaf trees (Japanese oak, birch, and painted maple), whose views change with the seasons. The Hokkaido Regional Forest Office classifies national forests into four of five types depending on their main functions to conduct appropriate managements: water resource conservation type; landslide prevention type; recreational use type and nature conservation type [42].
Satellite images are downloaded from the United States Geological Survey (USGS) website at a resolution of 30 m, including Landsat Thematic Mapper (TM) data, Landsat Enhanced Thematic Mapper Plus (ETM+) data and Landsat Operational Land Imager (OLI) data. To reduce the missing information caused by snow cover, for each year, images in warm seasons (from 1 May to 31 October) were integrated at each pixel position over time by taking the median. Detailed information is shown in Table 1.  Considering the territory dispute between Japan and Russia regarding the boundaries of Hokkaido (the four northern islands), the administrative shape file in the disputed areas was excluded to fit the subsequent driving force analysis. Therefore, the Hokkaido vector files were not obtained from the Japanese government GSI but rather from the DIVA-GIS website (www.diva-gis.org (accessed on: 12 September 2020)) which removes the above disputed islands from Hokkaido domain.
The change of land use is the result of the combined effect of natural factors and socioeconomic factors, static and dynamic. Therefore, the selected factors should adhere to both natural factors and socio-economic factors to achieve a more accurate simulation of land use change. The selection of driving factors for land use change mainly conforms to the following four principles: (1) accessibility of factor data, which is a key issue for driving factor research and a guarantee for the accurate operation of related models; (2) consistency of factor data: The data should be consistent in time and space; (3) factor quantification, which can be input into the model, and simulation prediction can be realized; and (4) comprehensiveness of factor selection [43]. Based on this, seven land use change factors were selected, including digital elevation model (DEM), slope, distance to road (DTROAD), distance to railway (DTRAIL), organic carbon in topsoil (TOC), and bulk density in topsoil (TBD), population density (PD). In addition, a protected area of national parks from GSI are displayed as a policy refrain on land use change. The detailed descriptions of each factor are illustrated in Table 2: The administration-based masked raster files of land use driving factors were processed at unified resolution at 150 m due to the model calculator ability, as shown in Figure 3.

Methodology
This study includes two main processes: land use classification and future land use simulation. Figure 4 indicates detailed data input, computing methods (algorithm and models) and specific steps in each process.

Land Use Classification
Google Earth Engine (GEE) was used for processing image data. Landsat images are selected at Tiers 1 level, which has been treated systematically, geometrically, and topographically. The study site was mosaicked by clipping the Hokkaido boundary shape file in ArcGIS Pro 2.0 software. Sequentially, this study conducted radiometric and atmospheric correction in GEE and Pixel Information Expert Engine (PIE) (https://engine.piesat.cn/ (accessed on: 26 May 2021)) before classification step. By considering the land use situation of Hokkaido from GSI database, the types of land use were divided into six categories: cultivated land, forest, construction, grassland, waterbody, and others. The descriptions of each land use types are showed in Table 3.
The random forest (RF) classifier was employed in the supervised classification process of the GEE algorithm database to extract the above various land use classes. The RF algorithm can avoid errors such as data noise and over-fitting, which has great advantages for the classification of ground observation image data. Compared with the maximum likelihood method, single decision tree and single layer neural network method, the random forest algorithm is more successful in processing high-dimensional data and has higher accuracy [44].
The result of supervised classification depends heavily on the input training samples. To distinguish features under different conditions, the sample capacity of the initial training data set of the classifier needs to be large enough, especially in com-plex regions [45]. By visual observation on false-color composites (RGB) image to re-strict unchanging area in these three periods, totally 1695 points were selected in training sample set. A precise testing sample set is important for accuracy assessment of classifier and the quality of classification process [46]; hence, this study drew 1071 polygon-based testing samples by referring to the land use maps released by GSI. The number of samples for each land use is showed in Table 4. Table 3. Description of land use types.

Land Use Types Description
Cultivated land Arable land that is intended to plow and sow and raise crops for farming.

Forest
Land dominated by trees, including deciduous forest and evergreen forest Waterbody Rivers, lakes, ponds, wetlands.

Construction
Balconies, terraces (with or without roof), mezzanine floors and other detachable habitable areas.
Grassland Grass, herb, and temporary meadows, including natural and artificial grass.

Others
Land covered by snow for long time and bare land in the top of mountains Five spectral features are applied to characterize typical land use categories, including normalized difference vegetation index (NDVI), normalized difference built-up index (NDBI), modified normalized difference water index (MNDWI), DEM and slope.

Land Use Simulation in 2040 3.2.1. Logistic Regression Analysis for Qualitative Analysis of Driving Factors
Because the FLUS model uses ANN to calculate the distribution probability of each land use type, it has good operability and explanatory properties. However, this method cannot intuitively reflect the positive or negative influences of driving factors on distribution probability of each land use. To reflect the qualitative impact of each driving factor on the distribution of land use, we used logistic regression analysis to calculate the regression constant of each driving factor on land use, thereby analyzing the positive or negative effects of each driving factor on each land type. After reclassing land use in 2019 into 6 single layers for each land use type, this study processed them with seven driving factor raster and its related autocorrelation factor raster into ASCII forms, and then computed results in Statistical Product and Service Solutions (SPSS, Version 25.0, IBM Corp, Armonk, NY, USA).

Autocorrelation Factor
A spatial autocorrelation factor is introduced as the input variable of the artificial neural network (ANN) model to detect the spatial autocorrelation effect to improve the simulation effect of the model. Spatial autocorrelation factors can be expressed as: where, Autocov p is the spatial autocorrelation variable of a land use type on the cellular p; j is a neighborhood cell p; y j indicates the occurrence of a class on the cell j (1 or 0); and w pj is the spatial weight coefficient between cellular p and j. When the distance between p and j is less than threshold distance w pj = 1; otherwise w pj = 0.

Coupled Markov-FLUS Model
This study used a coupled Markov and future land use simulation (FLUS) model to dynamically simulate and predict land use/cover types in Hokkaido. Its construction principle is as follows: The cellular automata (CA) module based on the adaptive inertia mechanism in the FLUS model can only process spatial distribution by handling interaction between cells [34], while the Markov model can be used to predict the quantitative change of land use types. The Markov model cannot compute the changes in the spatial pattern [31].
The Markov model is a stochastic model used to predict the quantitative changes of land use, which dominates in its stability and no aftereffect. No aftereffect indicates that in a process of random development, the future state represented by t + 1 is only related to the current state represented by t. Stability refers to the stability of Markov's process over a long period of time. Dynamic evaluation of land use/cover change features Markov's stability and no aftereffect nature [47]. Therefore, it can be used to predict the changing trend of land use in the study area in the future. Its mathematical expression is shown below: where, S (t+1) represents the state of land use type at t+1 time in the future; S (t) represents the state of the land use type at the current time t; P ij represents the state transition probability matrix of land use type. The FLUS model is a tool for simulating future land use changes under both natural and social elements [48]. It combines an artificial neural network (ANN) algorithm, adaptive inertia, and competition mechanism to improve the traditional CA model, which is suitable for multi-scenario simulation. The back propagation neural network (BPNN) in the FLUS model is a multi-layer feedforward neural network [34]. Its basic composition includes an input layer, one or more hidden layers and an output layer, where the neurons of the input layer and the input the driving factors of land use change correspond to each other, and each neuron in the output layer corresponds to each type of land use. This model further combines the occurrence probability with converting cost, neighborhood effect and a self-adaptive inertia coefficient to estimate the combined probability of each land use grid, which can be relatively input in the panel of FLUS software. Among them, neighborhood effect ranges from 0 to 1, and a larger value indicates a stronger land expansion ability; inertia coefficient can be computed by CA evolution iterations; converting cost indicates the conversion difficulty from the current land use type to the target type, which is another factor shaping the land use dynamics.

Land Use Classification Results
The land use change pattern in Hokkaido for three corresponding years are statistically and spatially illustrated in this section. The pixel-based data is resampled to 150 m resolution for the purpose of being compatible to model settings. The area of each land use and its changing trend are shown in Table 5 and Figure 5, respectively. The spatial land use distributions are disclosed in Figure 6.  As the indication of Figure 6, spatial change of land use in Hokkaido has not been significant during this period. Cultivated land is mainly located in areas dominated by plains, mainly in the center and southeast of Hokkaido. Except for some sporadic distribution the in southern and central parts of Hokkaido, grasslands are mainly distributed in the west of Hokkaido. During this period, part of the grassland and cultivated land expanded to forest, especially in the north of the Tokachi subprefectural Bureau and the west of Nemuro subprefectural Bureau. However, some cultivated land in Hiyama subprefectural Bureau has been converted into grassland and construction. Figure 7 shows some examples of distributed change between 2000 and 2019.    Table 8 present the accuracy assessment results of Hokkaido obtained from RF algorithm with user's accuracy (UA), producer's accuracy (PA), overall accuracy (OA) and Kappa value. Training points for accuracy assessment adopted stratified random sampling method, which depends on the proportional area of each class; therefore, land types with small area easily tend to get the lower accuracy. However, the overall accuracy for all the LULC classification maps exceeded 88%.    Table 9 reflects the regression results of regression constants, regression coefficients between land use distribution and driving factors, and ROC value for each group. From the regression results, soil quality TOC and TBD have a positive impact on the distribution of cultivated land, while DTROAD, DTRAIL and slope influence the distribution negatively. This is consistent with a previous study [49]. Forests can be more dramatically affected by slope. DTROAD and DTRAIL shows as negative drivers targeted to upbringing of waterbody and construction. However, it is showed that distribution of forest and DTRAIL have positive relationship, which may be due to the special topography trait of Hokkaido. In general, it is believed that the ROC value greater than 0.7 can indicate a reliable simulation result; therefore, this study can better explain the effects of land use driving factors with all ROCs over 79%.

Accuracy Assessment for Land Use Simulation
By importing the current land use data in 2010 in the CA module of the FLUS software, based on the occurrence probability map computed by inputting driving factors to ANN with and without autocorrelation variables, two 2019 land use simulation maps are generated for accuracy comparison. The setting parameters include the number of iterations 300; the size of the neighborhood 3; the model acceleration factor 1, and the simulated land quantity to the actual land use quantity in 2019. The optimal land-type transformation rules are shown as Table 10, where 0 means that it cannot be changed, and 1 means that it can be changed. This study mainly compared the overall accuracy (OA) under two different conditions, we found after inducing autocorrelation factor, the OA of simulation climb slightly from 0.831 to 0.834 and Kappa value increased from 0.665 to 0.671, therefore, we will apply this factor into simulation process in subsequent steps.

Scenario-Based Future Land Use Prediction
In natural development (ND) scenario, land type transformation rules in Table 10 have been applied in process. The national park is added as a restricted area, and the natural evolution forecast area based on Markov estimation in Table 11 is used as the land use demand area for 2040. Considering that the area of cultivated land under the ND scenario has been reduced to a certain extent, this study constructed a CP scenario by restraining this change. CP scenario is based on ND scenario, while it defines area of cultivated land in 2019 as a protected area which is added into restricted area. Similarly, FP scenario is based on ND scenario, while it defines forest area in 2019 as a protected area which is added into restricted area. The land use distributions under the three scenarios are shown in Figure 8.   The results of this study show that Hokkaido has had a significant decline in cultivated land in the past and will the decline will continue in the next 20 years, while the area of grassland will have a clear upward trend. This shows that during the period, the abandonment of cultivated land will become more serious. According to the simulation analysis based on the driving factor, the abandoned farmland will be converted into grassland (in this case, abandoned area) to a large extent, and part of it will be occupied by construction, which is consistent with the trend reflected in both Tables 6 and 7 and Figure 9. The abandonment of farmland in Hokkaido has become an issue affecting the sustainability of land use since the last 1980s [50], which is a direct reflection of the future reduction of cultivated land in Hokkaido. Previous studies have shown that the main factors for the abandonment of farmland in Hokkaido include natural and social factors: natural factors consist of geographic location, slope and elevation, soil quality, while social factors mainly include policy and demographic changes. It is worth mentioning that the impact of policies on farmland abandonment can be two contrast ways. For instance, the "Hokkaido Comprehensive Development Plan (HCDP)" aims to promote the relocation, concentration, and intensification of farmland, and to increase agricultural productivity and farmers' income [49], which encourages to abandoned barren and hard-to-use farmland. However, the "Direct Payment (DP) system" in Hokkaido slowed down the growth rate of abandoned farmland through financial subsidies to farmers [51]. Therefore, to maximize the sustainability of cultivated land use in Hokkaido, this study did not simulate the impact of certain policies. Instead, it directs the effects that policy can make, that is, to minimize the abandonment of cultivated land to construct CP scenario.
The demographic changes in Hokkaido will also have an impact on the abandonment of farmland. The results of logistic regression analysis in this study show that population density and the distribution of cultivated land are negatively correlated, which means that the more concentrated the population is in urban areas, the lower the probability of cultivated land. According to the government's employment statistics (www.pref.hokkaido.lg.jp (accessed on: 15 December 2020)) from 2000 to 2019, the proportion of people engaged in farming in Hokkaido has decreased and the number of farmers over 65 years old has consistently risen, while the urbanized population is increasing, and it can be claimed that the increase of urban population will promote the expansion of the urban construction [52][53][54]. This is consistent with the results of the decrease of cultivated land in Figure 8a and the increase of construction in Figure 8e.
Compared with the area changes of other land use types, waterbody and others did not change significantly from 2000 to 2019 and are not estimated to change much in the future. For waterbody, the changes in this stage are likely to be caused by the classification error of the cultivated land, especially in areas with paddy fields, which occurs frequently in the Ishikari Plain. As far as others, the main reason for its low accuracy is the uncertainty of snow distribution on the surface. Except for some high-altitude mountains, most of the snow-covered areas have a great relationship with the climatic conditions of the target year.

Comparison of Scenario-Based Land Use Situations
As the area of cultivated land in Hokkaido has a clear downward trend in the context of natural development (ND), This study achieves the sustainable use of cultivated land by restricting the transfer of cultivated land to other land ( Table 12). The area of cultivated land obtained in ND scenario is compared with the CP scenario and the FP scenario, respectively. An increase of about 25% in cultivated land in CP scenario indicates that the policy of protecting cultivated land can effectively improve the sustainable use of cultivated land. In addition, the amount and spatial distribution of cultivated land use in 2019 are almost the same as in the CP scenario, which shows that the demand for cultivated land in Hokkaido has been basically saturated, and it may only remain unchanged or decrease in the future. By comparing the ND scenario and the CP scenario, we can summarize the spatial change trends of the cultivated land in the two scenarios in the future, mainly showing the following three points. Firstly, in the ND scenario, the cultivated land is more likely to be converted into forest, grassland, or construction due to the concentration of agriculture. This includes the expansion of forest land from the center to the surrounding area and the increase of scattered buildings (Table 13 and   Due to the important position of forest land on Hokkaido, this study explored the possible changes of forest land in the future by constructing FP protection scenarios. In CP scenario, a reduction of 1310 km 2 was found in the forest compared to the ND scenario ( Figure 11), while in FP scenario, the forest has increased by 509.7 km 2 compared to the ND scenario ( Figure 12). The increase in forest land is mainly due to the conversion from abandoned farmland and the slowing down of grassland expansion, especially in mass graze-based grassland in east of Hokkaido (Figure 13a,c and Figure 13b,d).

Recommendations and Suggestions
Based on the comparative analysis of the characteristics of land use evolution, existing problems, and multi-scenario simulation results in this study, the following suggestions are put forward for the optimization and adjustment of future land use in Hokkaido: (1) Our results show that restricting land use areas can effectively enhance the sustainable use of cultivated land. Therefore, it is advisable to strengthen the protection of cultivated land, by regulating a basic cultivated land protection system and compensation system (such as Direct Payment) should be established to ensure the balance of cultivated land occupation and its quality. (2) Promote the transformation of government functions by strengthening land regulation and supervision, meanwhile restrict the trading of private cultivated land in Hokkaido. Public funds should be used to restore the use of abandoned arable land, and at the same time, through agricultural technology, to ensure that Hokkaido's food production continues to meet the minimum standards of national food security. (3) The results show a mass area of grassland replacing cultivated land in the future. Hence, by centralizing grazing activities, the expansion of Hokkaido's grassland into cultivated land and forest can be effectively restrained, especially at the boundaries between grassland and these lands. Simultaneously, the unused grassland converted from abandoned farmland should be conserved to ensure that its fertility is capable of being converted back into cultivated land in the future. (4) From the regression results, this study concludes intensive population in urban area is not beneficial to the development of cultivated land. Consequently, controlling the concentration of population in Hokkaido cities and encouraging people to relocate to the non-central areas is necessary for cultivated land protection. By improving transportation, medical care, education and other functions, the population distribution of Hokkaido can be dispersed.
The fifth basic land use plan (www.pref.hokkaido.lg.jp (accessed on: 22 May 2021)) of Hokkaido stipulates that the cultivated land in the basic farming area cannot be used for other purposes, which fits well with recommendations 1 and 2. The existing farmland will be protected to the maximum extent, and the fertile but abandoned cultivated land will be protected and appropriately improved. For cultivated land in non-agricultural areas, where adjustments are expected to conduct non-agricultural land use planning (such as urban planning), stakeholders will follow the adjusted planning older which is used to define the "Good Cultivated Land" or not, and this "Good Cultivated Land" will be properly preserved. However, the criteria for formulating this order are multifaceted, not only related to the land itself, but also to its geographic location and the attitude of the owner, which complicates the protection of cultivated land outside of agricultural areas. Policymakers are supposed to clearly define what is "Good Cultivated Land" but not depend on an obscure consultation with landlord.

Limitations and Future Works
In general, some shortcomings of the study can be concluded based on the following four points: (1) The remote sensing data for making land use classification are at a resolution of 30 m and the simulation process is at a resolution of 150 m, which indicates there are still some details requiring more accurate interpretation. (2) The classifications of cultivated and waterbody, cultivated land and grassland without considering the phenological law, are more likely to be confused. Moreover, due to the large area of Hokkaido, it is difficult to collect training sample sets based on field surveys, resulting in certain classification errors. (3) Driving factors have relatively limited impact on land use, and different types of driving factor weights cannot be assigned, that is, the influence weights of factors on land use cannot be grasped well. (4) Moreover, scenarios differ from absolute restricted area, but the actual situation is hard to develop towards these directions.
Accordingly, further studies are supposed to include the below solutions: (1) use of relatively high resolution satellite data to optimize land use classification process if it is technically and economically feasible; (2) to investigate the phenological laws of various vegetation types, especially paddy rice in cultivated land, for improving the classification accuracy; (3) to use appropriate weight analysis methods to rank the influence of land use driving factors; (4) to collect policy documents to support the establishment of target future land use scenarios.

Conclusions
This research used remote sensing classification methods based on a random forest algorithm to monitor the spatiotemporal changes in six land classes (cultivated land, forest, waterbody, construction, grassland, and others) on Hokkaido from 2000 to 2019, with an overall accuracy of classification over 88%. Based on the analysis of land use changes and driving factors including DEM, slope, PD, DTRAIL, DTROAD, TOC, TBD and autocorrelation factors added to improve simulation accuracy, the coupled Markov-FLUS model is applied by changing the restricted area to simulate spatial distributions of land use in natural development scenario, cultivated land protection scenario, and forest protection situation of Hokkaido in 2040. The land use classification results show from the period of 2000 to 2019, the area of cultivated land in Hokkaido shows a downward trend with a reduction of 642.1 km 2 , while construction has increased by 111.0 km 2 and grassland expanded by 2179.4 km 2 , while the forest and waterbody areas have not changed significantly. The simulation results of future land use changes show that under the natural development scenario, the cultivated land in Hokkaido will be reduced by about 25% by 2040, and the reduced cultivated land is more likely to be converted into grassland and forest land. In the cultivated land protection scenario, the area of cultivated land remains basically unchanged, which shows that the use of cultivated land in Hokkaido has reached saturation, and about 239.0 km 2 forest land will be invaded by other land types. In forest protection scenario, the area of forest in Hokkaido will increased by around 1580.8 km 2 , but the change trend of cultivated land is following natural development scenario. Grassland, especially in mass intensive graze-based area of eastern Hokkaido, will have an obvious expansion trend. This study provides a reference for future land use planning and management in Hokkaido, and it is expected to support the sustainable land use development with different targets.

Data Availability Statement:
The data presented in this study are cited within the article.