Modeling of Forest Communities’ Spatial Structure at the Regional Level through Remote Sensing and Field Sampling: Constraints and Solutions

This study tests modern approaches to spatial modeling of forest communities at the regional level based on a supervised classification. The study is conducted by the example of mapping the composition of forest communities in a large urbanized region (the Moscow Region, area 4.69 million hectares). A database of 1684 field descriptions is used as sample plots. As environmental variables, Landsat spectral reflectances, vegetation indices (5 images), digital elevation model and morphometric parameters of the relief, 54 layers in total, are used. Additionally, the Palsar-2 radar dataset is included. The main mapped units are formations and groups of associations identified on the basis of the ecological-phytocoenotic classification. Formations and groups of associations are similar in semantics and principles of allocation to units of forest typology. It is shown that the maximum entropy method has a wide range of applications, in particular, for mapping the typological diversity of forest cover. The method is used in combination with geographically structured spatial jack-knifing, spatial rarefication of occurrence data and independent testing of model feature classes and regularization parameters. Spatial rarefication is a critical technique when points are not evenly distributed in space. The resulting model of the spatial structure of forest cover is based on the integration of the best models of each thematic class of different types of forest cover into a single cartographic layer. It is shown that under conditions of uneven and sparse distribution of points, it is possible to provide an average point matching level of 0.45 for formations and 0.29 for association groups. Herewith, the spatial structure and the ratio of the formation’s composition correspond to the official data of the forest inventory. An attempt is made to identify and evaluate the distribution of more detailed syntaxonomic units: association groups. The necessary requirements for improving the quality of the forest cover model of the study area for 2 hierarchical typological units of forest cover are formulated. These include the additional sampling in order to equalize their spatial density, as well as to achieve equality of samples based on stratification according to the resulting map.


Introduction
The mapping of the spatial structure of forest communities is an integral part of biodiversity research and environmental planning at the regional level. There are three basic constraints on forest data: First, the data on forest spatial structure must be up-to-date and be able to be regularly updated with a step of 2-5 years. Second, the data should equally describe the parameters of biodiversity and might potentially improve the robustness of the model [22]. In the current study, the SDMtoolbox (Spatial Distribution Modeling toolbox 2.4, Durham, NC/Manhattan, NY/Carbondale, IL, USA) is chosen as a modeling tool. The SDMtoolbox is a python-based ArcGIS 10.7 toolbox (Redlands, CA, USA) for spatial studies of ecology, evolution and genetics. The SMDtoolbox is chosen because it includes the basic MaxEnt (Maximum Enthropy, Manhattan, NY, USA) algorithm and a number of additional tools necessary to control the autocorrelation of spatial data [23]. Among other methods, MaxEnt is shown as an effective tool for non-linear interactions between response and predictor variables, and is robust to small sample sizes [21]. MaxEnt is used not only for the Species Distribution Model, but also for a wide range of natural phenomena, for instance, for tree pests monitoring (Salento Peninsula, Italy) [24], to create predictive risk maps for soil-transmitted helminth infections in Thailand [25], to study invasive species (Korea) [26] and critically endangered Alaotran gentle lemur (Madagascar) [27]. MaxEnt is even applied in mineral prospectivity analysis (Nanling, China) [28] or as the model for the prediction of landslide patterns (Arno Basin, Italy) [29]. The geology, soil, climate and vegetation spatial data are used in the above-mentioned cases along with remote sensing data. It can be assumed that the maximum entropy method together with the MaxEnt software is the universal geographical tool for spatial modeling and mapping.
It is also a developing practice to apply spatial modeling to such natural phenomena as forest formations. For example, the habitat suitability modeling was performed for 10 types of forest formations in Europe. The Random Forest model based on 1 km climate environmental variables and more than 6000 field data forest inventory plots were used in that study. The overall accuracy of the final map was 76% [30]. Finite mixture model was applied to a national forest inventory of Italy consisting of 6714 plots with a measure of abundance for 27 tree species, and the map of potential forest types was produced also based on 1 km climate data supplemented with some geological and soil data [31]. According to the study of Panamanian tree species, their distribution appears to be primarily determined by dispersal limitation, then by environmental heterogeneity. This study used a permutation-based regression model computed on distance matrices and a hierarchical clustering of the tree composition to construct a predictive map of forest types of the Panama Canal. Fifty-three sample plots describing the floristic composition along with climatic data, elevation, geologic formation and slope are used in the referred study [32]. The study of southern Atlantic Rainforest formations (Brazil) aimed to verify the existence of indicator species and identify relationships among distributions of tree species with environmental and spatial variables. The study was based on 21 sample plots and altitude and climatic variables [33]. The forest formations are often referred to in the aforementioned studies as Floristic Patterns or Species Composition Patterns. The satisfactory results are shown for application of MaxEnt in land-cover classification and land change analysis. For example, the study in Trentino-South Tyrol, Italy, developed land cover and difference maps between 1976 and 2001 based on multispectral data and topographic variables. MaxEnt applied to land cover classes can provide reliable data, especially when referring to classes with homogeneous texture properties and surface reflectance [34].
The purpose of this study is to assess the utility of modern approaches in spatial modeling of forest communities for the example of the Moscow Region, based on field data obtained outside the state forest inventory. The tasks of this study are dictated by the need to develop and adapt optimal methods for managing the array of field descriptions unevenly distributed in space and between syntaxonomic units, as well as to develop a probabilistic cartographic model of forest communities at the regional level, as an alternative to the generalized official data of the forest inventory.

