Assessing the impact of land surface temperature on urban net primary productivity increment based on geographically weighted regression model

Urbanization had a huge impact on the regional ecosystem net primary productivity (NPP). Although the urban heat island (UHI) caused by urbanization has been found to have a certain promoting effect on urban vegetation NPP, the factors on the impact still are not identified. In this study, the impact of urbanization on NPP was divided into direct impact (NPPdir) and indirect impact (NPPind), taking Kunming city as a case study area. Then, the spatial heterogeneity impact of land surface temperature (LST) on NPPind was analyzed based on the geographically weighted regression (GWR) model. The results indicated that NPP, LST, NPPdir and NPPind in 2001, 2009 and 2018 had significant spatial autocorrelation in Kunming based on spatial analytical model. LST had a positive impact on NPPind in the central area of Kunming. The positively correlation areas of LST on NPPind increased by 4.56%, and the NPPind caused by the UHI effect increased by an average of 4.423 gC m−2 from 2009 to 2018. GWR model can reveal significant spatial heterogeneity in the impacts of LST on NPPind. Overall, our findings indicated that LST has a certain role in promoting urban NPP.


Results
The spatial features of NPP, IS and LST. In order to assess the spatial agglomeration, global Moran's indexes of NPP, land surface temperature (LST) and impervious surface (IS) were calculated using Geoda software 39 . As shown in Table 1, the global Moran's I of NPP, IS and LST were greater than 0 (p < 0.001), meaning that they are positive spatial autocorrelation in Kunming. In other words, NPP, LST and IS among the city areas tend to gather in space. However, the global Moran's I only shows whether there is agglomeration in the city, but it does not reveal in which area the agglomeration occurs. To further understand the evolution in local spatial distribution, the local indicators of spatial association (LISA) map of NPP, IS and LST was calculated (Fig. 1).
The local Moran LISA maps ( Fig. 1d-f) revealed that the NPP mostly were concentrated in the LL (low-low) and HH (high-high) clusters areas. Specifically, LL clusters of NPP distributed mainly inside the city and the Chenggong district, a new investment and development zone. From 2001 to 2018, the LL clusters of NPP gradually expanded to Chenggong and airport area. Interestingly, a shrinking trend of the LL clusters was observed in 2018, which had no more concentrated areas than in 2009, especially in Chenggong district. Because of the intensification of urban development, the value of regional NPP is similar to the surrounding grid points and shows insignificant characteristics. The clusters of NPP distributed mainly in north, northwest and northeast regions, where more vegetation zones. Correspondingly, the clusters of NPP could be an index for the urbanization. HH clusters indicated less vegetation inside the urban area, while LL clusters indicated more vegetation in the forest (Fig. 1a-c). In addition, HL (high-low) clusters appeared in or around the LL clusters, representing the Table 1. The global Moran's I of NPP, IS and LST in 2001. GMI denotes the globe Moran's index. Z is the standardized value which denotes the Z-Scores. P is the significance test level value which denotes the P-values. NPP is net primary productivity. IS is impervious surface. LST is land surface temperature. www.nature.com/scientificreports/ existence of vegetation in the city. HH clusters of IS ( Fig. 1g-i) mainly gathered in the city area corresponding to true color image (Fig. 1a-c). While the distribution of LL clusters of IS was approximately consistent with the HH clusters of NPP. LST is geographically concentrated in HH clusters areas ( Fig. 1j-l). And HH clusters of LST is mainly concentrated in the central, northeast and south of Kunming. The HH clusters of LST gradually expanded with the development of city through time.
In general, the HH clusters of IS, LL clusters of NPP and the HH clusters of LST were roughly distributed in city area of Kunming. The urban expansion of Kunming can be obtained from the HH clusters of IS ( Fig. 1g-i). From 2001 to 2018, the urban area gradually expanded to the airport, Chenggong district and the vicinity of Dianchi Lake. At the same time, urbanization has led to the decrease of vegetation in urban newly expansion   (Fig. 2) further illustrates the local spatial distribution geographically. The NPP dir mostly concentrated in the LL and HH clusters. LL clusters mainly concentrated in the newly expanded urban areas such as Chenggong district, representing the most obvious area of direct impact of NPP. From 2009 to 2018, the area of the LL clusters has expanded significantly and is consistent with the HH clusters of IS (Fig. 1), which proves that newly expanded urban areas mainly reflect the NPP dir because of the transformation of LULC. On the contrary, with the increase of urban intensification, NPP ind expanded from the central area to the surrounding area in the city. The indirect impact of NPP was mainly reflected in the better growth of urban vegetation within city, so the HH clusters of indirect impacts mostly are concentrated in urban areas. Corresponding to LISA cluster map of LST in 2009 and 2018, the HH clusters of LST and NPP ind have roughly similar distribution ranges.
Regression analyses of OLS model. LST is negatively correlated at the 1% significance level ( Table 3), indicating that UHI effect will lead to decreased indirect impact of NPP on the entire city area of Kunming. But it is inconsistent with our argument that UHI effect will promote the growth of urban vegetation. Because the traditional OLS regression method assumes that there are no heterogeneous differences in Kunming city, and estimates the "global" parameters of the explanatory variables. Specifically, the forest areas in the north, northwest and northeast of Kunming are also included, but the UHI effect in these areas is not obvious and may not promote vegetation growth, which lead to negative correlation between LST and NPP ind in the whole area. From 2009 to 2018, the correlation coefficient increased from − 6.019 to − 2.994. The negative correlation of LST on NPP ind became smaller in the whole areas.
Regression analyses of GWR model. We applied the GWR method to explore the spatial heterogeneity of the LST on NPP ind of Kunming. A summary of estimated coefficients is given in Table 4. The standard residual of coefficient estimation between − 2.5 and 2.5 is considered as the high reliability 37 . Most of the standard residual of coefficient estimates in the study area are reliable with the proportion is 98.8% and 97.9% in 2009 and 2018 respectively ( Table 4).
As displayed in Fig. 4, the estimated coefficients of the LST on NPP ind range from − 27.488 to 8.971, and from − 9.076 to 10.742 in 2009 and 2018. Respectively, implying that the impact of LST varies greatly from region to region. The positive coefficients are mainly concentrated in the central city shown in red and orange in Fig. 4, which is consistent with the areas of HH clusters of LST that characterizes UHI effect. The high temperature caused by the UHI effect has positively impact on increasing the urban vegetation NPP of the central city area in a certain extent. However, low estimated coefficients were observed in the surrounding area probably due to lower LST. From the perspective of spatial changes, the areas with positive coefficients gradually expanded from the main city area to the surrounding area from 2009 to 2018. A small part of the southeast of Kunming, the interior of Chenggong district, also showed a positive correlation until 2018. As shown in Fig. 3, the positive correlation area of the study area increased from 6.30 to 10.86% from 2009 to 2018. The rapid urbanization has led to an increase in the scope of LST, which can be seen in Fig. 1. In addition, due to the promotion of LST on NPP, NPP ind increased by an average of 4.423 gC m −2 from 2009 to 2018 in Kunming ( Table 5).

