Coupling ITO3dE model and GIS for spatiotemporal evolution analysis of agricultural non-point source pollution risks in Chongqing in China

To determine the risk state distribution, risk level, and risk evolution situation of agricultural non-point source pollution (AGNPS), we built an ‘Input-Translate-Output’ three-dimensional evaluation (ITO3dE) model that involved 12 factors under the support of GIS and analyzed the spatiotemporal evolution characteristics of AGNPS risks from 2005 to 2015 in Chongqing by using GIS space matrix, kernel density analysis, and Getis-Ord Gi* analysis. Land use changes during the 10 years had a certain influence on the AGNPS risk. The risk values in 2005, 2010, and 2015 were in the ranges of 0.40–2.28, 0.41–2.57, and 0.41–2.28, respectively, with the main distribution regions being the western regions of Chongqing (Dazu, Jiangjin, etc.) and other counties such as Dianjiang, Liangping, Kaizhou, Wanzhou, and Zhongxian. The spatiotemporal transition matrix could well exhibit the risk transition situation, and the risks generally showed no changes over time. The proportions of ‘no-risk no-change’, ‘low-risk no-change’, and ‘medium-risk no-change’ were 10.86%, 33.42%, and 17.25%, respectively, accounting for 61.53% of the coverage area of Chongqing. The proportions of risk increase, risk decline, and risk fluctuation were 13.45%, 17.66%, and 7.36%, respectively. Kernel density analysis was suitable to explore high-risk gathering areas. The peak values of kernel density in the three periods were around 1110, suggesting that the maximum gathering degree of medium-risk pattern spots basically showed no changes, but the spatial positions of high-risk gathering areas somehow changed. Getis-Ord Gi* analysis was suitable to explore the relationships between hot and cold spots. Counties with high pollution risks were Yongchuan, Shapingba, Dianjiang, Liangping, northwestern Fengdu, and Zhongxian, while counties with low risks were Chengkou, Wuxi, Wushan, Pengshui, and Rongchang. High-value hot spot zones gradually dominated in the northeast of Chongqing, while low-value cold spot zones gradually dominated in the Midwest. Our results provide a scientific base for the development of strategies to prevent and control AGNPS in Chongqing.


