Automatic rapid calculation method for ecological system function service evaluation

Abstract

The invention discloses an automatic rapid calculation method for ecosystem function service evaluation, and the method comprises the steps: presetting a spatial basic grid and a coordinate system covering a research region, obtaining meteorological text data, remote sensing data and other raster data, carrying out all kinds of preprocessing on the three kinds of data, and unifying the three kinds of data into the same data specification; interpolating or resampling the preprocessed meteorological text data, remote sensing data and other raster data, and unifying spatial numerical values into standard grid data; according to a timestamp carried in a standard grid data name, executing year-month addition or average calculation, unifying sampling time, deleting all spatial information, and converting into matrix data; for the matrix data, the water conservation amount, the water and soil conservation amount and the wind prevention and sand fixation amount are calculated according to the water conservation model, the water and soil conservation model and the sand prevention and sand fixation model respectively, and result matrixes are output respectively. According to the method, the operation process is simplified, manual operation links are reduced, the error probability is effectively reduced, and the accuracy and consistency of results are ensured.

Description

Technical Field

[0001] The present invention relates to the field of computer information technology, and particularly to an automated rapid calculation method for evaluating ecosystem function services. Background Art

[0002] Ecosystems, through their natural structures and processes, provide key services such as water conservation, soil and water conservation, and windbreak and sand fixation. These services are crucial for maintaining the ecological balance of the Earth and the sustainable development of human society. With the exacerbation of global environmental problems such as climate change, land degradation, and biodiversity loss, the protection and management of these services have become the focus of attention of governments and international organizations.

[0003] In recent years, with the progress of remote sensing technology and ecological theory, a series of large-scale, long-term climate remote sensing data, as well as various supporting ecological models, have emerged. This has not only improved the accuracy and coverage of assessments but also enhanced the capture and understanding of the spatio-temporal changes in ecosystem service functions. However, using these data and technologies for ecosystem service function assessment still has high difficulties, mainly reflected in the following three aspects: (1) Data format conversion and unification. Ecosystem service assessment requires integrating multiple data sources, such as meteorological data in NetCDF (.nc) format, MODIS (.hdf) remote sensing data, and various raster format geographic information data. These data usually have different formats, resolutions, projection methods, and sampling times, so complex format conversion and data unification are required. This process not only consumes time and effort but also requires a variety of professional software tools and technical knowledge to ensure data consistency and availability; (2) Computational complexity. The assessment process involves multiple models and formulas, such as the water balance equation, vegetation cover factor model, soil erodibility equation, modified wind erosion equation, etc. These models usually need to consider multiple input variables, look up tables for calculation, and process large-scale spatial data. This makes the calculation process very complex and error-prone; (3) Long process. The assessment process is cumbersome, including multiple steps such as data collection, preprocessing, format alignment, model calculation, etc. Each step requires a large amount of time for technical operations and professional judgment. Especially in the model calculation stage, many algorithm formulas are nested up and down rather than existing independently. This complex process not only requires the team to have a high level of organizational coordination ability and professional knowledge but also increases the human burden. Once there is a personal technical error, it will directly affect the accuracy of the overall result.

[0004] In order to improve the efficiency of ecosystem function service assessment, how to unify the calculation format, reduce the computational complexity, and shorten the process length has become a difficult problem to be solved urgently in the assessment. Summary of the Invention

[0005] The object of the present invention is to provide an automated and rapid calculation method for evaluating ecosystem function services. By simplifying the operation process and reducing the manual operation links, the probability of errors is effectively reduced, ensuring the accuracy and consistency of the results.

[0006] To achieve the above object, the present invention provides the following solutions:

[0007] An automated and rapid calculation method for evaluating ecosystem function services, comprising the following steps:

[0008] S1. Preset the spatial basic grid and coordinate system covering the study area, obtain meteorological text data, remote sensing data, and other raster data, and perform various preprocessings on the above three types of data to unify them into the same data specification;

[0009] S2. Interpolate or resample the preprocessed meteorological text data, remote sensing data, and other raster data to unify the spatial values into standard grid data;

[0010] S3. According to the timestamps carried in the standard grid data names, perform monthly or average calculations, unify the sampling time, delete all spatial information, and convert it into matrix data;

[0011] S4. For the matrix data, calculate the water conservation volume, soil and water conservation volume, and wind prevention and sand fixation volume according to the models of water conservation, soil and water conservation, and wind prevention and sand fixation respectively, and output the result matrices respectively;