Study Area
The Moscow Region is located in the central part of the East European (Russian) Plain: 35 • (Figure 1). The population is 20.4 million people (Moscow: 12.7 million people). The average population density is 2.29 people per square kilometer. The region has several important natural and phytogeographical boundaries. The significant border is the south edge of the natural range of the spruce forests in the broad-leaved spruce forests zone.
Region were recognized as a green zone, including a ban on industrial felling. However, up to the end of the 20th century, there was an increase in the pace of industrial development: the construction and operation of machine-building plants and related infrastructure, including enterprises of the energy complexes and oil refining complexes [35]. This has led to an increase in emissions of pollutants into the atmosphere and hydrosphere. At the beginning of the 21st century, the direction of influence changed to post-industrial. The growth of the population along with the construction of housing stock and transport infrastructure started in connection with the development of the financial sector of the economy [36,37]. In addition, deterioration in the volume of forestry activities aimed at maintaining the sustainability of forest plantations was recorded. The nature protection regime violations are as follows: unauthorized felling and household/industrial waste dumps. These violations result in large pest outbreaks, forest fires and degradation of species composition [17].

Design of the Study
1684 sample plots (including 1494 forest sample plots) were collected during 2010-2019 ( Figure  1). The sample plots' locations selection is based on the principle of representativeness for the main species prevailing in the forest stand, forest types, taking into account age groups, as well as the ecology of the habitat. The sampled vegetation communities are homogeneous in terms of the general floristic composition, the composition of the dominants of each layer, physiognomy (aspect, Due to the high soil fertility, the forest cover of the region has experienced an extensive anthropogenic impact (felling, plowing) for several centuries. At the beginning of the 20th century, and especially after the Second World War, there was a significant change in the direction of impact: the use of forestry practices on the site of the former arable land. In 1947, all forests of the Moscow Region were recognized as a green zone, including a ban on industrial felling. However, up to the end of the 20th century, there was an increase in the pace of industrial development: the construction and operation of machine-building plants and related infrastructure, including enterprises of the energy complexes and oil refining complexes [35]. This has led to an increase in emissions of pollutants into the atmosphere and hydrosphere. At the beginning of the 21st century, the direction of influence changed to post-industrial. The growth of the population along with the construction of housing stock and transport infrastructure started in connection with the development of the financial sector of the economy [36,37]. In addition, deterioration in the volume of forestry activities aimed at maintaining the sustainability of forest plantations was recorded. The nature protection regime violations are as follows: unauthorized felling and household/industrial waste dumps. These violations result in large pest outbreaks, forest fires and degradation of species composition [17].  species prevailing in the forest stand, forest types, taking into account age groups, as well as the ecology of the habitat. The sampled vegetation communities are homogeneous in terms of the general floristic composition, the composition of the dominants of each layer, physiognomy (aspect, community structure) and habitat conditions. According to the methodology, the sample plots were limited to 20 × 20 m and located at a distance of at least 200 m from each other. When localizing the sample plots, the representation from the main types of communities was taken into account. Lack of road infrastructure was a serious limitation in the uniform laying of test plots. The following properties of vegetation communities are recorded at each sample plot: • Composition and structure of the tree layer (projective crown cover, average height of mature trees and undergrowth).