Study area. Chongqing is located in southwestern China and is one of China's four municipalities directly
under the central government. The current fertilizer application level in Chongqing is 411 kg ha, while the internationally recognized safe upper limit is 225 kg ha. More seriously, the fertilizer use rate is only about 35%. Pesticide use intensity is 9.5 kg ha, with a use rate of only 30%. Livestock and poultry stocks are currently increasing significantly. In addition, the hilly and mountainous areas in Chongqing account for 94% of the total area, while the cultivated land area with a slope greater than 25 degrees accounts for 16% of the total cultivated land area, which is 11 percentage points higher than the national average level. Rainfall in Chongqing is large and concentrated, and the area subjected to soil erosion accounts for 48.6% of the total area 11 . Based on the schematic illustration of the land use types in Chongqing from 2005 to 2015, the artificial surface area displays significant external diffusion in the main urban area of Chongqing, the western region of Chongqing, Changshou County, Kaizhou County, and Wanzhou County. In contrast, the area of paddy fields shows a significantly decreasing trend, while the area of dry land is relatively stable; farmland is mainly distributed in the western region as well as in Wanzhou County, Kaizhou County, Liangping County, and Dianjiang County in the northeastern region. According to the data analysis results, the water area accounts for 1.59-1.96% of the coverage area of Chongqing. The area of paddy fields decreased from 14.8% in 2005 and 2010 to 7.64% in 2015, while the area of dry land is relatively stable at 25%. The artificial surface area has increased from 1.42 to 3.18%, while the proportion of forested land showed a steady increase (Fig. 1). Overall, the farmland area of Chongqing showed a significantly decreasing trend, with a reduction from 39.35% to 33.03%. Farmland is the main source of AGNPS. To sum up, the actual regional conditions of land use types, topography, climate, fertilizer and pesticide application lead to the aggravation of AGNPS in Chongqing. www.nature.com/scientificreports/ Research methods. The risk assessment of AGNPS in Chongqing was carried out by constructing the ITO3dE model, and the specific content was analyzed by using GIS technology. The overall framework of the study was as follows (Fig. 2).
ITO3dE model building. The ITO3dE model was built from three dimensions of "Input-Translate-Output" (Table 1). Here, the I dimension (A 1 ) includes fertilizer use intensity (I 1 ) 22 , pesticide use intensity (I 2 ) 23 , and livestock intensity (I 3 ) 24 , the T dimension (A 2 ) includes erosion caused by rainfall (I 4 ), slope length and gradient  www.nature.com/scientificreports/ (I 5 ), soil erodibility (I 6 ), sloping field (I 7 ), and distance from water area (I 8 ), and the O dimension (A 3 ) includes water quality (I 9 ) 25 , water capacity (I 10 ) 26 , water network density (I 11 ) 27 , and degree of paddy field retention (I 12 ). The computational formula of positive indices is as follows: The computational formula of negative indices is as follows: In the above equations, I i denotes the calculation result of a certain index of i, C i is the true value of the i index, and E i is the reference value. The calculation method, reference value and source of reference value involved in ITO3dE model calculation were described in Table 1. The calculation results for each index can be classified into five grades of no risk, low risk, medium risk, high risk, and extremely high risk, corresponding to the values of ≤ 0.7, 0.7-1.0, 1.0-3.0, 3.0-5.0, and ≥ 5.0. This grading standard refers to the relevant technical planning of China's agricultural sector 28 , and the grading standard of this study was written based on Chongqing local standard "Specification for evaluation risk of AGNPS in Chongqing".
Delphi method. The weights of the evaluation indices of AGNPS risks were determined by the Delphi method 18 . We distributed 120 questionnaires to the leaders of related departments engaged in agricultural and rural work, environmental work, and water conservancy programs, as well as to experts at colleges, universities, and research institutes, and to technical personnel at the grass-roots level; in total, we received 115 valid questionnaires. After the first round of weights assignment, we returned the calculation results to the experts for adjusting the index weights, representing the second round of weights assignment. According to the evaluation opinions of the experts, we again revised the evaluation indices and conducted the third round of weights assignment to obtain the final weights results. Based on these, we could obtain the multifactor comprehensive evaluation relational expressions as follows: In the above equations, the letter A represents the comprehensive risk result of AGNPS, while the other letters have the same meanings as in Table 1, and indicating the calculation results of each index.
Spatiotemporal transition matrix of risk index. Transition matrices are widely used in analyzing the spatiotemporal variations and can clearly exhibit the variations of risk indices at different periods 29 . In this study, we assigned the five risk grades to the values of 1, 2, 3, 4, and 5, respectively, and subsequently employed the following formula to analyze the risk status: where B 2005 , B 2010 , B 2015 is the operation layer and B M is the result layer. In the calculation results, "123" represents the risk transition process in a certain region from no risk in 2005 to low risk in 2010 and to medium risk in 2015; the remaining results can be interpreted in the same manner.
Kernel density analysis of risk index. Kernel density is mainly used to analyze the spatial concentration of an event and is widely used in the distribution of buildings, schools, and criminal activities 30 . The kernel density can reflect the concentration degree and agglomeration location of AGNPS above a high risk level. Kernel density estimation is a spatial analysis method based on nonparametric testing 30,31 . The basic idea of kernel density estimation for certain elements is to assume that there always exists an element intensity at an arbitrarily position within a particular region. The density intensity of geographical elements in a specific location can then be estimated by measuring the number of elements per unit area. By determining the element density at different locations and the spatial differences, the relative concentration degree of the spatial distribution of elements can be depicted, and the hotspot distribution regions can be identified. Subsequently, using the kernel density analysis tool in the ArcGIS software, we can analyze the grids at risk grades of high risk and extremely high risk and explore the extreme point position of AGNPS in our study area, Chongqing.
Getis-Ord Gi* analysis. The Getis-Ord Gi* analysis is widely used in crime analysis, epidemiology and economic geography, and is used to identify spatial gathering of high values (hot spots) and low values (cold spots) with statistical significance 32 . In Getis-Ord Gi* analysis, the z score, p values, and confidence intervals (Gi_Bin) are employed to create a new output class for each element in the input element class. Here, the z score and p values can help to judge whether the null hypothesis can be rejected, while the Gi_Bin field is used to identify statistically significant hot and cold spots. The elements in the confidence interval of [+ 3, − 3] have a statistical significance with a confidence level of 99%, while those in the confidence interval of [+ 2, − 2] have a statistical significance with a confidence level of 95%, and those in the confidence interval of [+ 1, − 1] have a statistical (1) www.nature.com/scientificreports/ significance with a confidence level of 90%; when the element gathering of Gi_Bin field is 0, there is no statistical significance.

