Open AccessArticle Spatiotemporal Trajectories and Divergent Drivers of Cropland Non-Grain Use: Evidence from the Changsha–Zhuzhou–Xiangtan Urban Agglomeration, China 1 College of Landscape Architecture and Art Design, Hunan Agricultural University, Changsha 410128, China 2 Hunan Institute of Agricultural Economics and Information, Hunan Academy of Agricultural Sciences, Changsha 410125, China * Author to whom correspondence should be addressed. † These authors contributed equally to this work. Land 2026, 15(6), 985; https://doi.org/10.3390/land15060985 (registering DOI) Submission received: 1 May 2026 / Revised: 1 June 2026 / Accepted: 2 June 2026 / Published: 4 June 2026 Cropland non-grain use has become an important challenge for food security and cropland governance in rapidly urbanising agricultural regions, yet its trajectory heterogeneity and the divergence between current spatial patterns and long-term-change mechanisms remain insufficiently understood. Taking the Changsha–Zhuzhou–Xiangtan (CZT) urban agglomeration in China as a case, this study quantified the cropland non-grain rate (NGR) on a 1 km grid for 2000, 2010, and 2020, classified grid-level transition trajectories, and developed three temporally structured eXtreme Gradient Boosting (XGBoost) models with spatial block cross-validation, Shapley additive explanations (SHAP) interpretation, and geographically explicit SHAP (GeoSHAP) local attribution. The results show that low-NGR and stable grids formed the dominant regional background, while recent NGR increases were mainly concentrated along the urban development corridor and metropolitan fringe. Current NGR status and long-term NGR change showed divergent explanatory structures. The current spatial pattern was mainly associated with terrain constraints and contemporary urban pressure, whereas long-term change was more strongly conditioned by baseline urbanisation and subsequent urban–environmental changes. Nonlinear dependence analysis further identified model-derived response zones related to slope, impervious surface conditions, hydrothermal change, and hydrological proximity. GeoSHAP mapping revealed that locally dominant mechanisms varied substantially across the study area, indicating that cropland non-grain use was shaped by spatially heterogeneous combinations of terrain, urbanisation, hydrothermal background, and hydrological context. These findings support a shift from aggregate status monitoring toward trajectory-specific and mechanism-differentiated cropland management in urban agglomerations. 1. Introduction This study focuses on the Changsha–Zhuzhou–Xiangtan (CZT) urban agglomeration in central China, a region that combines an important grain-production base in Hunan Province—one of China’s major rice-producing regions—with rapid urban corridor expansion and a diverse terrain backdrop of plains, hilly land, and low mountains. Rather than characterising non-grain use through a single cross-sectional model, we adopt a three-model framework that jointly analyses: (1) the 2020 NGR spatial pattern (Model A); (2) long-term ΔNGR over 2000–2020 using both baseline urbanisation and change predictors (Model B); and (3) ΔNGR using only change and static predictors (Model C). Trajectory classification across 2000–2010, 2010–2020, and 2000–2020 further characterises the process heterogeneity underlying cross-sectional patterns. This design allows the explanatory structures of current NGR status and long-term NGR changes to be explicitly compared, and makes it possible to evaluate whether historical urbanisation context provides additional explanatory information beyond concurrent change variables. SHAP-based global driver analysis, nonlinear dependence plots, interaction networks, and dominant local driver maps are then used to characterise mechanism heterogeneity across the three models. Specifically, this study addresses three questions: (1) What spatiotemporal patterns and trajectory types did NGR exhibit in the CZT urban agglomeration from 2000 to 2020? (2) How did the drivers of the 2020 NGR pattern differ from those of long-term ΔNGR, and what explanatory role did baseline urbanisation play? (3) What nonlinear responses, interaction structures, and locally dominant mechanisms characterised NGR drivers, and how can they inform spatially differentiated cropland management? By addressing these questions, this study shows that cropland non-grain use in the CZT urban agglomeration was dominated by a stable low-NGR regional background, while meaningful increases were concentrated mainly along urban fringe and development corridor areas. The modelling results further demonstrate that current NGR status and long-term NGR change were associated with different explanatory structures: the former was mainly related to terrain constraints and contemporary urban pressure, whereas the latter was more strongly conditioned by baseline urbanisation and subsequent urban–environmental changes. In addition, SHAP and geographically explicit SHAP (GeoSHAP) analyses revealed nonlinear response zones, coupled interaction structures, and spatially heterogeneous local dominant drivers, supporting a shift from aggregate status monitoring toward trajectory-specific and mechanism-differentiated cropland management. 2. Materials and Methods 2.1. Analytical Framework The analytical framework was designed to distinguish the current spatial pattern of cropland non-grain use from its long-term-change process and to identify locally dominant explanatory mechanisms ( Figure 1). The workflow consisted of four steps. First, multi-source crop and land cover, terrain, soil, hydrological, urbanisation, accessibility, and climate/hydrothermal datasets were harmonised to a common 1 km grid in EPSG:32649. Second, the major-grain layer (MGL) and the China Land Cover Dataset (CLCD) cropland mask were overlaid to calculate the cropland non-grain rate (NGR) for 2000, 2010, and 2020. Valid grids were then screened using a cropland-share threshold of 0.05, and grid-level trajectories were classified for 2000–2010, 2010–2020, and 2000–2020. Third, three temporally structured XGBoost models were developed after variance inflation factor (VIF) screening and evaluated using spatial block cross-validation. Model A explained the 2020 NGR pattern, Model B explained long-term ΔNGR using both baseline urbanisation and change variables, and Model C explained the same ΔNGR target after excluding baseline urbanisation variables. Finally, SHAP and GeoSHAP were used to interpret the model outputs, including global driver importance, nonlinear dependence, interaction networks, and dominant local drivers, thereby supporting trajectory- and mechanism-informed cropland management. 2.2. Study Area The CZT urban agglomeration is located in the central–eastern part of Hunan Province, China (111°53′–114°15′ E, 26°03′–28°41′ N). It comprises the three prefecture-level cities of Changsha, Zhuzhou, and Xiangtan and covers a total land area of approximately 28,000 km 2 ( Figure 2). The region has a subtropical humid monsoon climate, abundant precipitation, and a pronounced coincidence of heat and moisture. Its topography is dominated by an interwoven pattern of low mountains, hills, and river-valley plains. These favourable natural conditions have long supported the region as one of China’s major double-cropping rice production areas, providing a strong basis for grain production. The main stem of the Xiang River runs from south to north across the region, shaping contiguous agricultural hinterlands in the river valleys while also serving as a major axis for population movement and industrial concentration. As an important growth pole in the middle reaches of the Yangtze River, the CZT urban agglomeration is also facing intense land-use competition under continued urbanisation and spatial restructuring. High-quality cropland concentrated in river-valley plains overlaps substantially with corridors of active urban expansion, making the tension between cropland protection and construction development particularly pronounced. Beyond the direct encroachment associated with expanding built-up land, operational constraints imposed by hilly terrain and the urban–rural land-rent gradient under land marketisation further intensify the spatial differentiation of cropland use. The region therefore combines a strong agricultural base, substantial urbanisation pressure, and a complex topographic background, making it a representative case for examining the spatial evolution, nonlinear response zones, and locally dominant mechanisms of cropland non-grain use and NGR change. 2.3. Data Collection and Harmonisation To integrate datasets with different native spatial resolutions, all raster and vector datasets were projected or transformed to WGS 1984/UTM Zone 49N (EPSG:32649) [ 36] and summarised within a regular 1 km grid. Area-proportion variables, including cropland share, major-grain share, NGR, and impervious surface proportion, were calculated as class-area proportions within each grid. Continuous raster variables, including elevation, slope, soil organic carbon, climate variables, nighttime-light intensity, and population density, were aggregated using grid-level mean values. Long-term-change variables were calculated by subtracting baseline values from terminal-year values, while hydrothermal trend variables were derived from TerraClimate time-series trends during 2000–2020. Road density was calculated as the total length of major roads within each grid divided by grid area, and distance variables were measured as Euclidean distances from grid centroids to the nearest town centre, highway, river, or water body. Elevation and slope were derived from the SRTM DEM after projection, clipping, and NoData masking, and were aggregated using the same 1 km grid framework. 2.4. Measurement of NGR and NGR Change Trajectories NGR was calculated for 2000, 2010, and 2020 by overlaying the major-grain layer (MGL) with the CLCD cropland mask within each 1 km grid [ 27, 28]. These three observation years were selected because both MGL and CLCD provide comparable decadal data for these years, allowing a consistent grid-level NGR calculation across the study period. In addition, the three time points broadly represent the initial, intermediate, and recent stages of cropland-use transition in the CZT urban agglomeration during rapid urbanisation and subsequent strengthening of cropland protection. The two sub-periods, 2000–2010 and 2010–2020, were used to distinguish earlier and later transition processes, whereas the full 2000–2020 period was used to identify long-term trajectory types. For each grid cell, the total cropland area and the area occupied by major-grain crops were calculated using the CLCD cropland mask and the MGL dataset, respectively. NGR was defined as the proportion of cropland not occupied by major-grain crops: NGR i = NGR i , t = 1 − A i , t MGL A i , t CL (1) where A i , t MGL is the area of major-grain crops in grid i and year t , and A i , t CL is the total cropland area in the same grid and year. NGR i , t was stored as a proportion ranging from 0 to 1 in model fitting and was converted to percentages for descriptive statistics and figure interpretation. A value close to 0 indicates that cropland in the grid was mainly used for major-grain production, whereas a value close to 1 indicates a high proportion of non-grain use within the retained cropland. To avoid unstable NGR estimates caused by very small cropland areas, a cropland-share threshold was applied before trajectory classification and model fitting. Grids with cropland share greater than 0.05 were retained as valid cropland grids, corresponding to more than 5% of the 1 km grid area. For Model A, which explained the 2020 NGR pattern, the threshold was applied to cropland share in 2020. For trajectory classification and the ∆ NGR models, grids were retained only when cropland share exceeded 0.05 in both 2000 and 2020 and NGR values were valid for 2000, 2010, and 2020. Grids with zero or insufficient cropland area, undefined NGR values, NoData values, abnormal values, or missing predictor values were treated as invalid and excluded from the corresponding analysis. Percentages reported for NGR classes and trajectory types were calculated using valid grids as the denominator. NGR changes were calculated for 2000–2010, 2010–2020, and 2000–2020 as follows: ∆ NGR i , t 1 − t 2 = NGR i , t 2 − NGR i , t 1 (2) where ∆ NGR i , t 1 − t 2 denotes the change in NGR for grid i from the initial year t 1 to the terminal year t 2 . Because NGR was represented as a proportion, a threshold of 0.05, equivalent to 5 percentage points, was used to distinguish substantive changes from small fluctuations. For each two-period comparison, grids were classified as NGR-increasing when ∆ NGR > 0.05 , NGR-decreasing when ∆ NGR 5 were removed stepwise until all retained predictors satisfied the VIF criterion. Detailed VIF screening results are provided in Tables S1–S3. After VIF screening, Model A retained 14 predictors. PDSI2020 was removed because of high multicollinearity with other 2020 climate/hydrothermal variables. Model B retained all 13 candidate predictors, and Model C retained all 10 candidate predictors. The final model-specific driver systems are shown in Table 2. 2.7. XGBoost Modelling and Spatial Block Cross-Validation XGBoost was used to model the relationships between NGR or ∆ N G R and the selected drivers [ 20]. Model A used ΔNGR2020 as the response variable, whereas Models B and C used ΔNGR2000–2020. For a dataset containing n grid samples and m predictors, XGBoost predicts the response by additively combining K regression trees: y ^ i = ∑ k = 1 K f k ( x i ) , f k ∈ F (4) where y ^ i is the predicted response for grid i , x i is the vector of explanatory variables for the same grid, f k is the k -th regression tree, and F denotes the space of regression trees. The objective function combines the prediction error with regularisation terms that penalise model complexity: L = ∑ i = 1 n l ( y i , y ^ i ) + ∑ k = 1 K Ω ( f k ) (5) where y i is the observed response, l is the loss function, and Ω ( f k ) is the regularisation term for tree f k . This formulation allows the model to capture nonlinear relationships and variable interactions while controlling tree complexity. To reduce spatial leakage, model performance was evaluated using spatial block cross-validation rather than random cross-validation [ 23, 24]. Regular 10 km spatial blocks were generated in the projected coordinate system, and samples belonging to the same block were assigned to the same fold. A five-fold grouped cross-validation procedure was then used to obtain out-of-fold (OOF) predictions. Model performance was evaluated using coefficient of determination (R 2), root mean square error (RMSE), and mean absolute error (MAE) calculated from OOF predictions. Hyperparameters were optimised using Optuna with a Tree-structured Parzen Estimator search strategy and block-aware inner cross-validation [ 38]. A final model was then trained on the full modelling sample using the selected hyperparameters for SHAP-based interpretation. 2.8. SHAP Interpretation, Interaction Network, and Geographically Explicit SHAP (GeoSHAP) Mapping To interpret the trained XGBoost models, SHAP values were calculated for each predictor and each grid [ 21, 22]. SHAP decomposes model predictions into feature-specific additive contributions, allowing the direction and relative magnitude of each predictor’s contribution to be examined at the grid level. In this study, SHAP was used to identify global driver importance, nonlinear dependence patterns, approximate transition zones, pairwise interaction structures, and locally dominant drivers. Because SHAP values are derived from fitted machine learning models, they were interpreted as model-based explanatory signals rather than causal effects. For feature j in grid i , the SHAP value can be expressed as ϕ i , j = ∑ S ω ( S ) [ f i ( S , j ) − f i ( S ) ] , ω ( S ) = | S | ! ( | F | − | S | − 1 ) ! | F | ! (6) where F is the full feature set, S is a possible subset of features excluding feature j , ω ( S ) is the Shapley weight, and f i ( S ) is the expected model output for grid i when only the features in subset S are used. A positive SHAP value indicates that the feature contributes to a higher predicted response relative to the model baseline, whereas a negative SHAP value indicates that the feature contributes to a lower predicted response. In both cases, the sign and magnitude of the SHAP value are conditional on the fitted model, the other predictors, and the analysed sample. Therefore, SHAP values should be interpreted as model-based contributions rather than causal effects: a positive value does not imply that a driver universally increases NGR or ΔNGR, and a negative value does not imply that a driver universally reduces NGR or ΔNGR. For Model A, the response was NGR2020; for Models B and C, the response was ΔNGR2000–2020. Global feature importance was calculated using the mean absolute SHAP value: I j = 1 n ∑ i = 1 n | ϕ i , j | (7) SHAP dependence plots were used to examine nonlinear response patterns of dominant predictors. The smoothed dependence curves were used to identify approximate transition zones where SHAP contributions changed direction or intensity. These transition zones were treated as empirical, model-derived response ranges within the observed data rather than fixed regulatory thresholds or universal causal breakpoints. Pairwise interaction effects were summarised using SHAP interaction values: I j , k = 1 n ∑ i = 1 n | ϕ i , j , k | (8) where ϕ i , j , k denotes the interaction contribution between predictors j k for grid i . The resulting interaction network used node size and colour to represent individual predictor importance and edge width and colour to represent pairwise interaction intensity. The network was therefore used to compare whether the explanatory structure of each model was dominated by isolated predictors or by coupled driver relationships. Finally, SHAP values were spatially linked back to each 1 km grid to construct a GeoSHAP interpretation layer. For each grid, the locally dominant driver was defined as the predictor with the largest absolute SHAP contribution: D i = a r g m a x j | ϕ i , j | (9) where D i denotes the locally dominant driver in grid i . The dominant-driver map identifies the predictor with the largest model-derived contribution to each grid-level prediction. Because this definition is based on relative SHAP magnitude, it indicates local model dominance rather than exclusive physical control over NGR or ∆ N G R . 3. Results 3.1. Spatiotemporal Patterns of Cropland Non-Grain Use and NGR Trajectories From 2000 to 2020, NGR in the CZT urban agglomeration showed an overall spatial pattern characterised by widely distributed low values and locally clustered medium- to high-value patches ( Figure 3), while the trajectory classification further indicated that stable trajectories dominated the region and NGR-increasing trajectories were spatially concentrated in urban-fringe and corridor areas ( Figure 4). Across the three cross-sectional years, low-NGR grids remained dominant. Grids with NGR 30% accounted for 12.40% in 2000, declined to 9.54% in 2010, and then increased to 11.25% in 2020 ( Figure 5a). This temporal fluctuation suggests that high-level non-grain use did not expand uniformly across the whole region, but evolved through localised patch contraction and re-expansion. Spatially, these higher-value grids were more evident around the Changsha metropolitan fringe and extended southward along the core urban corridor. Compared with the broad agricultural hinterlands dominated by low NGR, the corridor and near-urban fringe areas showed more heterogeneous and unstable cropland-use patterns. The trajectory classification further revealed the process characteristics behind these cross-sectional patterns ( Figure 4). During 2000–2010, stable grids dominated, accounting for 77.08% of valid grids, while NGR-decreasing grids accounted for 17.02% and NGR-increasing grids accounted for 5.90% ( Figure 5b). This indicates that the first decade was mainly characterised by stability and partial decline in NGR, with non-grain increase concentrated in limited patches around the urban fringe. During 2010–2020, the proportion of stable grids remained high at 77.97%, but the share of NGR-increasing grids rose to 10.46%, while NGR-decreasing grids declined to 11.57%. This shift suggests that the second decade experienced a stronger tendency toward localised non-grain increase, especially around the central Changsha area and along the CZT corridor. Over the full 2000–2020 period, stable grids still formed the dominant trajectory type, accounting for 70.63% of valid grids ( Figure 5b). However, long-term transition patterns became more differentiated: NGR-decreasing grids accounted for 16.31%, NGR-increasing grids accounted for 7.40%, and reversal grids accounted for 5.66%. Spatially, NGR-increasing trajectories were concentrated mainly in the outward-expanding urban core belt and along the major development corridor, indicating that newly intensified non-grain use was closely associated with urban expansion zones. NGR-decreasing trajectories were more widely scattered in peripheral agricultural areas, implying partial recovery or weakening of non-grain use in some cropland grids. Reversal trajectories were distributed in fragmented patches, reflecting bidirectional fluctuations rather than a simple monotonic transition process. 3.2. Predictive Performance of the Three XGBoost Models The predictive performance of the three XGBoost models differed across modelling targets ( Table 3; Figure S1). Model A, which explained the 2020 cross-sectional NGR pattern, achieved the highest predictive performance, with an in-sample R 2 of 0.466 and a spatial block OOF R 2 of 0.433. Its OOF RMSE and MAE were 0.057 and 0.040, equivalent to 5.7 and 4.0 percentage points, respectively. The relatively small difference between in-sample and OOF performance indicates that Model A captured the spatial variation in NGR2020 with moderate predictive stability under spatial block cross-validation. The two long-term ΔNGR models showed lower R 2 values than the cross-sectional model, reflecting the greater difficulty of explaining grid-level NGR changes over two decades. Model B, which incorporated both baseline urbanisation and change variables, achieved an in-sample R 2 of 0.303 and an OOF R 2 of 0.259, with OOF RMSE and MAE values of 0.042 and 0.028, respectively. Model C, which excluded baseline urbanisation variables, showed a lower in-sample R 2 of 0.246 and an OOF R 2 of 0.208, with OOF RMSE and MAE values of 0.044 and 0.029, respectively. Although the explanatory power of the two ΔNGR models was lower, their prediction errors remained small in absolute terms because most grid-level NGR changes were limited in magnitude. Model B consistently outperformed Model C across R 2, RMSE, and MAE, indicating that baseline urbanisation variables provided additional explanatory information for long-term NGR change. Accordingly, Model A was used to interpret the current spatial pattern of NGR, Model B served as the main model for explaining long-term ΔNGR mechanisms, and Model C was retained as a stricter change-based comparison model for evaluating the contribution of baseline urbanisation context. 3.3. Global Drivers of NGR Patterns and Changes The global SHAP results indicate that the dominant drivers differed substantially among the three models ( Figure 6). In Model A, which explained the spatial pattern of NGR in 2020 ( Figure 6a), slope was the most important predictor, accounting for 29.34% of the total mean absolute SHAP contribution. It was followed by Imp2020 (23.66%), Tmean2020 (9.69%), Pre2020 (9.50%), and SOC (4.08%). At the category level, terrain/soil factors contributed the largest share (37.02%), followed by 2020 urban pressure (30.99%) and 2020 climate/hydrothermal conditions (22.28%). This suggests that the 2020 NGR pattern was mainly associated with terrain constraints, contemporary urban development pressure, and climatic background conditions. In Model B, which explained ΔNGR during 2000–2020 using both baseline and change variables ( Figure 6b), the most important predictor was Imp2000, contributing 33.39% of the total SHAP importance. The next most important variables were SRad_TS (18.74%), dImp (11.87%), slope (9.41%), and Dist_Riv (8.02%). At the category level, baseline urbanisation had the largest contribution (38.32%), followed by climate/hydrothermal change (22.70%) and urban change (18.25%). Compared with Model A, this result shows a clear shift from present-state conditions to historical background and long-term process variables, indicating that NGR changes were not only related to recent urban expansion but also conditioned by the initial urbanisation context. In Model C, which excluded baseline urbanisation variables and used only static and change variables to explain ΔNGR ( Figure 6c), dImp became the most important predictor, contributing 22.47%. It was followed by slope (21.86%), SRad_TS (18.99%), Dist_Riv (13.85%), and dNTL (7.21%). At the category level, urban change contributed the largest share (31.79%), followed by terrain/soil factors (28.74%) and climate/hydrothermal change (24.63%). The persistence of dImp and dNTL among the leading predictors in Model C suggests that urban expansion and nighttime-light intensification remained important explanatory signals even after removing baseline urbanisation variables. Overall, the three SHAP panels show a clear difference between the drivers of current NGR patterns and those of long-term NGR changes. The 2020 NGR pattern was mainly associated with terrain/soil background and contemporary urban pressure, whereas long-term ΔNGR was more strongly linked to baseline urbanisation, urban expansion, and climate/hydrothermal changes. This distinction supports the separation of cross-sectional pattern analysis and transition-process modelling in the subsequent interpretation. 3.4. Nonlinear Dependence of Dominant NGR Drivers The dependence plots further reveal the nonlinear effects of the dominant predictors identified by the global SHAP analysis ( Figure 7). For each model, the four most important variables were selected because they accounted for most of the total individual SHAP contribution. Specifically, the top four predictors contributed 72.19% in Model A, 73.41% in Model B, and 77.17% in Model C. Therefore, the following interpretation focuses on these dominant variables to clarify their model-derived response transitions and nonlinear effects. For Model A, which explained the spatial pattern of NGR in 2020, the four dominant variables were slope, Imp2020, Tmean2020, and Pre2020 ( Figure 7a–d). Slope showed a clear positive transition around 5.88° ( Figure 7a). Below this approximate transition point, its SHAP values were mostly negative, whereas above it the contribution became positive and increased rapidly, indicating that higher terrain constraints were associated with higher NGR2020. Imp2020 showed a sharply nonlinear response ( Figure 7b). Its SHAP value was positive at very low impervious surface proportions but shifted rapidly to negative after approximately 0.01. Because NGR was calculated only within retained cropland grids, this negative segment should be interpreted as a conditional model response within the cropland mask. It indicates that, among valid cropland grids, higher Imp2020 values after this point were associated with lower predicted NGR2020 relative to the model baseline, rather than showing that impervious surfaces generally reduce non-grain use. Tmean2020 showed an approximate transition near 18.12 °C ( Figure 7c); values below this point generally contributed positively, while higher values produced negative SHAP contributions. By contrast, Pre2020 exhibited a positive transition around 1469.30 mm ( Figure 7d), after which its contribution to NGR2020 became positive and gradually increased. For Model B, which represented the baseline-change explanation of ΔNGR during 2000–2020, the top four predictors were Imp2000, SRad_TS, dImp, and slope ( Figure 7e–h). Imp2000 showed a positive transition at approximately 0.01 ( Figure 7e). When baseline impervious surface was very low, its contribution was negative, but it became positive as Imp2000 increased, implying that initial urbanisation background was a key condition associated with subsequent NGR change. SRad_TS also showed a distinct transition near 4.38 ( Figure 7f), with negative SHAP values below this point and positive contributions above it. The effect of dImp was more complex ( Figure 7g). Although the curve shifted rapidly after very small positive changes in impervious surface, most SHAP values remained negative at the global level, suggesting that the association between impervious surface expansion and ΔNGR varied across different transition pathways rather than acting uniformly across all grids. Slope displayed a non-monotonic pattern ( Figure 7h), with positive contributions at low values, negative effects over the middle range, and a gradual recovery toward zero or positive values after approximately 9.91°. Overall, Figure 7 shows that the dominant drivers of NGR patterns and changes were characterised by nonlinear and model-dependent responses rather than simple linear effects. The cross-sectional NGR pattern was mainly associated with response transitions in slope, impervious surface, temperature, and precipitation, whereas long-term NGR changes were more sensitive to baseline impervious surface, impervious surface change, solar-radiation change, and hydrological proximity. These transition points should be interpreted as empirical response zones derived from the fitted XGBoost–SHAP models, rather than as deterministic regulatory thresholds. They provide a basis for interpreting pathway-specific mechanisms and interaction effects, but their transferability requires further validation across different regions and temporal settings. 3.5. Interaction Network of NGR Drivers The interaction networks showed clear differences in coupled driver structures among the three XGBoost models ( Figure 8). In Model A, the strongest interaction occurred between Slp and Imp2020, which were also the two leading predictors in the global SHAP ranking ( Figure 8a). This indicates that the 2020 NGR pattern was associated not only with terrain constraints or contemporary impervious surface conditions separately, but also with their joint variation. The relatively strong links among Imp2020, Tmean2020, Pre2020, and Slp further suggest that current NGR was shaped by a combined terrain–urban–hydrothermal structure. In Model B, the interaction network shifted toward the baseline urbanisation context ( Figure 8b). The strongest link was Imp2000–SRad_TS, and Imp2000 was also connected with dImp, Slp, and Dist_Riv. This Imp2000-centred pattern indicates that the effect of initial urbanisation background on long-term ΔNGR was conditioned by hydrothermal change, subsequent impervious surface expansion, terrain, and hydrological proximity. Compared with Model A, the Model B network therefore highlights the role of historical urbanisation context in structuring long-term NGR change. In Model C, after baseline urbanisation variables were excluded, the network became more concentrated around SRad_TS, Slp, dImp, and Dist_Riv ( Figure 8c). The strongest interaction was SRad_TS–Slp, while dImp also showed relatively strong links with SRad_TS and Slp. This pattern suggests that, once the historical urbanisation signal was removed, long-term ΔNGR was explained more by the coupling between ongoing urban expansion, terrain constraints, hydrothermal-background change, and hydrological proximity. Overall, the interaction networks indicate that the explanatory structure changed with the modelling target. Model A mainly reflected current terrain–urban coupling, Model B emphasised baseline-conditioned long-term change, and Model C shifted toward the coupling of urban-change and environmental-background variables. These patterns should be interpreted as SHAP-based conditional association structures within the fitted XGBoost models rather than as direct causal interactions. 3.6. Spatial Heterogeneity of Dominant Local Drivers To further identify where different factors exerted the strongest local influence, the dominant driver in each grid was defined as the variable with the largest absolute SHAP value ( Figure 9). The resulting dominant-driver maps revealed clear differences among the three model specifications, indicating that the local explanatory structure of cropland non-grain use varied with the modelling target. For Model A, which explained the 2020 NGR pattern, local dominance was highly concentrated in two variables ( Figure 9a,d). Slope was the dominant factor in 10,733 grids, accounting for 54.0% of valid grids, followed by Imp2020 in 6776 grids, accounting for 34.1%. Together, these two variables dominated nearly 90% of valid grids. This pattern indicates that the current NGR pattern was mainly structured by the combination of terrain constraints and contemporary impervious surface conditions. Spatially, slope-dominated grids were widely distributed across the study area, whereas Imp2020-dominated grids were more concentrated around urbanised and peri-urban zones, consistent with the terrain–urban coupling identified in the interaction network. For Model B, which explained long-term ΔNGR using both baseline and change variables, the dominant-driver pattern shifted strongly toward the initial urbanisation background ( Figure 9b,e). Imp2000 dominated 11,488 grids, accounting for 62.6% of valid grids, far exceeding other predictors. SRad_TS ranked second, dominating 4694 grids, accounting for 25.6%, while dImp, Dist_Riv, and slope showed more limited but spatially visible dominance. This pattern suggests that long-term ΔNGR was strongly conditioned by the spatial legacy of baseline urbanisation, while hydrothermal change, subsequent impervious surface expansion, hydrological proximity, and terrain conditions played secondary but locally differentiated roles. For Model C, after baseline urbanisation variables were removed, the local dominance structure changed substantially ( Figure 9c,f). SRad_TS became the most frequent dominant factor, covering 8881 grids, accounting for 48.4% of valid grids, followed by slope, Dist_Riv, and dImp. This shift indicates that when baseline urbanisation was excluded, the explanation of ΔNGR relied more heavily on broad hydrothermal-background change, terrain constraints, hydrological proximity, and urban expansion intensity. Given the coarser resolution of the TerraClimate product, this pattern more likely reflects a regional hydrothermal-background signal than grid-level climatic control over NGR change. Overall, Figure 9 confirms that dominant local mechanisms differed between the current-status and long-term-change models. The 2020 NGR pattern was mainly characterised by slope and contemporary impervious surface dominance, whereas long-term ΔNGR was more closely associated with baseline impervious surface conditions in Model B and shifted toward hydrothermal, terrain, hydrological, and urban-change signals in Model C. These spatial contrasts support the interpretation that cropland non-grain use is not governed by a single uniform driver, but by locally varying combinations of terrain, urbanisation, hydrothermal background, and hydrological context. 4. Discussion 4.1. Stable Regional Background and Localised NGR Trajectories of Cropland Non-Grain Use The NGR-decreasing (16.31%) and reversal (5.66%) trajectories should be interpreted in relation to the specific spatial development context of the CZT urban agglomeration. The Changsha–Zhuzhou–Xiangtan Metropolitan Area Development Plan emphasises a spatial structure organised around the Xiangjiang development axis, major transport corridors, three-city linkage, clustered urban development, and joint protection of the ecological green core. Under this corridor-based metropolitan structure, cropland in fringe and intercity areas is not exposed to a single directional pressure. Instead, it is repeatedly affected by urban development expectations, infrastructure-led accessibility change, agricultural profitability, cropland-use regulation, and ecological or farmland-protection constraints. Therefore, NGR-decreasing trajectories do not necessarily indicate a simple or complete “return to grain”; they may reflect the withdrawal of temporary non-grain uses, adjustment of market-oriented crop choices, farmland consolidation, or partial recovery of major-grain cultivation where regulatory and production incentives became stronger. Reversal trajectories may likewise occur where land-use expectations fluctuate between urban expansion, agricultural adjustment, and protection-oriented management. In this sense, decreasing and reversal pathways are not residual categories, but important evidence that cropland use in CZT evolved through corridor-sensitive, policy-constrained, and partly reversible adjustment processes. By classifying grid-level trajectories over a 20-year period rather than relying only on net ΔNGR values, this study extends the analytical focus from static status mapping to dynamic pathway characterisation. This distinction is important for cropland governance because stable low-NGR agricultural hinterlands, corridor-related NGR-increasing areas, NGR-decreasing areas, and reversal areas imply different management priorities. In particular, decreasing and reversal trajectories should be verified together with local land-management records and field evidence before being translated into policy categories, because they may represent different combinations of production recovery, temporary land-use adjustment, urban fringe expectation change, and data-product uncertainty. 4.2. Historical Conditioning and Mechanistic Divergence Between NGR Status and Long-Term NGR Changes A central finding of this study is that the explanatory structures of the 2020 cross-sectional NGR pattern and long-term ΔNGR differed substantially. Model A was mainly associated with terrain/soil background and contemporary urban development pressure, with slope and Imp2020 as the dominant individual predictors. In contrast, Model B was more strongly conditioned by baseline urbanisation, especially Imp2000, which provided considerably more explanatory information than concurrent impervious surface expansion. This divergence indicates that explaining where NGR is currently high requires a different analytical lens to explaining where and by how much NGR has changed over time. Therefore, cross-sectional driver analyses cannot straightforwardly substitute for longitudinal NGR change modelling. Recent studies applying interpretable machine learning to cropland non-grain or non-agricultural dynamics in China also suggest that driver interpretation can vary with the definition of the response variable and the temporal framing of land-use change [ 3, 6]. The comparison between Models B and C further supports this interpretation. The stronger predictive performance and different importance structure of Model B indicate that baseline urbanisation variables contributed additional information for explaining long-term ΔNGR. At the same time, Model C shows that urban expansion intensity, terrain, hydrothermal change, and hydrological proximity still retained explanatory signals after baseline urbanisation variables were removed. Long-term NGR changes should therefore be understood as a joint outcome of historical conditioning and ongoing urban–environmental change processes, rather than as a response to either baseline context or concurrent changes alone. The repeated importance of SRad_TS in Models B and C suggests that long-term NGR change was associated not only with urbanisation and terrain, but also with broader hydrothermal-background differences. This association should be read at the regional scale because SRad_TS was derived from the 4 km TerraClimate product, whereas NGR was analysed on 1 km grids. Thus, SRad_TS is better used to distinguish areas with different long-term environmental settings than to identify parcel- or grid-level climatic controls. Its main value is diagnostic: it indicates where hydrothermal conditions may interact with urbanisation and terrain processes and where higher-resolution data or field verification would be needed before management translation. 4.3. Nonlinear Response Zones, Coupled Interactions, and Local Mechanism Heterogeneity The SHAP dependence and interaction analyses reveal that NGR responses to dominant predictors were characterised by nonlinear transitions, conditional coupling, and spatial variation in local mechanism dominance—collectively indicating that non-grain use and NGR change in CZT are context-dependent processes rather than a uniformly driven regional trend. The interaction results further show that these nonlinear responses were embedded in different coupled-driver contexts. In the current-status model, the strong Slp–Imp2020 link suggests that contemporary urban pressure was expressed through terrain-dependent conditions. In the long-term-change models, the Imp2000-centred network in Model B indicates that historical urbanisation formed a background through which later environmental and urbanisation changes were expressed, whereas Model C shifted toward the coupling of dImp, Slp, SRad_TS, and hydrological proximity after baseline urbanisation variables were removed. This contrast reinforces the need to distinguish current NGR status from long-term NGR change: the former reflects present terrain–urban constraints, whereas the latter depends more strongly on accumulated urbanisation legacy and subsequent change processes [ 12, 42]. The dominant-driver maps ( Figure 9) demonstrate that different modelling questions generate different spatial interpretations of cropland non-grain use and NGR change: the current-status model highlights terrain and contemporary urban pressure, whereas the NGR change models reveal historical urbanisation legacy and hydrothermal-background signals. In Model A, slope and Imp2020 jointly dominated approximately 90% of valid grids. In Model B, Imp2000 dominated 62.6%, reflecting a pervasive historical conditioning imprint across the study area. In Model C, SRad_TS became the most spatially extensive dominant factor (48.4%), with slope second (27.5%). However, because SRad_TS was derived from a coarser-resolution climate product, its spatial dominance should be interpreted as a broad hydrothermal-background signal rather than as evidence of a precise local mechanism. A governance approach based on a single indicator or static threshold may therefore overlook important spatial differences in dominant mechanisms across the landscape, underscoring the analytical value of mechanism-specific spatial diagnostics for differentiating management priorities. 4.4. Trajectory- and Mechanism-Informed Cropland Management First, stable low-NGR agricultural hinterlands should be distinguished from high-risk transition zones. These areas retained relatively stable major-grain functions over the study period and should not be treated as equivalent to rapidly changing urban fringe cropland. For these grids, the main priority is to maintain grain-production incentives, agricultural infrastructure, irrigation and drainage conditions, and the continuity of cultivated land management. Excessively intensive non-grain control in these relatively stable areas may misallocate governance attention