[0012] S5. Add the spatial information of the preset coordinates to the result matrix and convert it into GeoTif format for output.

[0013] Preferably, in S2, when interpolating or resampling the preprocessed meteorological text data, specifically including:

[0014] Obtain the meteorological text data in NetCDF format, convert the meteorological text data in NetCDF format into JSON format, extract meteorological parameters, and structurally store them as NCDS = [x, y, time, value], where x and y are the longitude and latitude coordinates of the data points, time is the timestamp, and value is the climate value of the data points; use the datetime module to unify the timestamp into daily granularity to generate a daily climate text data file NCD; spatialize the NCD file into a set of coordinate points with values, and use the inverse distance weighted interpolation method to interpolate the meteorological point data to the spatial basic grid to generate a daily climate grid data set.

[0015] Preferably, in S2, when interpolating or resampling the preprocessed remote sensing data, specifically including:

[0016] Obtain remote sensing data in MODIS format, parse the file names of the remote sensing data, and extract the product code, shooting time, and satellite orbit row and column numbers; according to the satellite orbit row and column numbers and time intervals of the study area, automatically combine the names of the remote sensing data to generate a name dictionary; batch read the remote sensing data according to the name dictionary, splice them according to the satellite orbit row and column numbers, and generate a daily remote sensing image set; use the nearest neighbor interpolation method to resample the daily remote sensing image set to the spatial basic grid to generate a regular grid data set.

[0017] Preferably, in S2, interpolate or resample the preprocessed other raster data, specifically including:

[0018] Obtain other raster data in tif format, and resample the tif format data to the spatial basic grid using the nearest neighbor interpolation method to generate a regular grid data set.

[0019] Preferably, in S3, according to the timestamps carried in the standard grid data names, perform monthly addition or average calculation, specifically including:

[0020] The standard grid data includes daily climate grid data sets, regular raster data sets, and regular grid data sets. Convert the daily climate grid data sets, regular raster data sets, and regular grid data sets into Numpy matrices and remove the spatial information; according to the timestamps in the data names, use the datetime module to calculate the monthly average, monthly sum, annual average, or annual sum to generate a matrix data set with a unified time granularity.

[0021] Preferably, in S4, for the matrix data, perform relevant calculations according to the water conservation model, including:

[0022] The surface runoff factor is calculated as follows:

[0023] R = P × α

[0024] Where, R represents the annual surface runoff, P represents the annual precipitation, and α is the average surface runoff coefficient; the water conservation amount is calculated as follows:

[0025] TQ = (P - R - ET) × A × 10 3

[0026] Where, TQ is the water conservation amount, P is the annual precipitation, ET is the actual evapotranspiration, R is the surface runoff, and A is the grid area.

[0027] Preferably, in S4, for the matrix data, perform relevant calculations according to the soil and water conservation model, including:

[0028] The rainfall erosivity is calculated as follows:

[0029]

[0030] In the formula, R is the average annual rainfall erosivity, P i is the total monthly precipitation, and P is the total annual precipitation; i is the month of the precipitation used, i = 1, 2,..., m; the soil erodibility factor is calculated as follows:

[0031]

[0032] K = (-0.01383 + 0.51575K EPIC ) × 0.1317

[0033] In the formula, K EPIC represents the soil erodibility factor before correction, and K represents the soil erodibility factor after correction. m c , m silt , m s and orgC are the proportions of clay, silt, sand, and organic carbon respectively; the soil and water conservation is calculated as follows:

[0034] A c = A p - A r = R × K × L × S × (1 - C)

[0035] In the formula, A c is the soil and water conservation amount, A p is the potential soil erosion amount, A r is the actual soil erosion amount, R is the rainfall erosion factor, K is the soil erodibility factor; L is the slope length factor, S represents the slope factor, and C is the vegetation cover factor.

[0036] The present invention also provides a non-transitory computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, it implements an automated rapid calculation method for evaluating ecosystem function services as described in any one of the above.

[0037] According to the specific embodiments provided by the present invention, the following technical effects are disclosed:

[0038] (1) The present invention proposes an automated rapid calculation method for evaluating ecosystem function services, which supports the automatic parsing, format conversion, and spatial alignment of NetCDF meteorological data, MODIS remote sensing data, and TIF raster data, and solves the problems of data format incompatibility and projection inconsistency in traditional methods; uses inverse distance weighted interpolation (IDW) and nearest neighbor interpolation to ensure the unification of data with different resolutions to a standard grid; uses standardized matrix calculation to unify all data formats, avoiding manual adjustment, and improving calculation accuracy and consistency;