Statement.
We confirm that all experimental protocols were approved by College of Resources and Environment of Southwest University, all methods were carried out in accordance with relevant guidelines and regulations of questionnaires, and confirming that informed consent was obtained from all subjects or, if subjects are under 18, from a parent and/or legal guardian.

Results
Results of risk assessment by the ITO3dE model. The results in the I dimension show that, overall, the distribution was high in the west and low in the northeast and the southeast in all three periods ( Fig. 3, I, II, III; Table 2), and this tallies with the topography of Chongqing. The northwestern and central regions of Chongqing are mainly hilly and slightly mountainous, while the southeastern and northeastern regions represent the Dabashan Mountain system and the Daloushan Mountain system, respectively. Thus, farmland in Chongqing is mainly distributed in the western regions as well as in regions with extensive flat areas, such as Dianjiang and Liangping. Some regions in Dianjiang, Yongchuan, Dazu, Shapingba, Wansheng, and Jiangbei show relatively high risks, but the risk level is still medium. Hence, it can be concluded that the risk level in the I dimension during 2005-2015 is, overall, not high. Considering there are too many single-factor graphs, we omitted these graphs, but provide the following description: Among the three single factors, I 1 has the highest value, and I 1 and Liangping, showed a certain decrease in 2015, but the risk levels of some regions such as Pengshui, Qianjiang, and Xiushan showed an increasing trend. The risk grade of I 2 was relatively lower than that of I 1 , but overall, the spatiotemporal variation trend was consistent with that of I 1 , except for the increasing trend of the risk level of www.nature.com/scientificreports/ Qianjiang. Basically, the risk grade of I 3 was zero; only the risk level of Bishan was in the medium risk status, while those of Hechuan and Fengdu were low. Spatially, the results in the T dimension presented, overall, an opposite distribution pattern when compared to the I dimension, that is, with low levels in the western regions and high levels in the northeastern and southeastern regions (Fig. 3, IV, V, VI; Table 2). The annual differences in the T dimension data are mainly determined by the variations in the factors I 4 and I 7 , which showed relatively higher risk levels in all three periods. The values of I 4 in the years 2005, 2010, and 2015 were 1.42-5.78, 0.84-6.12, and 0.14-6.93, respectively, while those of I 7 were 0, 0-5.38, and 0-5.06, respectively. Because Chongqing is a typical mountainous city with purple soil 33 , high-risk and extremely high-risk regions, I 5 and I 6 , are widely distributed across the city. In addition, due to the introduction of the factor I 8 , the water areas had a higher risk level, which is consistent with the actual situation of AGNPS.
The results in the O dimension showed a smaller interannual variation, with a low overall risk level (Fig. 3, VII, VIII, IX; Table 2). The O dimension levels were mainly affected by the spatial changes in the paddy field area. As mentioned above, during the 10 years, the area of paddy fields in Chongqing was nearly reduced by half, which led to the decrease in the spatial distribution of I 12 and an increased risk in counties such as Kaizhou, Fengjie, Liangping, and Changshou. Spatially, Yongchuan, Shapingba, Bishan, Dianjiang, Changshou, and Kaizhou showed higher risk levels, and the risk levels of Kaizhou, Fengjie, Wanzhou, Liangping, and Changshou showed a significantly increasing trend. The high risk values of I 9 were mainly distributed in Yongchuan, Shapingba, Jiangbei, Changshou, Dianjiang, and Liangping, with Shapingba showing the highest value of 3.75, while Chengkou, Wushan, Fengjie, Shizhu, and Xiushan had lower values. The high risk values of I 10 were mainly distributed in the western regions and were below the medium risk levels. The risk values in 2010 were higher than those in 2005 or 2015, but did not surpass 3.0, and the high values were mainly distributed in the western regions as well as in Dianjiang, Wanzhou, and Liangping. The risk values of I 11 were all below 3.0, and the highest value of 2.78 was found for Fengjie; higher values were mainly distributed in the northeastern and southeastern counties. The high risk values of I 12 were mainly distributed in the northeastern and southeastern counties, which mostly have only small areas of paddy fields.   www.nature.com/scientificreports/ Spatiotemporal change results of risk by transition matrix analysis. By assigning no risk, low risk, and medium risk levels with 1, 2, and 3, respectively, in GIS, we can obtain the spatiotemporal transition matrix according to the formula of the transition matrix. Figure 5 shows the spatiotemporal transition situation of the AGNPS risk evaluation in Chongqing. Basically, high levels show no changes, and the proportions of 'norisk no-change' , 'low-risk no-change' , and 'medium-risk no-change' situations were 10.86%, 33.42%, and 17.25%, respectively, accounting for 61.53% of the total area of Chongqing. Among these, the 'no-risk no-change' situation was mainly distributed in Rongchang, the east of Nanchuan, Shizhu, Pengshui, and Qianjiang; the 'low-risk no-change' situation was widely distributed in Wulong, the southeast of Fengdu, the south of Nanchuan, and the northeastern counties of Chongqing, while the 'medium-risk no-change' situation was mainly distributed in Shapingba, Yongchuan, Dianjiang, the north of Nanchuan, and Kaizhou. During 2005-2015, the proportions of risk increase, risk decline, and risk fluctuation were 13.45%, 17.66%, and 7.36%, respectively. Risk increases mainly occurred in central Jiangjin, central Fengdu, Pengshui, Qianjiang, the midwest of Yunyang, central Liangping, Wuxi, Wushan, and Chengkou, while risk declines were mainly observed for the main urban area of Chongqing, northern Tongliang, Dazu, Youyang, and Xiushan. Risk fluctuation was concentrated in Jiangjin, Bishan, Fuling, and Youyang. Figure 6 shows the kernel density analysis results of the medium-risk regions. As seen in the figures, the peak values of the kernel density at these three periods were all around 1,110, suggesting that the maximum gathering degree of medium-risk pattern spots basically showed no changes. The spatial distribution of kernel density at these three periods showed a consistent trend, but the distribution differences at different periods were significant. In 2005, medium-risk regions were mainly concentrated in Shapingba, southern Dazu, central Yongchuan, eastern Beibei, Dianjiang, central Kaizhou, northwestern Shizhu, northern Nanchuan, central Wanzhou, southwestern Zhongxian, and southeastern Xiushan, while in 2010, such regions mainly occurred in Shapingba, eastern Jiangjin, southeastern Beibei, northern Nanchuan, northeastern Changshou, Dianjiang, northern Fuling, northern Fengdu, northeastern Shizhu, northeastern Liangping, central Kaizhou, Wanzhou, northeastern Pengshui, and eastern Xiushan. In 2015, medium-risk regions were mainly concentrated in Shapingba, Yongchuan, central Jiangjin, northwestern Nanchuan, northeastern Beibei, Dianjiang, Liangping, the junction of Fuling and Fengdu, central Kaizhou, northern Yunyang, eastern Pengshui, southeastern Qianjiang, and central Xiushan.