•
Complete species composition of shrub, grass-dwarf shrub and moss layers, with an estimate of the cover in percent.

•
Species saturation of the plants of the ground layers, estimated as the average number of species per unit area (to assess the species diversity).
The previous study of Chernenkova [38] has shown that the ecological-phytocenotic classification is best for analysis of communities in this region. The techniques allows for great differentiation of communities across the region [38]. The use of units of ecological-phytocenotic classification in the rank of formation and association groups is explained by a number of reasons: (1) good correspondence of typological and mapped units, (2) compliance with Russian units of forest typology and (3) consideration of rare types of forest areas, as well as secondary derivative communities, which is important from an environmental point of view.
Formations are identified according to the dominant forest species. Each formation is represented by communities with a different combination of common dominants in the lower stories. The classifications of syntaxons at the level of association groups is carried out according to the dominant ecological and morphological groups of plants of subordinate stories [38]. The association groups are connected with the features of the herb-dwarf shrub layer, trophic conditions and moisture conditions. Comparison of the formations' composition and association groups allows one to assess the direction of successional development, the degree of anthropogenic impact and the stability of forest communities.
For mapping the spatial structure, the maximum entropy is chosen [39]. The choice of method is based on the nature of the field data collection. There are two main ways to collect data on the distribution of natural objects [40], which require the use of various modeling algorithms: One is collection of only occurrence points. With this method of analysis (it is called "presence only" or "presence/background" presence analysis), one can use GPS data, materials of collections and herbaria, publications describing the places of species registration. Another is separate collection of occurrence and absence points. This is called a "presence absence" analysis. The possibility of collecting correct information about the absence of a natural phenomenon in any territory is debatable and may turn out to be false. The maximum entropy method is designed specifically for processing "presence/background" data [9].
SDMtoolbox software includes additional tools, namely calibration of multiple models, their testing and final selection of the best model using the quality indicators [41]. An important part of the tool is the spatial jack-knifing (or geographically structured k-fold cross-validation) [42]. This tool tests evaluation performance of spatially segregated, spatially independent localities. The tool also allows for testing different combinations of model feature class types (linear, quadratic, hinge, product and threshold) and regularization multipliers to optimize MaxEnt model performance.
Landsat 8 and Landsat 5 spectral reflectances and spectral indices as well as Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) and morphometric variables are used as the environmental variables [43,44]. Palsar-2 25 m resolution images are also included: ortho and slope-corrected backscattering coefficient (horizontal-horizontal (HH) and horizontal-vertical (HV) polarization) for 2019 [45]. To manage the autocorrelation, the variables that have 95% correlation  Table A1. The Global Forest Watch dataset is utilized to prepare forest/non-forest masks as well as the loss year mask and water mask [46].
The modeling of forest cover is performed within two hierarchical levels: (1) forest formations and (2) association groups. Forest formations are aggregated syntaxons that are warranted by statistically sufficient and homogenous training samples. This makes forest formations more robust for spatial modeling. Association groups are more divisional, in that they have heterogenous training samples and thus they are more sensitive natural objects for modeling. The overall quality of modeling is evaluated by two confusion matrices.
On the first upper level, forest formations are modeled. Using the Global Forest Watch dataset, the territory is divided into two strata based on 30% forest cover threshold and forest loss/gain data: forest stratum and non-forest stratum. Raster masks are used to apply strata during mapping. Multiple approaches are tested for selection of bias layers [40][41][42]: no bias, and 10 and 25 km biases around occurrence points. The best results (highest area under the curve and 1-omission error) are obtained using forest layer as bias for forest formations and non-forest layer as bias for open lands and agriculture.
The systematic sampling approach that showed the best results along with bias layer, clustering, splitting and background restriction approaches is utilized [21]. To provide balance between sample sizes and sample equality, the sample size for each formation is systematically decimated to around 100 sample plots, because MaxEnt allows to include all non-linear feature class types only when the size is over 80 samples. To achieve 100 sample plots as well as to equalize sample sizes of formations and to reduce spatial autocorrelation, the following compromising approach is used. The spatial rarefication of occurrence data is applied.
Spatial jack-knifing is performed by three regions, five model feature class types are used (linear, quadratic, hinge, product and threshold) and three regularization multipliers (0.5, 1, 2). Final models of formations are integrated into one map by the method of highest position in ArcGIS. The confusion matrix is calculated between the final map and the initial sample plots.
On the second level, the modeling of association groups is performed. The number of sampling points varies significantly from 9 to 154: average 48.96, median 35, standard deviation 39.95. For this reason, no transformation and no spatial rarefication of the number of sample plots is applied. For every association group, the formation mask is applied. The confusion matrix is also calculated. Table 1 shows the hierarchical structure of formations (in columns) and association groups (in rows) of the Moscow Region, in accordance with the previously published results of the ecological-phytocenotic classification [43]. The cells indicate the number of field sample plots for the association groups.