Comparison of OLS and GWR .
In order to prove whether the GWR model is better than the OLS model for the results in this study, we compared the R 2 and AIC statistics respectively in the GWR and OLS model. The GWR model revealed a higher R 2 (0.496 and 0.486 in 2009 and 2018 respectively) than that of the OLS model (0.143 and 0.388 in 2009 and 2018 respectively), indicating that the GWR model is better fitted. Then we compared the AIC statistics. When the AIC value differs by more than 3, the lower the AIC value, the better the www.nature.com/scientificreports/  www.nature.com/scientificreports/ fit 40 . As shown in Table 6, the AIC value for the GWR was approximately 4314 less than the AIC value of OLS in 2009, and 5312 less in 2018. To sum up, the GWR model is better fitted than the OLS model. The analysis of the estimated coefficients of the explanatory variables in the previous section also shows that the wider the range of variable coefficients, the greater the spatial variability of the contribution and influence of LST on NPP ind , which illustrates the appropriateness of the GWR model in providing a better explanation for the local estimation.

Discussion
The urbanization has led to dramatic carbon flux fluctuations. Most of the previous studies focused on the reduction of NPP caused by urban expansion [18][19][20][21] . A typical example is that Wen et al. 41 quantified the large urban expansions occurred in most provinces of China and were accompanied with huge NPP losses. They have only considered the total NPP when study about the driving factors of NPP. If the indirect impact was not considered, we could only make the conclusion that the influencing factors brought by urbanization have reduced the urban NPP. But the NPP loss just consider land-cover replacement is inaccurate, because compensation effect  www.nature.com/scientificreports/ of climate and artificial factors will reduce the NPP loss. As shown in Fig. 1d-f, the LL clusters of NPP in cities have expanded with urban development which was in good agreement with the results of previous studies. In fact, the reduced NPP was the result of the area of vegetation being replaced mainly in newly expanded areas of city, but the remaining vegetation was much more productive after the urbanization. This phenomenon has been discovered in early studies. For example, Gregg et al. 42 reported that enhanced biomass of vegetation in New York City compared with that of rural sites. Takagi and Gyokusen 24 also found that the higher photosynthetic rate in the urban core due to the air pollutant. They all began to realize the urban vegetation had much more productive affected by urbanization, in essence, this is due to the indirect impact of urbanization on urban vegetation NPP. Then there are some studies to quantify this indirect impact. And the indirect impact was found to be able to offset the direct loss of vegetation NPP caused by replacing vegetated surfaces in the urban area 16,17 . In this study, the HH clusters of NPP ind has expanded significantly with the urbanization from the perspective of spatial distribution (Fig. 2). It showed that NPP ind in Kunming mainly existed in urban areas, especially in the old city, and it gradually expanded to the Chenggong district. The mean value of NPP ind increased by 4.423 gC m −2 from 2009 to 2018 (Table 5). This result is similar to that of Guan et al. 17 that the indirect impact can promote NPP more in the old city region. In addition, separating the indirect impact of NPP can contribute to identify the real relationship between UHI effect and NPP. The total NPP of vegetation cannot represent the growth status of urban vegetation. Because of the direct impact, even though the vegetation growth is promoted by UHI effect, urban NPP will decrease with urbanization. All these conclusions could only be drawn by separating indirect impacts of urbanization on NPP. Thus, it is necessary to separate the indirect impact of NPP, in order to figure out the real relationship between UHI effect and urban vegetation NPP. The driving factors of NPP have been examined through previous research 14,17 . Due to the complexity of urban spatial heterogeneity, the coefficients of different factors varied by land use and region under different urbanization intensities. In order to divide the areas with different urbanization intensity, Tian et al. 14 divided Beijing into four eco-regions including the Capital Core Functional Region (CCFR), the Capital Extended Functional Region (CEFR), the New Developing Functional Region (NDFR), and the Ecological Reservation and Developing Functional Region (ERDFR). Guan et al. 17 divided Kunming into old city (OC), expansion area (EA), sub-urban area (SA) and non-urban area (NA). They have taken into account the spatial patterns of the urbanization impact by dividing the buffer zone, and did regression analysis by region to analyze the indirect impact of climate factors on NPP. However, these studies rely on traditional regression models, and have failed to reveal the spatial changes in each buffer zone. The regional linear regression will result in cliff-jumping changes in the boundaries between regions and also make the relationship between different regions unable to be expressed. Compared www.nature.com/scientificreports/ with the regression analysis, GWR can accurately display the regional variation of coefficient 36 . As the results of the impact of LST on NPP ind , the results from the GWR (Fig. 4) revealed more specific information which cannot be revealed by OLS model (Table 3). GWR divided Kunming into 10,984 fishnet grids with 500 × 500 m resolution to obtain the regression coefficients of each grid locally taking into account the spatial weight matrix.
Since GWR considers the spatial information of each point, it has higher accuracy compared to sub-regional regression. The spatial distribution of regional correlation coefficients revealed the specific regions where LST promoted NPP ind . Pei et al. 20 indicated that the increases of NPP might be correlated with the effects of UHI caused by urban land development around some regions that experienced rapid urbanization, as well as the arid regions in northwest China. In this study, the NPP ind also found to be positive correlated with LST in the city of Kunming, where has experienced rapid urbanization in recent years. Guan et al. 17 also mentioned that climate condition in the old city was years ahead of that in the expansion area, as a result of higher temperature in the old city. In Fig. 4, the vegetation in the old city also has more productivity affected by the UHI effect. Compared with the previous studies 18,20 , GWR specifically clarify the areas of LST promotes NPP ind over time rather than dividing urbanization intensity regions artificially. Overall, GWR has reference significance for analyzing the driving factors of NPP in urban areas with spatial heterogeneity.
The major uncertainties of this study mainly come from CASA model and other factors that may affect NPP. First, although the CASA model is widely used for NPP inversion 4,11,17,43 , the data processes and the parameter value definition could bring some bias to the NPP results 17 , which also influencing NPP dir and NPP ind . As an input data of CASA model, the monthly total solar radiation obtained by Kriging interpolation may have error caused by inadequate meteorological data sites. Subsequent research can consider other improved models to overcome the limitation of fewer sites 44,45 . And as an important static parameter of CASA model, the maximum light use efficiency ε max has been shown to be different among different vegetation types, ecosystems and seasons. But the ε max values in this study may be more suitable for the whole country, which adopted from the study of Zhu et al. 46 . So further adjustment of ε max for urban vegetations is necessary. Then, due to the lack of field measurement of NPP data in Kunming, we used the MODIS NPP product (MOD17A3) to assess the accuracy of the NPP from CASA model. The mean estimation values of CASA and the MODIS product were extracted from 500 random points except the districts in the central city, because the NPP on building land was not estimated in the MODIS product 11 . The validation results of these points indicating that the CASA model is suitable for estimating the NPP in Kunming city (R = 0.73, P < 0.01). Secondly, this study separated the impact of urbanization on NPP into NPP ind and NPP dir , so it would be meaningful to consider other influencing factors of NPP. For the total NPP, at regional scale, variation of precipitation was found to dictate most of the inter-annual variation of NPP of the tropical region. And in the mid northern latitudes the variation of NPP seems to be attributed to the relative variability and mean of temperature 47 . So Kunming NPP may be more related to temperature. Therefore, this study provides a certain understanding of the driving mechanism of LST on NPP ind in urban areas. However, other factors that may affect the NPP, such as the urban dry island, urban rain island and urban CO 2 concentration should also be considered in subsequent studies 20,22,48 . The NPP was due to the vegetation replacement in the process of urbanization, but the changes of the type, form and maintenance of the vegetation should also be explored in further study. For example, among the vegetation in different urban structures, the vegetation in old urban areas is older and have longer growth times, which is more adaptable to the urban environment. While vegetation in new urban areas is mostly newly transplanted trees with younger age. In addition, this study divided the study area into 500 × 500 m grids, so it would be meaningful to divide into finer grids. And how to increase the time point and expand the scope of study area to 32 major cities in China are also interesting directions for future work, which will be of great significance in advancing our knowledge of urbanization and terrestrial ecosystems. In general, although some conclusions were obtained in this study, we believe that the conclusions obtained in this study are credible and valuable. But further efforts are still required for an in-depth exploration of the indirect impact of urbanization on NPP.