Results of risk concentration degree by Kernel density analysis.
To further explore the distribution of regions with the high-risk gathering zones (Table 3), we conducted a separate analysis on the regions with kernel density values higher than 1,000 (the kernel density values of these regions were divided into 10 grades with equal intervals, and the 10th grade had values from 1,000 to 1,110).  www.nature.com/scientificreports/

Results of hot and cold spots by Getis-Ord Gi* analysis.
Applying Getis-Ord Gi* analysis is helpful to clearly identify high-value hot spots (Hot Spot-99% Confidence) and low-value cold spots (Cold Spot-99% Confidence). Figure 7 shows the Getis-Ord Gi* analysis results; the overall variation trends of high-value hot spots and low-value cold spots were consistent in all periods, with significant distribution differences. The regions located in the high-value hot spot zones in all three periods were Yongchuan, Shapingba, Dianjiang, Liangping, northwestern Fengdu, and Zhongxian, while those located in the low-value cold spot zones were Chengkou, Wuxi, Wushan, Pengshui, and Rongchang. Throughout the 10 years, the high-value hot spot zones showed significant diffusion in Fengjie, Yunyang, Kaizhou, central Qianjiang, and northern Nanchuan, while the low-value cold spot zones showed significant diffusion in some parts of the midwestern counties such as central Fuling and southern Yubei. These high-value hot spots or low-value cold spots were mainly distributed in the above-mentioned regions and their surrounding areas and showed significant "gathering trends". The spatiotemporal variation trend of the distribution of these high-value hot spots or low-value cold spots can reflect the variation tendencies of hot spots or cold spots in different regions. Over time, the high-value hot spot zones gradually migrated towards the northeastern counties of Chongqing, while the low-value cold spot zones in the midwestern counties presented an obvious diffusion trend. The low-value cold spot zones in the northeastern regions gradually decreased, while those in the southeastern regions tended to become more fragmented. These results indicate that the high-value hot spot zones gradually dominated the northeastern regions, while the lowvalue cold spot zones gradually dominated the midwestern regions.