Pre-Processing of Samples
Formation A: The composition of spruce forests is complex (combinations of spruce with birch, aspen, pine and broad-leaved species) and it is similar to the composition of the zonal primary coniferous broad-leaved forests. The proportion of silviculture is high (mainly monodominant spruce forests) [44]. The species composition of the subordinate stories (grass-dwarf shrub and moss stories) is represented by the full spectrum of transitions from boreal to nemoral types. A relatively small number of community types is noted in the composition of boreal spruce forests (small herb and small herb-green moss groups). While, subnemoral (small herb-broad herb) and nemoral (broad herd) groups of spruce forests have a higher coenotic diversity, due to the higher participation of other tree species, and the variety of combinations of dominant land cover species. Formation B: spruce-aspen/birch forests. This group of formations includes communities where spruce and birch are represented in equal proportions (with a small admixture of aspen). Spruce-small-leaved forests are interpreted by many researchers [47,48] as a short-term stage of spruce forests that form at the felling site as a result of both spontaneous succession and the development of spruce silviculture. The composition of the vegetation of the ground stories is close to that of spruce forests (Formation A).
Formations C and D: pine-spruce forests, and pine forests with spruce and birch, locally with linden, oak and hazel. Pine and pine-spruce forests on uplands are not completely indigenous communities and represent a successional stage in transition to mature forest communities. The absence of pine regeneration in automorphic habitats indicates a derivative origin of pine forests after fires and fellings, as well as in the composition of silvicultural forests. In one case, succession is accompanied by active recovery of spruce forests (Formation C). In the other case, in habitats with nutrient-rich soils, the broad-leaved succession is observed (Formation D). The broad-leaved species displace the pine and pine-spruce communities after a few decades. Since the ecological range of the pine is wide, the pine forests can be found on soils with different texture and moisture regime. Respectively, the pine forests' typological diversity is higher compared to spruce forests (Table 1).
Formation E: oak forests with linden, spruce and birch. This group characterizes broad-leaved forests with a predominance of oak in the first story of the stand, the participation of spruce in the first and second sub-stories and species of the nemoral group in the lower stories of communities.
Formation F: broad-leaved-spruce forests. Forests with a mixed composition of spruce and broad-leaved species (oak, linden, maple) and species of the nemoral group in the lower stories of communities. This is primary (indigenous) communities, which are usually replaced by linden or spruce forests during felling.
Formation G: linden with oak, locally with spruce and birch. Broad-leaved forests with a predominance of linden in the first or second sub-stories and nemoral species in the lower stories. Spruce is occasionally represented in the upper stories and in undergrowth. Oxalis and some other boreal species are involved in the grass story. These communities are primarily only on the slopes of ravines and river valleys and on the uplands in the central and northern parts of the region. Here, they are derivatives of coniferous-deciduous forests [49,50]. Indigenous linden nemoral grass forests are represented only in the southern part of the region.
Formation H: birch forests with spruce and aspen, locally with oak and linden, and Formation I: aspen forests with spruce and oak. The predominance of small-leaved birch and aspen forests in the region is associated with the formation of young forests in fellings. In recent decades, it is also associated with a regenerative succession on massively abandoned tillage. Most of the community types are secondary and can develop in a wide range of habitat conditions. The typological diversity of small-leaved forest communities of birch and aspen is associated with a wide ecological tolerance and the ability to grow on soils of different texture, moisture and nutrient richness. However this forest formation creates conditions for the restoration of conditionally primary communities.
Formations J and K: gray and black alder. These communities are more often referred to as primary forests, preferring either waterlogged or rich brook habitats. However, there is another point of view, according to which gray alder forests are considered derivative forests, developing on the site of broad-leaved spruce communities [51,52] (Table 1). Each forest formation is represented by one or multiple association groups: • Dwarf shrubs-small herb-green moss (DShG), Spatial rarefication for formations is performed in accordance with the methodology: A (10 km), B (10 km), C (1 km), D (5 km) and H (10 km). Two pairs of ecologically similar formations are merged: broadleaf formation (E and F) and alder formation (J and K). This made it possible to group the sample plots in accordance with the ecological similarity of syntaxa and achieve a relatively equal sample size: (82-112 sample plots in each sample).
The following non-forest habitats' field data were collected during field surveys: small leaf scrub (L), meadows (N), open marshy habitats (O) and agricultural fields (P). The following non-forest habitats are taken from the Global Forest Change dataset: cuts (M) and water objects (Q) [46]. Settlements (R) are taken from Openstreetmap (OSM) layers [53]. These habitats are mapped and merged with the final map through non-forest mask ( Table 2).