[0039] (2) The present invention incorporates three core models for water conservation, soil and water conservation, and windbreak and sand fixation. Among them, the water conservation model calculates the water conservation volume of each grid automatically through built-in formulas based on data such as annual precipitation, annual evapotranspiration, surface runoff, and ecosystem types. The soil and water conservation model calculates variables such as rainfall erosivity and soil erodibility factor, combines factors such as land cover and vegetation cover, and automatically executes corresponding mathematical formulas to calculate soil erosion. The windbreak and sand fixation model combines climate factors (such as wind speed and wind force factor) and surface soil characteristics (such as soil moisture and vegetation cover), and automatically evaluates the wind and sand control effect using the wind erosion calculation formula. By combining these three core ecological models, the present invention realizes the fully automated calculation of ecosystem function service evaluation. All calculations in the model are based on a standardized matrix, and the data connection between several complex scientific formulas and models in the model has been set up, eliminating the need for manual step-by-step input of formulas and avoiding data loss or calculation errors in the middle. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required in the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0041] Figure 1 It is a design framework diagram of an automated rapid calculation method for ecosystem function service evaluation provided by the present invention;

[0042] Figure 2 It is a spatial basic grid map covering the research area in an embodiment of the present invention;

[0043] Figure 3 It is a MODIS preprocessing interface diagram in the software to which the rapid calculation method in an embodiment of the present invention is applied;

[0044] Figure 4 It is an NC preprocessing interface diagram in the software to which the rapid calculation method in an embodiment of the present invention is applied;

[0045] Figure 5 It is a water conservation calculation interface diagram in the software to which the rapid calculation method in an embodiment of the present invention is applied;

[0046] Figure 6 It is a soil and water conservation calculation interface diagram in the software to which the rapid calculation method in an embodiment of the present invention is applied;

[0047] Figure 7 It is a windbreak and sand fixation calculation interface diagram in the software to which the rapid calculation method in an embodiment of the present invention is applied;

[0048] Figure 8 This is a partial calculation result diagram in an embodiment of the present invention;

[0049] Among them, (a) is the preprocessing result of MODIS remote sensing data, (b) is the preprocessing result of NC meteorological data, (c) is the wind prevention and sand fixation result, (d) is the soil and water conservation result surface, and (e) is the water source conservation result. Specific implementation manners

[0050] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0051] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below in conjunction with the accompanying drawings and specific implementation manners.

[0052] As Figure 1 shown, an automated rapid calculation method for evaluating ecosystem function services provided by the present invention includes the following steps:

[0053] S1. Preset the spatial basic grid and coordinate system covering the research area, obtain meteorological text data, remote sensing data, and other raster data, and perform various preprocessings on the above three types of data to unify them into the same data specification;

[0054] S2. Interpolate or resample the preprocessed meteorological text data, remote sensing data, and other raster data to unify the spatial values into standard grid data;

[0055] S3. According to the timestamp carried in the standard grid data name, perform monthly or average calculation, unify the sampling time, and delete all spatial information, and convert it into matrix data;

[0056] S4. For the matrix data, calculate the water source conservation amount, soil and water conservation amount, and wind prevention and sand fixation amount respectively according to the models of water source conservation, soil and water conservation, and sand prevention and fixation, and output the result matrices respectively;

[0057] S5. Add the spatial information of the preset coordinates to the result matrix and convert it into GeoTif format for output.

[0058] Specifically, 1. The basic principle of the rapid calculation method design is as follows:

[0059] The basic principle of the design of the automated rapid calculation program for ecosystem function service assessment is as follows: preset the basic spatial grid and coordinate system, perform various format conversions and interpolations on meteorological text data NetCDF (.nc), remote sensing data MODIS (.hdf), and other tif data, put the values into the regular grid, uniformly convert them into a Numpy matrix, then use the datetime module to unify the sampling time, execute relevant model calculations, and after obtaining the final numerical results, spatially output them as GeoTif. The specific implementation principle is described as follows:

[0060] 1) Preset the basic spatial grid that can cover the study area, set the grid size (default is 250m) and coordinate system (CGCS2000, EPSG:4490), and set the default value of the grid as np.nan (this format will only return nan for any mathematical operation unless overwritten. This value can keep the administrative region boundary unchanged during the calculation because it is nan outside the boundary), that is