Discussion
The I, T, and O dimensions could better comprehensively analyze the risk situation of AGNPS. There were many achievements in AGNPS risk research using different methods. For instance, Huang has analyzed the risk of NPS in Taihu Lake from multiple perspectives, and the method was only applicable to small-scale studies 35 . In addition, it focused on the distribution of pollution sources, and on waste disposal, but did not consider the process from pollution sources to water bodies. Blankenberg has analyzed the relationship between wetland and pesticide concentration in streams from the perspectives of pesticide spraying and changes in pesticide concentration in water bodies 36 . The research focused on the relationship between wetland and pesticide concentration in small watersheds, while the land use types, topographic conditions, rainfall, and other factors that affect concentration changes were not considered. This research therefore comprehensively considers input factors, translate factors and output factors. The chemical fertilizers, pesticides, livestock and poultry factors selected in the I dimension of this study were recognized as the main sources of AGNPS [37][38][39] . In addition, the factors of the T dimension such as rainfall, slope length and gradient, soil erodibility, sloping field, and distance from water area have gradually increased in recent Table 3. Distribution of regions with high-risk gathering zones.

Years
The high-risk gathering zones www.nature.com/scientificreports/ years. For instance, Zhang has analyzed the effects of rainfall intensity and slope on the loss of suspended solids and phosphorus in runoff 40 , which showed that the slope difference of sloping farmland was also an important factor affecting pollutant migration 41 , and the distance from the water reflected the difficulty degree of pollutants entering the water body 42 . The influence factors of the O dimension increased the water network density and the degree of paddy field retention factor compared with other studies. Most of the studies mainly focused on risk the analysis of AGNPS from water quality and water capacity 43,44 . In Wang et al. 34 used the minimum cumulative resistance model to calculate the risk of ANSP of cultivated land in the Three Gorges Reservoir area, and believed that the topography, hydrology, soil and vegetation had a great impact on the risk of ANSP. In fact, high water network density and paddy field ratio could effectively reduce the pollution. In this study, the I, T, and O dimensions were considered comprehensively, and factors affecting AGNPS risk were added to the large-scale framework.
Geographical methods were well applied in the field of agriculture. The use of geographic methods in this study makes the results of AGNPS risk analysis more intuitive and reflects the temporal and spatial changes of research results from different perspectives. The spatiotemporal transition matrix quantifies the variation in the different risk levels and can accurately identify the regions with increasing, decreasing, or unchanged risk levels. The spatiotemporal transition matrix analysis was one of the most widely used methods in geographic research and often used to study the change of the same spatial event over time. For instance, Shawul and Chakma have used this method to analysis the conversion of land use types in different periods 45 . The kernel density analysis could effectively identify the spatial agglomeration area of an event and was widely used in geographic, medical event analysis and case analysis 46 . The advantage of the Getis-Ord Gi* analysis is to judge the hot spots of high and low values at the same time and to analyze the evolution trend of hot spots of high and low values. This method was usually applied in studies that require the analysis of high-value and low-value changes simultaneously 47 . Zhang et al. 48 used the hot spot analysis method in geography to analyze the spatial and temporal pattern of ANSP in Chongqing section of the Three Gorges Reservoir area, and achieved good results. In our research, we need to consider both high-risk and low-risk changes in AGNPS.