Modeling of Formations
Ten different methods of grouping formations were tested (Table 3). In methods 1 and 2, all sample plots are used and each formation is modeled individually, and in method 2, the DEM is excluded from environmental variables. In methods 3-6, formations E, F and G (oak, linden, broadleaf forests) are combined, in method 3, the DEM is excluded from the environmental variables and in method 4, only spectral reflectances of July and September are left. Method 5 combines the H and J formations (birch and gray alder forests). Methods 7-9 apply spatial rarefication of occurrence data. In method 7, only the formation points A (spruce) are rarefied. In methods 8 and 9, different rarefication distances are applied for formations A, B, C, D and H. Two pairs of ecologically similar formations are merged: E with F and J with K. In method 9, the DEM is excluded from environmental variables. The results of the eighth method of grouping of formations are closest to the ratio of tree species given in the Forest Plan (FP) [17]. The modeling is re-run with settings of the eighth method and addition of the Palsar-2 dataset (Figure 1). The results of modeling the dominant formations are in good agreement with the data of the Forest Plan: birch 30.1% (39.6% in the FP), spruce 25.1% (24.4% FP), pine 18.5% (20.7% FP), aspen 5.1% (8.4%) and oak with broad-leaved (E + F) 2.9% (1.78% FP). The proportions of linden (G) and alder (J + K) are overestimated by 7.3 and 3.1 times respectively, compared to the data of the Forest Plan (Figure 2). Several possible reasons for this discrepancy are suggested. First, there is a lack of field data. Naturally, in the absence of preliminary stratification, such formations will rarely be encountered during field routes. At the same time, despite the fact that the number of points of the rare formations has been conditionally brought to the level of 80-100, they nevertheless still have a significant role of autocorrelation: often these points are located in spatial clusters, which reduces the quality of models. However, another factor is also important. In the course of field work, it was repeatedly noticed that the forests indicated on the forest inventory maps as birch are in fact alder forests. According to the confusion matrix (Table 4), the classification accuracy is 0.46. The best classification quality is for alder forests. This is noteworthy in the context of the above difference with the official data of the Forest Plan. A low level of matching with field data is found for the formations of aspen (I) 0.   According to the confusion matrix (Table 4), the classification accuracy is 0.46. The best classification quality is for alder forests. This is noteworthy in the context of the above difference with the official data of the Forest Plan. A low level of matching with field data is found for the formations of aspen (I) 0.