[0061] 2) Read all meteorological text data NetCDF (.nc), convert them into JSON format at one time, select the corresponding meteorological data (such as air temperature temp, rainfall pre, etc.) from them, and change the unordered text data into regular text NCDS = [x, y, time, value], where xy are the longitude and latitude coordinates of the data point, time is the timestamp, and value is the climate value of the data point;

[0062] 3) The time in the NCDS data is a continuous and irregular timestamp of "year:month:day hour:minute:second". Use the datetime module to truncate the long text NCDS into an NCD file for each day;

[0063] 4) Spatialize the series of NCD files into a set of coordinate points with values, and finally use the inverse distance weighted interpolation method to interpolate the point set into Net_empty. After interpolation, obtain a daily climate grid data set {Net_NCD};

[0064] 5) Use the (Operating system, OS) to read the names of MODIS (.hdf) data, and parse these file names (for example: MOD16A2GF.A2023001.h26v04.061.2024024042743.hdf) into the product code (MOD16A2GF) + remote sensing shooting time (2023001, 001 represents the 1st day of 2023, which also needs to be converted into a standard date like January 1, 2023) + satellite orbit row and column number (26 / 4);

[0065] 6) Automatically combine the MODIS data names according to the satellite orbit row and column numbers (e.g., [26,04], [26,05], [27,04], [27,05]) and time interval (1 day) in the research area to obtain a name dictionary:

[0066] MODIS_Dict = {day1: [MdName1, MdName1..], day2: [MdName3, MdName4..]...}, where day represents the remote sensing shooting time at day intervals and serves as the key of the dictionary, and MdName represents the MODIS (.hdf) data names included on this day and serves as the value of the dictionary;

[0067] 7) According to the keys in MODIS_Dict, batch obtain the MODIS (.hdf) data names that meet the composite conditions, find the corresponding files in the folder, read them into multiple raster image data, and automatically splice them according to the orbit row and column numbers to form a remote sensing image set {MODIS_TIFS} with one image per day;

[0068] 8) Batch resample the raster images in the remote sensing image set {MODIS_TIFS} to Net_empty using the nearest neighbor method to form a regular grid data set {MODIS_Net};

[0069] 9) Similarly, resample other tif data to Net_empty using the nearest neighbor method to obtain a regular grid {Others_Net} of other input data;

[0070] 10) Remove the projected coordinate system from all of {Net_NCD}, {MODIS_Net}, and {Others_Net}, and convert them into a series of Numpy matrices, Nums = {Net_name, MODIS_name...}, where Net_name, etc. are the original image names of these matrices, such as rainfall data on a certain day, soil sand particle proportion data in a certain year, vegetation coverage in a certain year, etc.;

[0071] 11) For the matrices in Nums, according to the four calculation methods of time, monthly average, monthly sum, annual average, and annual sum, use the datetime module to read the timestamps on the data names to obtain the data of the current month or year (the reason for using the datetime module is that the acquisition time of remote sensing images is irregular and cannot be directly divided by 31 or 365), and then perform addition or average calculations to obtain a unified time-granularity Nums_time data set;

[0072] 12) For the matrix in Nums_time, perform calculations such as the water balance equation, vegetation cover factor model, and soil erodibility equation. After a series of previous treatments, the Numpy matrix has corresponding time, space, rows, and columns, and can directly participate in mathematical formula calculations without other settings. The calculation result is also a Numpy matrix;

[0073] 13) For the output Numpy matrix result, attach the preset spatial coordinate system (CGCS2000, EPSG:4490) and convert it to GeoTIF for output.

[0074] 2. Scientific algorithms included in the program:

[0075] The calculations of water conservation, soil and water conservation, and sand fixation are also important parts of this method. Continuing from the previous description of fast calculations and implementation steps, this section will elaborate on the theoretical formulas of these scientific calculations and how to use matrices in a unified format to implement the calculations of these formulas in the program.

[0076] 1) Calculation formula for water conservation - surface runoff factor

[0077] R = P × α