Land variations influence AGNPS risks. The analysis of land use changes indicates that the land area
in Chongqing has changed greatly during the 10-year-period from 2005 to 2015. In particular, the paddy field area was reduced by nearly 50%. Our observations are in agreement with previous findings 49 . The regions with declining farmland areas were mainly located in the urban areas of Chongqing and in the surroundings of the town, which goes hand in hand with the economic development and the urban expansion of Chongqing in recent years, reflecting the direct influences of urban development on AGNPS. This tallies with the decrease in chemical fertilizer use intensity and pesticide use intensity, as farmland reduction inevitably lead to a decrease in the application of these substances 50 , complying with the "one regulatory, two reduction, three basic" policy promoted by the Chinese government. A previous study has analyzed the relationship between the riparian forest buffer zone and pollution mitigation in Chesapeake bay watershed and found that the wide buffer zone could effectively reduce water pollution 51 . Other authors have analyzed the relationship between climate, land use change, and water quality in the Pike river watershed and found that land use changes had an impact on water quality 52 . The water quality factors in this study also reflected that the water quality in the southeast and northeast areas with high forest coverage was better, indicating that land use change had a great impact on water quality, especially on the cultivated land, which was closely related to the use of fertilizers and pesticides.
In this study, we analyzed the spatiotemporal variation in AGNPS in three dimensions of "Input-Transition-Output", with the conclusion that the risk levels in the Input and Output dimensions were lower, while that in the Translate dimension was higher. Our conclusion accords with the characteristics of Chongqing as a mountainous city; the widely distributed purple soil, the larger differences in slope gradient and slope length, and the widely distributed sloping fields are indeed conducive to the generation of AGNPS 53 . PENG et al. 54 analysis of the driving factors of ANSP in the Three Gorges Reservoir Area (Chongqing Section) showed that land use type, agricultural production and living intensity have a greater impact on ANSP, which was consistent with our research conclusion. From 2005 to 2015, the counties with higher risk levels in all three periods were basically those national basic farmland demonstration counties (Jiangjin, Dazu, Tongnan, Tongliang, Liangping, Dianjiang), with higher risk levels in the Input and Output dimensions. These areas therefore require planting activities under the conditions of ensuring soil safety, water quality safety, and food security.