Modeling of Association Groups
Despite the uneven spatial distribution of points, modeling of association groups is performed without spatial rarefication. Table 5 shows the results of modeling of association groups. The average level of points matching between all association groups is 0.29 and it ranges from 0.03 to 0.69. On average, a low percentage of points matching is typical for coniferous groups of associations (A-D) and broad-leaved conifers (F). Oak-spruce (E), linden (G), birch (H) and aspen (I) have an average level of recognition. The best quality is for gray alder (J) and black alder (K) formations.
According to the assessment of the spatial distribution of the identified association groups of forest communities (Table 5), the largest area (12.5%) is occupied by communities of derivative birch and aspen forests (Formations H and I) with a predominance of the mesotrophic and hydromorphic series (Bh, MhBh, Gm). The proportion (4.7%) of birch (H) grass-marsh forests (Gm), which are distinguished by strong recreational disturbance, is also high. Within the pine formation (D), the maximum proportion (about 5%) is occupied by communities (DShG and Sh), which tend to succession towards boreal spruce forests. The same pattern is observed in the formation of mixed spruce-small-leaved forests (B), occupying an area (4.72%) where their successional dynamics are also directed towards the restoration of spruce communities of groups (DShG and Sh). Another part of these communities of formation B (3.3%), in terms of the composition of the ground layer (ShBh and Bh), has a tendency to succession towards broad-leaved communities with a nemoral composition of the ground layer. It is obvious that this group is of artificial origin. The area of communities within the spruce forests (A) is small and varies in the range of 0.47-1.37%. Table 6 shows the proportion of the area occupied by non-forest habitats.