[0078] In the formula, R represents the annual surface runoff (mm), P represents the annual precipitation (mm), and α is the average surface runoff coefficient. The value of α needs to be obtained according to the attribute information of the ecosystem type and by referring to Table 1. In the calculation of this part, the program reads the annual surface runoff matrix, annual precipitation matrix, and ecosystem type matrix (using numbers to replace the Chinese codes of the ecosystem types) simultaneously. The program then generates an empty matrix and superimposes the above three matrices. Each grid is traversed sequentially from top to bottom and from left to right. For each grid, the program first queries the corresponding Chinese code type through the ecosystem type matrix, then queries the corresponding average surface runoff coefficient α value through the mean of the surface runoff coefficients of various ecosystems in Table 1, and finally uses this α value and the values of the annual surface runoff matrix R and annual precipitation matrix P at the corresponding position to calculate the surface runoff value R according to the formula and assign it to the corresponding position in the empty matrix.

[0079] Table 1

[0080]

[0081] 2) Water conservation - ecosystem area factor

[0082] Calculate the grid area (m 2 ) according to the ecosystem type data, which is the default grid area of 250m * 250m.

[0083] 3) Water conservation - calculation of water conservation volume

[0084] TQ=(PR-ET)×A×10 3

[0085] Where TQ is the water conservation capacity, P is annual precipitation (mm), ET is actual evapotranspiration (mm), R is surface runoff (mm), and A is the grid area (km²). In this calculation, the grid area is set to 250 m × 250 m by default. The program reads the annual precipitation matrix, the actual evapotranspiration matrix, and the surface runoff matrix and directly applies the formula to calculate the water conservation capacity matrix TQ.

[0086] 4) Soil and Water Conservation - Calculation Formula for Rainfall Erosivity

[0087]

[0088] Where R is the average rainfall erosivity over many years (MJ·mm / hm2·h·a); P i is the monthly total precipitation, P is the annual total precipitation, and i is the precipitation month used, i = 1, 2, ..., m. In this calculation, the program reads 12 monthly rainfall matrices and one annual precipitation matrix at a time, then uses the above formula to calculate the multi-year average rainfall erosivity matrix R.

[0089] 5) Soil and Water Conservation / Wind and Sand Control - Soil Erodibility Factor Calculation Formula

[0090]

[0091] K=(-0.01383+0.51575K EPIC )×0.1317

[0092] Where K EPIC represents the soil erodibility factor before correction, K represents the soil erodibility factor after correction, m c 、m silt 、m s and orgC are the proportions of clay, silt, sand, and organic carbon, respectively (between 0 and 1). In this calculation, the program first checks the numerical ranges of the soil clay fraction matrix, soil silt fraction matrix, soil sand fraction matrix, and soil organic carbon fraction matrix. Once they are within the interval [0, 1], the above formula is applied to obtain the modified soil erodibility factor matrix K.

[0093] 6) Calculation formula for soil and water conservation-vegetation cover factor

[0094] C 旱 =0.221-0.595logc1

[0095] Where C 旱C is the vegetation cover factor for dryland (ecosystem type), and c1 is the vegetation coverage in decimal form. For paddy fields, wetlands, and urban areas as ecosystem types, C is assigned 0, 0, and 0.01 respectively. For the remaining ecosystem types, after superimposing the vegetation coverage and ecosystem type data and calculating, the assignment is made according to Table 2. In this part of the calculation, the program will simultaneously read the vegetation coverage matrix and the ecosystem type matrix, superimpose the two matrices, and traverse each grid from top to bottom and from left to right in sequence. For each grid, the program first queries the corresponding Chinese coding type through the ecosystem type matrix. If the coding is paddy field, wetland, or urban area, the C value is assigned 0, 0, and 0.01 respectively. If the coding is forest, shrubbery, grassland, arbor garden, or shrub garden, the vegetation coverage at the corresponding grid position is read, and the C value of this grid is obtained according to the assignment of the vegetation cover factor for different ecosystem types in Table 2. Finally, the C values of all grids are calculated using the above formula 旱 to form a vegetation cover factor matrix.

[0096] Table 2

[0097]

[0098] 7) Soil and water conservation - Soil and water conservation calculation formula

[0099] A c = A p - A r = R × K × L × S × (1 - C)

[0100] In the formula, A c is the soil and water conservation amount (t / hm2·a); A p is the potential soil erosion amount; A r is the actual soil erosion amount; R is the rainfall erosion factor (MJ·mm / hm2·h·a); K is the soil erodibility factor; L and S are topographic factors, L is the slope length factor, S represents the slope gradient factor, and C is the vegetation cover factor. The rainfall erosion factor, soil erodibility factor, slope length factor, slope gradient factor, and vegetation cover factor are all processed matrices and can be directly used in the formula to calculate and obtain the soil and water conservation amount matrix A c .