Conclusion
Urbanization has a huge influence on regional NPP. In this study, first, the direct and indirect impact of urbanization on NPP were separated. Then spatial variations of the correlation coefficient between LST and NPP ind were analyzed both globally and locally with the support of OLS and GWR models. The main conclusions can be summarized as follows: (1) From 2001 to 2018, the results of spatial autocorrelation analysis showed that the study area has experienced an accelerated urbanization process. Most of the urban sprawl was concentrated in Chenggong district and the airport area of northeast of Kunming. With the expansion of the city and the decrease of vegetation coverage, NPP has shown a corresponding decreasing trend and LST is higher in the urban areas. (2) After dividing the impact of urbanization on NPP into NPP ind and NPP dir , the spatial autocorrelation analysis showed that NPP ind and NPP dir differ in the spatial agglomeration area. In the process of urbanization, the scope of LL clusters of NPP dir were mainly appeared in the newly expanded urban areas, which was due to the reduction of NPP caused by land cover replacement in the new urban areas. However, the HH clusters of NPP ind were concentrated in the central area of Kunming city and gradually expanded to the surrounding areas with the urban development, which was similar to the area affected by higher LST. , as the capital of Yunnan, was selected as the study area. Kunming city located at the northern of the Dianchi Basin in Yunnan, which is a representative of low-latitude plateau cities. As shown in Fig. 5, the city has a total area of 2602.46 km 2 , an average elevation of 1900 m, and is surrounded by mountains on three sides. Affected by the circulation of the southwest monsoon and the adjustment of the water surface of Dianchi Lake, a natural environment with spring seasons, abundant rainfall, long sunshine and perennial southwest wind has been formed 49 . Located in the northern subtropical monsoon climate zone, Kunming has typical temperate climate characteristics. However, the climate shows great seasonal heterogeneity, with humid summers and dry winters, and most of the precipitation occurs during the growing season (April to October). In order to distinguish the spatial patterns of the urbanization impacts, the total urban area was shown as the Wuhua district, Panlong district, Xishan district, Guandu district and Chenggong district (Fig. 5). The Wuhua district and Panlong district areas are the earliest developed old city in Kunming, which has the highest degree of impervious surface. Then Xishan District and Guandu District have a higher degree of urbanization. And the establishment of Chenggong New District in recent years has greatly promoted the rapid development of Kunming.
Data source and preprocessing. The Landsat and moderate-resolution imaging spectroradiometer (MODIS) images used in this study were selected based on the least amount of cloud (cloud cover < 10%). Landsat Thematic Mapper (TM)/Operational Land Imager (OLI) were downloaded from the Geospatial Data Cloud (http:// gsclo ud. cn/) and the United States geological survey (USGS) science center (http:// glovis. usgs. gov/). The MODIS products were also obtained from the USGS (http:// glovis. usgs. gov/). The monthly meteorological data used in this study include average temperature, total precipitation and the total solar radiation. The climate data, were derived from 29 climatological stations around Kunming city, which were provided by the China Meteorological Administration (http:// data. cma. cn/). The digital elevation model (DEM) data was obtained from LPDAA of NASA (http:// www. gdem. aster. ersdac. or. jp/), and the vector data was extracted from the GADM database (https:// www. gadm. org/). All Landsat and MODIS images were rectified to the universal transverse mercator (UTM) projection and world geodetic system 1984 (WGS84) datum.  www.nature.com/scientificreports/ NPP estimation framework, first, to obtain the normalized difference vegetation index (NDVI), the Landsat TM/OLI images were used to extract NDVI calculated by Landsat infrared (IR) band and near-infrared (NIR) band. MCD12Q1 at spatial resolution of 500 m covering the study area closer to the observation time of Landsat images were selected as land cover classification data. To control the CASA model, the land cover types data from MCD12Q1 were further reclassified into different vegetation types, including evergreen broad-leaved forest (EBF), deciduous broad-leaved forest (DBF), evergreen needle-leaved forest (ENF), deciduous needle-leaved forest (DNF), mixed forest (MF), shrub, urban land, cropland, grass-land, unused land and water area. Then the cover classification data were resampled to 30 m. Because it is difficult to acquire meteorological data including monthly average temperature, total precipitation and the total solar radiation due to the limited number of meteorological stations, the meteorological data were used to interpolated into 30 m raster images using Kriging method 53 . Kriging method is the most widely used interpolation method. It thinks that any attribute of spatial continuously change is stochastic, and semi-variation function is used to analyze data 43 . The root-mean-square standardized of Kriging is close to 1, indicating that the standard error is accurate. In this study, the root-meansquare standardized of monthly average temperature, total precipitation and the total solar radiation are 1.8525, 0.7344 and 0.9671 respectively. The mean standardized of Kriging is close to 0, indicating that the results are unbiased. And mean standardized of monthly average temperature, total precipitation and the total solar radiation are − 0.0653, − 0.0510 and − 0.0636 respectively. These Kriging statistical errors illustrate the reliability of the interpolation results. Finally, according to the study of Zhu et al. 46 , the static parameters such as the maximum light use efficiency ε max (gC MJ −1 ) 54 were obtained.
The direct and indirect impact of urbanization on NPP. In order to obtain the direct and indirect impact of urbanization on NPP, a specific method for separating and analyzing NPP dir and NPP ind was detailed as follows 17 . First, assuming that every pixel has ideal full vegetation coverage before urbanization which defined as NPP fv . And assume that NPP fv does not change over time. Then, according to SMA model, the urbanization intensity (β) is represented by the percentage of IS of every pixel. NPP h is hypothetical NPP assuming that there is no indirect impact during the urbanization period, which is determined by the percentage of the non-urban surface (1 − β) and the NPP fv together, and are described in detail as the following: Thus, based on this concept (1), under ideal conditions, the hypothetical NPP h after urbanization in time t can be expressed as: where NPP h (x, t 1 ) is the NPP in pixel x at time t. And it is the hypothetical NPP value after urbanization, just considering the direct impact of land-cover changes; (x, t) is the urban expansion intensity; and NPP fv (x, t) is the NPP value in pixel x at time t when it has full vegetation cover. The hypothetical NPP h in time t 0 can be expressed as: And the hypothetical NPP h after urbanization in time t 1 can be expressed as: Then the change in NPP due to urbanization from t 0 to t 1 is formula (4)-(3): Since the ideal full vegetation cover pixel NPP fv does not change over time: Then, we can calculate the direct impact at time t 1 as follows: where NPP dir (x, t 1 ) is the direct impact NPP in pixel x at time t 1 . The difference between NPP h (x, t 1 ) and the NPP at t 1 should be the indirect impact of urbanization on NPP. The indirect impact at time t 1 as follows: where NPP ind (x, t 1 ) is the indirect impact NPP in pixel x at time t 1 . As NPP x, t 1 is the NPP value estimated from CASA model after urbanization in pixel x at time t 1 , bring formula (3) into formula (8) to get the indirect effect NPP: According to the formula (9)  www.nature.com/scientificreports/ Kunming city were resampled to 500 × 500 m fishnet grid. Then, using Geoda software 39 to obtain the global Moran's index and local Moran's index of the above indicators. While, global Moran's I reflects the global spatial autocorrelation between different geographical regions and local Moran's I reflects the local similarities and differences between neighboring areas 55 . The two variables of the z-score and the p-value respectively refer to the spatial correlation level between neighboring areas and its corresponding significance level 35 .
(1) Global spatial autocorrelation Global spatial autocorrelation is used to measure the spatial autocorrelation. This study used the global Moran's I and its statistical test to analysis the spatial dependence of NPP, LST, IS, NPP dir and NPP ind . It was first proposed in 1950 by the Australian statistician Parker Moran 56 . The global Moran's I statistics of spatial autocorrelation can be expressed as: where n is the total number of observation area. As Y i is the Observation value in the i-th region. W ij is the spatial weight matrix, j = 1, 2, …, n. n is the number of cross-sectional observation units. The global Moran's I values between − 1 and 1, if the Moran's I is greater than 0, it means that this attribute value has a positive correlation in the study area. In contrast, if the Moran's I is less than 0, it means that the attribute value has a negative correlation. The larger the value, the greater the degree of spatial autocorrelation.
If the calculated global Moran's index shows that this attribute value has spatial autocorrelation, the Z value can be used to test its significance. Z value calculation formula is as follows: Z value can be used to test the spatial correlation. When Z > 0 and significant, there have spatial autocorrelation. When Z = 0, there have random distribution. When Z < 0 and significant, there have negative spatial correlation.
(2) Local spatial autocorrelation Compared with global spatial autocorrelation, local spatial autocorrelation is used to calculate the spatial correlation degree between each spatial object and its neighboring objects in a region. The calculation of local Moran's I is similar to the global Moran's I. Anselin proposed the definition of the local Moran's I in 1995 57 . The calculation of local Moran's I is shown as follow: The definitions of n, Y i , W ij and n are the same as the global Moran's index in the previous section (i = 1, 2, …, n, j = 1, 2, …, n).
The LISA cluster map can reflect where the indicators are clustered in the space. Local Moran's I and LISA cluster map can be divided into four space-related patterns, namely high-high cluster (HH), low-low cluster (LL), low-high cluster (LH), high-low cluster (HL). HH clusters indicates that the high-value region is surrounded by other high-value regions. LL clusters represents a low-value region is surrounded by other low-value regions. Similarly, LH clusters indicates that the low-value region is surrounded by other high-value regions, and HL clusters represents that a high-value region is surrounded by other low-value regions 36,55 . Geographically weighted regression model. The geographically weighted regression model (GWR) can be used to examine the spatial heterogeneity of the effect of LST on NPP ind . The geographically weighted regression model (GWR) is a local spatial technique that can be used to examine the spatial variabilities of regression parameters and model performance. Compared with the ordinary least squires (OLS) model, the GWR model carries out separate regressions at each location, considering only other observations within a specific distance to that location 35 . The model can be expressed as follow 37 : where y is the dependent variable. x is the independent variable. k is the number of dependent variables. u i , v i represents the geographic coordinates of the i-th unit. 0 u i , v i is a constant term, k u i , v i is the local