Discussion
A wide range of statistical models exist to model the spatial structure of natural features and phenomena: Resource Selection Function, Generalized Linear Models, Artificial Neural Networks, Maximum Entropy and Classification and Regression Trees [21]. Formerly for the Moscow Region, it was shown that Linear Discriminant Analysis demonstrates satisfactory results for modeling the spatial structure of forests based on field data and Landsat 8 spectral reflectance supplemented with digital elevation model and their derivative parameters [22]. However, linear functions of multispectral data demonstrate limited capabilities, which is demonstrated in our study. The average proportion of matching between field data and model was 52%. It varied from 20% to 100% for 38 association groups based on 1025 field sample plots. However, the number of sample plots per each association group varied from 4 to 114, which is unlikely in terms of sample size [54].
Previous works have shown the advantages of ecological-phytocenotic classification over ecological-floristic for the purpose of mapping. The relative quality of the discriminant analysis of the identified syntaxa within the ecological-floristic classification demonstrated a lower accuracy of the ecological-floristic classification (69.7%) compared to the ecological-phytocenotic classification (78.6%) [55].
The assessment of the spatial diversity of forest cover is connected with a number of limitations. Within the current study, the limitations may be generalized into three groups. The group of natural factors makes the most irregular and heterogeneous input into classification uncertainty. Foremost, the presence of multidirectional processes with degression dynamics (recreational impact, road and construction infrastructure) and restoration dynamics (tillage abandoning, forest silviculture). This limitation significantly disrupts the natural composition of coenotic types. The study area is located in the zone of Eastern European deciduous-coniferous forests, characterized by a mixed polydominant composition, which is difficult to analyze the species composition of communities and their classification, which has also been noted by other authors [56,57]. With regard to deciduous forests, an increase in their proportion due to climate warming, as well as the warming effect of the megalopolis, cannot be ruled out [58]. The polydominance of the tree layer [49], transitional succession status of most derivative forests [59], anthropogenic disturbance, as well as a high proportion of forest silviculture in the region [60], together with the anthropogenic impact mentioned above, make it one of the most complicated regional study areas.
Another group of classification uncertainty factors is generated by the quality and properties of environmental variables, especially Landsat spectral reflectance. Even excluding uncertainties of nadir angle, radiometric correction and atmospheric transparency, each of them varies by 5% to 7.5% [61], and there is still one important issue: the area of the Landsat pixel (0.09 hectares) slightly exceeds the area of the sample plot (0.04 hectares). In terms of statistics, it means that within each sample plot, one cannot evaluate statistical parameters of each sample of pixels and reflectance data are used with all potential extremes and outbreaks. Simply speaking, the Landsat dataset is too coarse for modeling typical 20 × 20 m sample plots [62]. A few solutions might be discussed in this context. One obvious solution is use of higher resolution imagery-Sentinel-2. Not to be overlooked is filtering of spectral reflectance, i.e., median filter [63].
The third group of uncertainty factors is considered in the Introduction Section, and it concerns the uneven spatial distribution of field data, typical for the regions of the Russian Federation, which is a critical factor affecting the quality of the models. This group might be characterized as of human and organizational origin. The method applied in the current study made it possible to obtain a map of the spatial structure of the formations corresponding to the official data for most of the species (birch, spruce, pine, oak, broad-leaved species and linden) and slightly overestimated results for some rarely distributed species (alder and aspen). The applied set of methods allowed to reach overall accuracy of 0.46 for forest formations.
Modeling more detailed syntaxa of forest cover, association groups, for which spatial rarefication of points is not used, emphasizes the negative contribution of uneven field data. This is especially important for coniferous and broad-leaved coniferous formations. However, the spatial pattern is quite plausible, the results are consistent with the previously developed forest cover models for part of the Moscow Region based on discriminant analysis [43,44], and also with a map of the vegetation cover of the Moscow Region [18].
An attempt at large-scale mapping is promising for assessing biodiversity and forest dynamics, but it has limitations in the area of study with the characteristic physical and geographical diversity of the territory. The overall technology holds promise, but still, the uncertainty of classification is rather low and one shall look for utilizing the higher spatial resolution datasets along with filtering approaches. The previously performed work in the southwestern part of the Moscow Region demonstrated higher quality of the cartographic model (78.6%) of the distribution of 15 types of forest communities [55]. Thus, for large regions with a complex natural structure and anthropogenic history, it might be useful to perform modeling within the individual landscape structures.
The results obtained underline the need to use the resulting map as a stratification matrix and to carry out additional field research, systematic and optimally justified. Additional field research should be aimed at achieving the following objectives: • Creation of the set of field descriptions, evenly distributed in space and taking into account rare and remote habitats.

•
Bringing the minimum number of descriptions of association groups to at least 50 (additional 494 descriptions), and in the long term, to 80 (1240 additional descriptions).

Conclusions
The use of MaxEnt nonlinear modeling together with additional tools (geographically structured spatial jack-knifing, spatial rarefication of occurrence data and independent testing of model feature classes and regularization parameters) can be used to manage the problem of uneven distribution of field data and to attempt to create a probabilistic cartographic model of forest formations at the regional level. The results of our modeling correspond well to the official data of forest inventory despite the high level of modeling uncertainty.
The main limitations of identifying and assessing the spatial distribution of types of forest communities at a more detailed typological level, in the rank of groups of associations, additionally to those mentioned above, including a series of studies at the sub-regional level within territorial units of natural zoning, were formulated. The need to utilize higher spatial resolution datasets along with filtering was emphasized.
The resulting cartographic model of the groups of associations can be used to stratify the study area and plan the optimal number and placement of field routes necessary for the final statistically valid model of forest communities.
Author Contributions: I.K. designed the study and performed the modeling; I.K. and T.C. did research and the data-analysis; T.C. performed ecological-phytocenotic classification of forest communities. All authors have read and agreed to the published version of the manuscript.

Funding:
The Russian Science Foundation (project no. 18-17-00129) supported this study. This study is also conducted in the framework of the Institute of Geography RAS (project no. 0148-2019-0007) in terms of studying the composition of forest communities.

Conflicts of Interest:
The authors declare no conflict of interest.