[0101] 8) Windbreak and sand fixation - Climate factor calculation formula

[0102]

[0103] In the formula, WF is the climate factor (kg / m); Wf is the wind force factor, μ i is the monthly average wind speed (m / s) in the i-th month, N dis the number of days with wind speed > 5 m / s in the i-th month; ρ is the air density (1.29 kg / m3); g is the acceleration due to gravity (9.80 m / s 2 ); SW is the annual average soil moisture, dimensionless, ET is the total annual actual evapotranspiration (mm / yr), PET is the total annual potential evapotranspiration (mm / yr); SD is the snow cover factor, a dimensionless constant, N SC is the number of times the annual SC (snow cover data observed by Modis every 8 days) equals 200, and N is the number of observations (365 days divided by observations every 8 days). In this part of the calculation, the program reads the 12-month average wind speed matrix and the matrix of the number of days with wind speed > 5 m / s in 12 months at one time (statistical daily wind speed data to obtain the number of days with wind speed > 5 m / s for each grid), and the rest of the data are sorted matrices. As for the snow cover data (observed once every 8 days), the number of grids equal to 200 in this matrix is counted as the constant value N SC , imported into the SD formula for calculation, and finally the climate factor matrix WF is obtained.

[0104] 9) Calculation formula for windbreak and sand fixation - soil crust factor

[0105]

[0106] In the formula, cl is the soil clay content (%), and OM is the soil organic matter content (%). In this part of the calculation, the program will first check the value ranges of the soil clay proportion matrix and the soil organic matter content matrix, and after determining that they are within the interval [0,1], execute the above formula to obtain the soil crust factor matrix SCF.

[0107] 10) Calculation formula for windbreak and sand fixation - vegetation coverage factor

[0108] Vegetation is divided into six vegetation types: forest, shrub, grassland, farmland, bare land, and desert. The vegetation coverage factors of each type are calculated according to different coefficients.

[0109]

[0110] In the formula, SC is the vegetation coverage (%) and a i is the coefficient for different vegetation types, which are: forest 0.1535, grassland 0.1151, shrub 0.0921, bare land 0.0768, sandy land 0.0658, farmland 0.0438. The program will read the vegetation coverage matrix and the ecosystem type matrix at the same time, and superimpose the two matrices. Each grid is traversed sequentially from top to bottom and from left to right. For each grid, the program first queries the corresponding Chinese coding type through the ecosystem type matrix and obtains the coefficient value a of the specified type by referring to the coefficient conversion table i, combined with the vegetation coverage SC at the corresponding grid position, calculate the vegetation coverage factor for each grid to obtain the vegetation coverage factor matrix C.

[0111] 11) Calculation formula for windbreak and sand fixation - surface roughness factor

[0112]

[0113] In the formula, K r is the ridge roughness, cm; C rr is the random roughness factor, taking 0, unit cm; ΔH is the altitude difference within the distance L, and L is the terrain undulation parameter, L = 250m. In this part of the calculation, as long as the altitude difference matrix is imported into the program and combined with the above formula, the surface roughness factor matrix K' can be obtained.

[0114] 12) Calculation formula for windbreak and sand fixation - windbreak and sand fixation

[0115] The modified wind erosion equation is used to calculate the windbreak and sand fixation amount, and the formula is as follows:

[0116] S = 150.71·(WF × EF × SCF × K' × C) -0.3711

[0117] Q MAX = 109.8[WF × EF × SCF × K' × C]

[0118]

[0119] Q MAX潜 = 109.8[WF × EF × SCF × K']

[0120] S 潜 = 150.71(WF × EF × SCF × K') -0.3711

[0121]

[0122] SR = S L潜 - S L

[0123] In the formula, SR is the windbreak and sand fixation amount (t / km2·a); S L潜 is the potential wind erosion amount (t / km2·a); S L is the actual wind erosion amount (t / km2·a); Q MAX is the maximum transfer amount (kg / m), S 潜 is the potential sand amount (kg / m), Q MAX潜is the potential maximum transfer amount (kg / m); Z is the maximum wind erosion occurrence distance (= 250 m); WF is the climate factor (kg / m); K' is the surface roughness factor; EF is the soil erodibility factor; SCF is the soil crust factor; C is the vegetation cover factor. In this part of the calculation, all input data are preprocessed matrices or result matrices calculated previously. According to the above formula, the wind and sand fixation amount result matrix SR can be obtained through calculation.

[0124] 3. Method efficiency:

[0125] To test the efficiency of this method, a comparison was made between the fast calculation program and the traditional multi-software combination calculation method. The experimental data used were the meteorological data NetCDF (.nc) in Beijing in 2023, remote sensing data MODIS (.hdf), and other tif data (annual potential evapotranspiration, soil silt content data, etc.), a total of 28 input data, with a total data volume of 2.8G, covering the basic spatial grid of the study area such as Figure 2 , the grid size is 250m, and the coordinate system is CGCS2000.

[0126] The traditional multi-software combination calculation method is extremely cumbersome. First, it is necessary to use Panoply to preprocess the meteorological data in NetCDF format (.nc), convert it into multiple tif image files one by one, and further crop and project them in ArcGIS. For the remote sensing data MODIS (.hdf), it must be preprocessed through the MRT software, then converted into tif format, and then these images are divided into pairs of four according to the row and column numbers. For each corresponding set of four tif images, they are mosaicked, cropped, and projected in ArcGIS. Other input tif files also need to be calibrated for spatial position, projection, and cropping one by one. All processed image data also need to be resampled one by one, and then classified according to the date of each tif image (note to unify the time format of different source tifs), and the grid calculator is used to calculate the monthly average, monthly total, annual average, and annual total. Finally, since the grid calculator does not support simultaneous calculation of multiple formulas, in the calculation process of complex factors, each formula must be input separately, and the intermediate data must be saved manually to continue the next step of calculation. Each step of the traditional calculation method involves different software, with repetitive and cumbersome operations, greatly increasing the operation complexity, time cost, and error rate.

[0127] The operable interface of this fast calculation program is as Figures 3 - 7 shown. The software divides the main functions of ecosystem function service assessment into five modules: MODIS preprocessing, NC preprocessing, water conservation calculation, soil and water conservation calculation, and wind and sand fixation calculation, which can be selected through the tabs at the top of the software ( Figures 3 - 7Red indication part), it should be noted that both the MODIS preprocessing and NC preprocessing modules have TIF preprocessing functions, and the five modules can respectively output intermediate process data for backup. Enter the data path in the corresponding position and click to calculate, and all the output intermediate process data will be in a unified format.

[0128] The evaluation and calculation of the entire ecosystem function service is divided into six parts for efficiency comparison. From the efficiency comparison in Table 3, it can be found that the fast calculation program significantly shortens the processing time when processing various data and executing calculation tasks. Compared with the traditional combined calculation of multiple software, the efficiency of this method has increased by up to 89% to 99%, greatly reducing the time required for calculation (some results are shown as Figure 8 (a)-(e) in the figure). This speed advantage not only benefits from the new algorithms built into this program, but also reflects in the simplification of the operation process. During the traditional combined calculation of multiple software, each operation requires switching between different software. Whether it is data conversion, cropping, or multi-step formula calculation, data needs to be manually input and processed, thus increasing the risk of errors. This method effectively reduces the error probability and ensures the accuracy and consistency of the results by simplifying the operation process and reducing human operation links. In addition, the traditional method requires users to be familiar with the operation methods and interface logics of multiple software, and the learning and exploration costs are high. However, this method integrates complex operations into a unified system, greatly reducing the learning cost of users and improving the overall work efficiency.

[0129] Table 3

[0130]

[0131] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on such an understanding, the above technical solutions, in essence, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disc, etc., and includes several instructions to enable a computer device (which can be a personal computer, server, or network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0132] In this article, specific examples are used to elaborate on the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method and its core idea of the present invention; at the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation methods and application scopes. In summary, the content of this specification should not be construed as a limitation to the present invention."

Claims

1. An automated and rapid calculation method for evaluating ecosystem function services, characterized in that, It includes the following steps: S1. Preset the spatial basic grid and coordinate system covering the research area, obtain meteorological text data, remote sensing data and other raster data, and perform various preprocessing operations on the above three types of data to unify them into the same data specification; S2. Interpolate or resample the preprocessed meteorological text data, remote sensing data and other raster data to unify the spatial values into standard grid data; S3. According to the timestamps in the standard grid data names, perform monthly addition or average calculation, unify the sampling time, and delete all spatial information to convert it into matrix data; S4. For the matrix data, calculate the water conservation capacity, soil and water conservation capacity, and wind prevention and sand fixation capacity according to the models of water conservation, soil and water conservation, and wind prevention and sand fixation respectively, and output the result matrices respectively; S5. Add the spatial information of the preset coordinates to the result matrix and convert it into GeoTif format for output.

2. The automated rapid calculation method for evaluating ecosystem function services according to claim 1, wherein In S2, when interpolating or resampling the preprocessed meteorological text data, it specifically includes: Obtain the meteorological text data in NetCDF format, convert the meteorological text data in NetCDF format to JSON format, extract meteorological parameters, and structurally store them as NCDS = [x, y, time, value], where x and y are the longitude and latitude coordinates of the data points, time is the timestamp, and value is the climate value of the data points; use the datetime module to unify the timestamp to daily granularity to generate a daily climate text data file NCD; spatialize the NCD file into a set of coordinate points with values, and use the inverse distance weighted interpolation method to interpolate the meteorological point data to the spatial basic grid to generate a daily climate grid dataset.

3. The automated rapid calculation method for evaluating ecosystem function services according to claim 1, wherein In S2, when interpolating or resampling the preprocessed remote sensing data, it specifically includes: Obtain the remote sensing data in MODIS format, parse the file name of the remote sensing data, extract the product code, shooting time and satellite orbit row and column numbers; automatically combine the names of the remote sensing data according to the satellite orbit row and column numbers and time interval of the research area to generate a name dictionary; batch read the remote sensing data according to the name dictionary and splice them according to the satellite orbit row and column numbers to generate a daily remote sensing image set; use the nearest neighbor interpolation method to resample the daily remote sensing image set to the spatial basic grid to generate a regular raster dataset.

4. The automated rapid calculation method for evaluating ecosystem function services according to claim 1, wherein In S2, when interpolating or resampling the preprocessed other raster data, it specifically includes: Obtain the other raster data in tif format, and use the nearest neighbor interpolation method to resample the tif format data to the spatial basic grid to generate a regular grid dataset.

5. The automated rapid calculation method for evaluating ecosystem function services according to claim 1, characterized in that In S3, according to the timestamps in the standard grid data names, performing monthly addition or average calculation specifically includes: The standard grid data includes a daily climate grid dataset, a regular raster dataset, and a regular grid dataset. Convert the daily climate grid dataset, the regular raster dataset, and the regular grid dataset into Numpy matrices and remove the spatial information; according to the timestamps in the data names, use the datetime module to calculate the monthly average, monthly sum, annual average or annual sum to generate a matrix dataset with a unified time granularity.

6. An automated rapid calculation method for evaluating ecosystem function services according to claim 1, characterized in that, In S4, for the matrix data, the relevant calculations based on the water conservation model include: The calculation of the surface runoff factor is as follows: R = P × α Wherein, R represents the annual surface runoff, P represents the annual precipitation, and α is the average surface runoff coefficient; the calculation of the water conservation amount is as follows: TQ = (P - R - ET) × A × 10 3 Wherein, TQ is the water conservation amount, P is the annual precipitation, ET is the actual evapotranspiration, R is the surface runoff, and A is the grid area.

7. An automated rapid calculation method for evaluating ecosystem function services according to claim 1, characterized in that In S4, for the matrix data, the relevant calculations based on the soil and water conservation model include: The calculation of the rainfall erosivity is as follows: where R is the average annual rainfall erosivity, P i is the total monthly precipitation, and P is the total annual precipitation; i is the month of the precipitation used, i = 1, 2,..., m; the soil erodibility factor is calculated as follows: K = (-0.01383 + 0.51575K EPIC ) × 0.1317 where K EPIC represents the soil erodibility factor before correction, K represents the soil erodibility factor after correction, m c , m silt , m s and orgC are the proportions of clay, silt, sand and organic carbon respectively; the soil and water conservation calculation is as follows: A c = A p -A r = R × K × L × S × (1 - C) In the formula, A c is the soil and water conservation amount, A p is the potential soil erosion amount, A r is the actual soil erosion amount, R is the rainfall erosion factor, K is the soil erodibility factor; L is the slope length factor, S represents the slope gradient factor, and C is the vegetation cover factor.

8. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements an automated rapid calculation method for evaluating ecosystem function services as described in any one of claims 1 to 7.

Images