Recommendations.
Focus on risk increase area in the analysis of spatiotemporal transfer matrix. In 62% of our study area, the risk levels remained unchanged, indicating that the influencing factors of AGNPS in most regions of Chongqing are relatively stable. For about 14% of the study area, the risk levels showed an increasing trend, and these regions, especially food security zones such as Jiangjin and Liangping as well as important ecological shelter zones of the Three Gorges Reservoir areas such as Fengdu, Yunyang, Wuxi, and Wushan, need to be considered in management programs. AGNPS in these areas might seriously threaten the ecological security of the Yangtze River economic belt. Regional managers should control the use of chemical fertilizers and pesticides and keep them within the international limits 55 . As Kesavan stated in their research, the use of chemical fertilizers and pesticides results in enhanced productivity over short periods, but leads to the degradation of soil health, freshwater, and biodiversity in the long term. In addition, the layout of farms in prohibited and restricted aquaculture areas should be strictly controlled, and strict total quantity control should be implemented in suitable areas. The prohibited and restricted livestock breeding areas in all districts and counties have been revised by Chongqing Environmental Protection Bureau in 2019, aiming to better manage the layout of livestock farms and reduce the risk of water environment pollution. www.nature.com/scientificreports/ Focus on high-risk agglomeration area in nuclear density results. The high values of kernel density changed only slightly over time, indicating that the maximum scope of gathering zones showed no significant changes. The spatial position variation of kernel density at different periods is closely related to the changes of the factors in the three dimensions. The main gathering zones of medium risk levels showed some spatiotemporal differences in the three periods, but there are high-risk gathering zones in Shapingba, Yongchuan, Nanchuan, Dianjiang, Liangping and Kaizhou, and Xiushan. These phenomena are closely related with the agricultural development degree of these counties. The analysis results of such zones after further grading more clearly display the distribution of high-risk gathering zones, which, to a certain extent, reflects the distribution of medium-risk pattern spots in the vicinity of the high-risk gathering zones. These regions, such as Youting town in Dazu, Baijia town and Chengxi town in Dianjiang, the junction of Gaofeng town and Gangjia town, and Linjiang town in Kaizhou, require significant attention from the respective prevention and control departments. The cultivation system in high-risk agglomeration areas should be rationally adjusted. The crop rotation system could be adopted to restore soil fertility while controlling the amounts of chemical fertilizers and pesticides. In addition, increasing the vegetation coverage of ridges could be considered, and previous studies have proved that these measures could effectively reduce AGNPS 56 . Chongqing government has constructed river regulation works and ecological corridors in several districts and counties of the Yangtze River Basin, with the purpose of reducing soil erosion near water bodies, increasing the interception capacity of vegetation near water bodies to AGNPS, and improving the prevention and control ability of AGNPS.
Focus on high-value hot spots area in Getis-Ord Gi* analysis. The results of the Getis-Ord Gi* analysis clearly reflect the spatiotemporal evolution situation of high-value hot spot zones and low-value cold spot zones in the three periods in Chongqing. Over time, the dominance of high-value hot spot zones in the Midwest gradually became lower than that of the low-value cold spot zones, with a concentrated distribution in Dianjiang, Zhongxian, Fengdu, Wanzhou, Liangping, Kaizhou, Yunyang, and Fengjie. In the Midwest, the coverage areas of high-value hot spot zones and low-value cold spot zones were generally neck and neck; by contrast, in the southeast, low-value cold spot zones dominated. This trend indicates that zones with high AGNPS in northeastern Chongqing have the tendency of further gathering and diffusing, and the high-risk areas migrate toward northeastern Chongqing. Against this background, these regions deserve special attention by decision makers. Due to the considerable amount of sloping farmlands in these areas, the cultivation of such farmland should be strictly controlled 57 . In addition, soil and water conservation measures should be strengthened in these areas. In recent years, Chongqing Planning and Natural Resources Bureau has changed a large number of slope farmland into non cultivated land in order to reduce the carrying capacity of soil and water loss to pollutants, and reduce pollution risk from the source and transmission process of AGNPS. AGNPS risk assessment is different from pollution load measurement. The accuracy of the risk results does not depend on the experimental or measured data, and may not be positively correlated with the measured data. This is similar to the multi index system assessment of NPS carried out by Huang for Taihu Lake. And according to the risk values calculated from multi-angle indicators, the risk differences in various regions are analyzed and countermeasures are proposed 35 . For example, if the I dimension value of an area is high, but the T dimension is low, that is, the measured data value is high, this does not automatically mean a high risk. The results of risk assessment are based on the reliability of the data source, factor weight and grading reference value of risk factors. It is based on the historical data to judge the trend of regional pollution risk, so as to prevent and control the risk, and this is the advantage of risk assessment in pollution prevention. Therefore, in the future, we should increase investment in data monitoring and acquisition of I, T and O dimensions and scientifically analyze the differences in the reference values of each factor in different regions.

Conclusions
We built an ITO3dE model, analyzed the land use change during 2005-2015 in Chongqing by the combination of various spatial analysis functions including the spatiotemporal transition matrix in GIS, kernel density analysis, and Getis-Ord Gi* analysis, and clarified the influences of farmland change on AGNPS risks. The spatiotemporal variation of the pollution risk indicates that the T dimension had the highest risk level in different regions of Chongqing. This phenomenon was closely related with the widely distributed purple soil, the large differences in slope gradient and slope length, and the widely distributed sloping fields. The risk evaluation results obtained via spatiotemporal transition matrix, kernel density analysis, and Getis-Ord Gi* analysis specifically quantified the regions that require close attention in the prevention and control of AGNPS from the perspectives of risk level variation, the distribution of risk highly-concentrated zones, and the shift of high-value hot spot zones and low-value cold spot zones; we therefore provide data to support land use planning strategies as well as measures to prevent and control AGNPS.
Taken together, there would be significant influences on the Input and Output dimensions of AGNPS that promote the "one regulation, two reductions, and three basics" policies promulgated by the Chinese government, indicating areas suitable for livestock breeding in Chongqing and implementing regulatory requirements in 2018. The implementation of related policies would effectively reduce the risk levels of these two dimensions. In this study, the results suggest that there is a higher risk level in the Translate dimension. Among its influencing factors, field slope is important and can be improved by land use planning or artificial handling. Hence, related departments should focus on the farmland protection in the high-risk gathering regions, avoid farmland reclamation as far as possible, and reduce the proportion of sloping fields. At the same time, the migration of pollutants into water bodies is mainly affected by the degree of soil erodibility, so there would a positive impact that strengthening the control measures of regional soil erosion on reducing AGNPS.