Developing an Automated Python Surface Energy Balance System (PySEBS) Software for Calculating Actual Evapotranspiration-Software Development and Application Case in Jilin Province, China

: In this study, Python Surface Energy Balance System (PySEBS) software was developed in the Python 2.7 programming language for continuous calculation of actual evapotranspiration (ET a ) at regional scales. The software is based on the Surface Energy Balance System (SEBS) model, which uses basic meteorological data, MODIS remote sensing data, and Digital Elevation Model (DEM) data as the original input data and ﬁnally outputs daily-scale ET a in the form of raster data with a spatial resolution of 1 km × 1 km. To verify the reliability of the PySEBS model, the ET a of spring maize during the growing season in Jilin Province, China, from 2001 to 2020 was calculated and analyzed in this study and compared with the results of similar studies by others. The ﬁndings showed that the PySEBS model has a reasonable accuracy in estimating ET a within ± 15% and is a robust model that can achieve the continuous calculation of ET a at a regional scale. Therefore, PySEBS software is a useful tool for regional irrigation scheduling and water resources management.


Introduction
Actual evapotranspiration (ET a ) is the combined term for ground evaporation and vegetation transpiration. ET a is an important component of water balance for agricultural fields and is widely used in many areas, such as agronomy, hydrology, ecology, meteorology, and environmental science [1][2][3][4]. The influencing factors include meteorological factors such as solar radiation, air temperature, relative humidity, and wind speed. In addition, it is also related to vegetation type, growth condition, and soil properties [5]. The concept closely related to ET a is the reference crop evapotranspiration (ET o ) [6,7]. The FAO Penman-Monteith approach was developed by defining the ET o as evapotranspiration from a broad surface of highly uniform, vigorously growing, well-watered green grass, assuming a surface resistance of 70 s m −1 , a height of 0.12 m, and an albedo of 0.23 [7].
Since the establishment of Dalton's theorem in 1802, different methods for estimating ET a have emerged. Bowen et al. [8] proposed the Bowen ratio-energy balance method for ET a calculation in 1926, according to the energy balance equation. Thornthwaite and Holzman [9] proposed an aerodynamic method for calculating ET a based on the boundary layer similarity theory. In 1948, Penman proposed the concept of "evaporation force" and established the joint ET a equation [10]. Monteith [11] derived the famous Penman-Monteith equation in 1965, which opened up a new way to study evaporation from an unsaturated subsurface. Allen et al. [7] used the method of assigning a reference crop canopy resistance and a fixed height to give a careful description of the use of the P-M formula under different environmental conditions. Subsequently, many studies of stratified models of ET a were carried out based on the energy balance principle. The singlelayer model mainly uses the energy balance residual method, which does not distinguish between soil evaporation and vegetation transpiration in the estimation process. It uses a holistic approach to estimate ET a , but this method usually results in overestimation under sparse vegetation. Therefore, Shuttleworth and Wallace [12] proposed a two-layer model based on the theoretical principle of the P-M formula, in which soil evaporation and vegetation transpiration are calculated separately. Based on Shuttleworth's two-layer model, Dolman [13] proposed a three-layer model, and Choudhury et al. [14] put forward a four-layer model. Norman et al. [15] developed a two-source model for row crops. All the above methods for calculating ET a are mainly based on point-source data, which is relatively accurate on small scales. However, in practice, it is often necessary to know the ET a on a large scale, and the traditional methods can hardly meet the requirements. The gradual rise of remote sensing technology in the late 1960s, especially the application and development of thermal infrared remote sensing, can achieve the conversion from point to surface, which makes the large-scale study of ET a possible.
A series of remote sensing models for accurate estimation of ET a have emerged, such as the triangle method model, Surface Energy Balance Algorithm for Land (SEBAL), Surface Energy Balance Systems (SEBS), Mapping Evapotranspiration with Internalized Calibration (METRIC), Two source Trapezoid Model for Evapotranspiration (TTME), Google Earth Engine SEBAL (geeSEBAL), etc. [16][17][18][19][20][21]. Du et al. [22] used the SEBAL model, combined with remote sensing data and meteorological data, to estimate the water productivity and specific water consumption of rice in the Songnen Plain. Bhattarai and Liu [23] developed a Matlab-based toolbox named LandMOD ET mapper to run SEBAL and METRIC on Landsat and MODIS images for mapping field to regional scale ET a . The SEBALIGEE model developed by Mhawej and Faour [24] is an automated, user-friendly, flexible, and openaccess 30 m ET retrieval system. Its main features are the estimations of ET and K c in diverse climatic regions with readily available Landsat-8 30-m satellite images. Ramirez-Cuesta et al. [25] developed a new ArcGIS toolbox for implementing the METRIC energy balance model, called Metric-GIS, and tested it in a semi-arid environment. Ellsasser et al. [26] developed the QWaterModel tool based on the DATTUTDUT energy balance model using land surface temperature maps as input, making ET predictions available to a wider audience. Wu et al. [27] introduced an ETWatch Cloud service, which is a model collection migrated from a local PC to the cloud platform for ET generation using readily available RS and standard ground-based meteorological observations. Laipelt et al. [21,28] developed the geeSEBAL software in the Google Earth Engine environment based on the SEBAL algorithm and achieved better results in the calculation of evapotranspiration in the Brazilian region. Huang et al. [29] and Ren et al. [30] used the SEBS model, combined with MODIS data, to estimate the water productivity and water consumption of maize. The results showed that the improved SEBS model could be well applied to the study of ET a in dryland farming.
Although there are some tools for calculating ET a at a regional scale, such as ILWIS (Integrated Land and Water Information System) software, the operation process is complicated for users [31]. The purpose of this study is to introduce the Python Surface Energy Balance System (PySEBS) software for calculating ET a at a regional scale to simplify the complex process of ET a calculation and to achieve multi-year, continuous, automatic, and efficient estimation of ET a . Moreover, a case study is performed, and the results of calculated ET a using this software are described, analyzed, and compared to similar studies.

Installing the PySEBS and Software Requirements
The PySEBS is software that can be used to calculate actual evapotranspiration (ET a ) on a spatial scale based on the Surface Energy Balance System model (SEBS). The software is written in Python 2.7.14 due to Python (2.7) being the latest version supported by ArcGIS Remote Sens. 2022, 14, 5629 3 of 22 10.1-10.6. Since ET a is calculated on a spatial scale, the software operation utilizes some of the built-in tools in ArcGIS, and the proper operation of these tools depends on the calling of the Arcpy package in ArcGIS. Therefore, the user needs to ensure that the ArcGIS software is installed successfully.
The PySEBS software is compressed into a PySEBS V1.0.0.zip file. The schematic diagram of the PySEBS file structure is shown in Figure 1. The schematic box labeled. /PySEBS/represents a folder named PySEBS in the working directory, where the PySEBS 1.0.0.zip file will be stored after extraction. The PySEBS folder contains two parts: the Py-SEBS V1.0.0.exe program; and the site-packages folder containing the code libraries (arcpy package, numpy package, pandas package, Pyshp package, etc.) and the configuration Desktop.pth file.
The PySEBS is software that can be used to calculate actual evapotranspiratio on a spatial scale based on the Surface Energy Balance System model (SEBS). The so is written in Python 2.7.14 due to Python (2.7) being the latest version suppor ArcGIS 10.1-10.6. Since ETa is calculated on a spatial scale, the software operation some of the built-in tools in ArcGIS, and the proper operation of these tools depe the calling of the Arcpy package in ArcGIS. Therefore, the user needs to ensure t ArcGIS software is installed successfully.
The PySEBS software is compressed into a PySEBS V1.0.0.zip file. The schem agram of the PySEBS file structure is shown in Figure 1. The schematic box l /PySEBS/represents a folder named PySEBS in the working directory, where the P 1.0.0.zip file will be stored after extraction. The PySEBS folder contains two pa PySEBS V1.0.0.exe program; and the site-packages folder containing the code li (arcpy package, numpy package, pandas package, Pyshp package, etc.) and the c ration Desktop.pth file.

Graphical User Interface Introduction
The PySEBS software is a concise and straightforward scientific computing contains one main interface and three sub-interfaces ( Figure 2). The main interface c of three button controls: the remote sensing data processing module button, the m logical data processing module button, and the PySEBS calculation and processin ule button. Each button control accesses the corresponding module (sub-interface) ule 1-Remote Sensing Data Processing interface, Module 2-Meteorological Da cessing interface, and Module 3-PySEBS Calculation Procedure interface.
The main functions of Module 1 and Module 2 are to provide properly process formatted input data for use in Module 3. Module 3 is used to calculate ET daily/monthly/annual scale, or for any period of time within a year and across based on the input data (raster data) prepared by Modules 1 and 2. There is no pro priority between Modules 1 and 2. Moreover, there are common properties configu all three modules, including the time setting of the start and end year, file and folde way settings, selection of processing method, and information output.

Graphical User Interface Introduction
The PySEBS software is a concise and straightforward scientific computing tool. It contains one main interface and three sub-interfaces ( Figure 2). The main interface consists of three button controls: the remote sensing data processing module button, the meteorological data processing module button, and the PySEBS calculation and processing module button. Each button control accesses the corresponding module (sub-interface): Module 1-Remote Sensing Data Processing interface, Module 2-Meteorological Data Processing interface, and Module 3-PySEBS Calculation Procedure interface.
The main functions of Module 1 and Module 2 are to provide properly processed and formatted input data for use in Module 3. Module 3 is used to calculate ET a on a daily/monthly/annual scale, or for any period of time within a year and across years, based on the input data (raster data) prepared by Modules 1 and 2. There is no processing priority between Modules 1 and 2. Moreover, there are common properties configured for all three modules, including the time setting of the start and end year, file and folder pathway settings, selection of processing method, and information output.
Each module is described below, and detailed descriptions are available in the PySEBS V1.0.0 Actual Evapotranspiration (ET a ) Calculator User's Guide in the Supplementary Materials: Module 1 is used to process remotely sensed data and generate the daily data files in .asc format for input into Module 3. Module 1 can only process data on an annual basis, and therefore the user must arrange the downloaded remote sensing files in specific folders by data type and by year. The downloaded remote-sensed data needs to be pre-processed before it can be used in the PySEBS program. For demonstration, remotely sensed data were obtained from the National Aeronautics and Space Administration (NASA, http://modis-land.gsfc.nasa.gov/ (accessed on 9 March 2019)). In this example, images h26v04 and h27v04 were downloaded, and the MRT (MODIS Reprojection Tool) projection tool was applied to preprocess the data for merging, resampling, format conversion, etc., and Cygwin software was used for batch processing (see Table 1 for specific parameters). The purpose of pre-processing is to convert the remote sensing data from the hierarchical data format (.hdf) to raster files (.tif). The location of the .tif file will be in the same folder as the .hdf file, and the name of the preprocessed .tif file must contain the date information string (seven digits) of the corresponding data; an example file name is MOD09A1_2019001.sur_refl_b01.tif. Module 1 is used to process remotely sensed data and generate the daily data files in .asc format for input into Module 3. Module 1 can only process data on an annual basis, and therefore the user must arrange the downloaded remote sensing files in specific folders by data type and by year. The downloaded remote-sensed data needs to be pre-pro-  The function of Module 2 is to process the meteorological data and generate the daily file in the .asc format that is required by Module 3. The meteorological data used in this study is China's surface climate data daily value data set (V3.0), which is derived from the National Meteorological Information Center (http://data.cma.cn/site/index.html (accessed on 30 May 2018)); the platform contains the daily meteorological data from 824 reference weather stations. The meteorological data in the dataset includes atmospheric pressure, air temperature, precipitation, evaporation, sunshine hours, surface temperature, wind speed and direction, and relative humidity at each station. The downloaded original meteorological data is in .txt format and includes all meteorological stations in each folder that need to be preprocessed. Once the downloaded V3.0 weather datasets are placed in the same folder, Module 2 can be used to process the raw data directly. If meteorological datasets from other sources are used, they need to be organized in .csv format and stored in the 1RawData folder. An example of the corresponding path may be: .\PySEBS\ETa\2MeteorologicalData\1RawData\2019.csv.
Module 3 uses the daily input data (in .asc format) generated by Module 1 and Module 2 to calculate ET a on a daily/monthly/annual basis, or for a specified period of time. The PySEBS model requires 15 necessary parameters. Module 1 provides eight essential parameters, including seven parameters from MOD09A1 (Band1, Band2, Band3, Band4, Band5, Band7, and Szen), and one parameter from MOD11A2 (Lst). Module 2 provides seven essential parameters, including the Digital Elevation Model (DEM) parameter and six meteorologically related parameters (Tem, Ws, Pa, Roh, Rs, and Rnl). Moreover, there are two non-essential parameters, Pre and ET o , which are not involved in the PySEBS calculation process.

PySEBS Software Theoretical Foundation
The SEBS model, as proposed by Su, is a single-layer model based on the surface energy balance principle [18], which consists of four components: (1) The determination of the land surface physical parameters based on the remotely sensed spectral reflectance and radiance; (2) The estimation of the roughness length for heat transfer; (3) The calculation of sensible heat flux; and (4) The calculation of latent heat flux. Compared with other single-layer models based on the energy balance principle, the PySEBS model has the advantage that the data in each grid is calculated independently, so that even if some grids have missing data on a certain day due to the influence of rain and clouds, the calculation in other grids will not be affected; this can maximize the utilization of the remote sensing data.
The net earth's surface radiation is the basis of various energy exchanges and can be divided into three parts: (1) the soil heat flux that warms the soil; (2) the sensible heat flux that warms the atmosphere; and (3) the latent heat flux of evapotranspiration (the latent heat of soil evaporation and crop transpiration). The main equations for surface energy balance are as follows, and the detailed derivation of the equations can be found in the relevant literature reviewed [7,18,[31][32][33][34][35][36][37][38][39][40][41][42] and in the Supplementary Materials: where R n is the net radiation (W m −2 ), G 0 is the soil heat flux (W m −2 ), H is the sensible heat flux (W m −2 ), and λE is the latent heat flux (W m −2 ), in which λ is the latent heat of evaporation (J kg −3 ) and E is the actual evapotranspiration (mm d −1 ). After calculating each parameter of the surface energy balance separately, the actual daily ET a can be obtained according to the constant evaporation ratio within one day. Calculations of net radiation R n : where α is surface albedo, a ratio of radiation reflected by the surface to global radiation (a sum of direct solar radiation and scattered radiation); R swd and R lwd (W m −2 ) represent downward solar radiation and downward long-wave radiation, respectively; ε is surface emissivity and can be classified into ε land (surface emissivity of land) and ε water (surface emissivity of water) according to different underlying surfaces; σ is the Stephan Boltzmann constant equal to 5.67 × 10 −8 (W m −2 K −2 ); and T 0 is the surface temperature (K).
Calculations of soil heat flux G 0 : where Γ c is the ratio of soil heat flux to net radiation in the condition of total vegetation cover, (Γ c = 0.05), and Γ s is a ratio of soil heat flux to net radiation in the condition of bare soil (Γ s = 0.315); ƒ c is vegetation coverage inversed through remote sensing data.
Calculations of sensible heat flux H: where ρ for air density (kg m −3 ); C p for constant-pressure specific heat (J kg −1 K −1 ); θ 0 for latent surface temperature (K); θ a for latent air temperature (K) at an elevation of z; r a for external resistance; k for Karman constant (k = 0.4); u * for friction velocity (m s −1 ); z for reference measured elevation (m), generally between 2 and 10 m; d 0 for zero plane displacement height (m); z 0h for scalar roughness height for heat transfer (m); Ψ h for stability correction coefficients of heat exchange MOS; and L for Obukhov length (m). Calculations of daily evapotranspiration E daily : The previously calculated net surface radiation, soil heat flux, and sensible heat flux are based on the instantaneous values at the moment of satellite transit. In contrast, in the PySEBS model, the daily evapotranspiration is calculated from the evaporation ratio with the daily average net radiation and daily average soil heat flux. The evapotranspiration ratio remains stable throughout the day [43][44][45][46], and assuming that the evapotranspiration ratio is constant throughout the day, the daily evapotranspiration is calculated by the following equation: In the above equation, E daily stands for the actual daily evapotranspiration (mm day −1 ); 8.64·10 7 stands for a unit conversion coefficient; 24 Λ 0 stands for the daily mean evaporation ratio that can be substituted by the evaporation ratio Λ figured out by the PySEBS model; R n stands for the daily mean net radiation flux (W m −2 ); and G 0 stands for the daily mean soil heat flux (W m −2 ), considering that the amount of heat absorbed and released by soil is approximately equal during the day, the daily mean soil heat flux is set to 0. Additionally, λ refers to the latent heat of vaporization of water and λ = 2.45·10 6 (J kg −1 ); ρ w refers to the density of water and ρ w = 1000 (kg m −3 ); α is the surface albedo; K ↓ 24 refers to incident global radiation (W m −2 ), which is obtained through the unit conversion of the incident solar radiation (MJ m −2 day −1 ) expressed in the Penman-Monteith equation; ε is surface emissivity; and L 24 refers to daily net outgoing long-wave radiation (W m −2 ), which is obtained through the unit conversion of the net long-wave radiation (MJ m −2 day −1 ) expressed in the Penman-Monteith equation.

Initial Input Data
Input data required by Modules 1, 2, and 3 can be divided into three parts: the first part is the basic data, which includes the information file (V3.0_station_jls.csv) of weather stations for the V3.0 meteorological dataset within the study area, the study area DEM file (dem_jls.asc), and the geographic coordinate system file (WGS 1984.prj); the second part is the remote sensing data, which includes MOD09A1 and MOD11A2; the third part is the meteorological data, which includes the raw data of the V3.0 meteorological dataset. Meteorological data from other sources can also be used; the data needs to be prepared by the users following the format of 2019.csv and located in the corresponding folder.
The input data or files described above are all saved in the designated folders. Each module is connected to the file or folder so that the data can serve as the input for the module. Specifically, the basic data is connected in the form of files; remotely sensed data and meteorological data are connected in the form of subdirectory.
The structure diagram of the path to the initial input data is shown in Figure 3.

Output Results
The ./PySEBS/ETa/ file path is a user-defined folder to save the output results, which will be created automatically if such a folder does not exist in the project directory. The folder 1RemoteSensingData/ is used to save the intermediate files generated during the remote sensing data processing in Module 1. The folder 2MeteorologicalData/ is used to save the intermediate files generated during the meteorological data processing in Module 2. The folder 3PySEBSInputData/ is used to save the final processing results of remote sensing and meteorological data. In addition, the folder 3PySEBSInputData/ also serves as the path for input files for Module 3. The 4PySEBSOutputData/ folder will store the daily ETa results, and the 5sum_ETa/ folder will store the monthly and annual ETa results calculated by Module 3 (Figure 4). The 6sum_ETo/ folder and the 7sum_Pre/ folder are used to store cumulative reference evapotranspiration and precipitation, respectively.
As mentioned earlier, if meteorological data derived from other sources need to be used, the data must be prepared in the required format ( Figure 5) and stored in the file path 2Meteorologicaldata/1RawData/ to replace the V3.0 meteorological data. The files in the file path 2MeteorologicalData/2SourceData/ (e.g., ETo_2019.csv) used the data in the

Output Results
The ./PySEBS/ETa/ file path is a user-defined folder to save the output results, which will be created automatically if such a folder does not exist in the project directory. The folder 1RemoteSensingData/ is used to save the intermediate files generated during the remote sensing data processing in Module 1. The folder 2MeteorologicalData/ is used to save the intermediate files generated during the meteorological data processing in Module 2. The folder 3PySEBSInputData/ is used to save the final processing results of remote sensing and meteorological data. In addition, the folder 3PySEBSInputData/ also serves as the path for input files for Module 3. The 4PySEBSOutputData/ folder will store the daily ET a results, and the 5sum_ETa/ folder will store the monthly and annual ET a results calculated by Module 3 (Figure 4). The 6sum_ETo/ folder and the 7sum_Pre/ folder are used to store cumulative reference evapotranspiration and precipitation, respectively.   As mentioned earlier, if meteorological data derived from other sources need to be used, the data must be prepared in the required format ( Figure 5) and stored in the file path 2Meteorologicaldata/1RawData/ to replace the V3.0 meteorological data. The files in the file path 2MeteorologicalData/2SourceData/ (e.g., ETo_2019.csv) used the data in the 2Meteorologicaldata/1RawData/ folder as input and calculated the ET o and its process variables; the purpose of these files is to provide the six required meteorological parameters for Module 3. The calculation of ET o is also a strength of this software, if the users need to use ET o and other parameters in the calculation process, the information can be found in the folder 2MeteorologicalData/2SourceData/.   The 3PySEBSInputData/ folder stores the final processing results calculated from Module 1 and Module 2 and the input files from Module 3. The files are stored in the form of one folder per day and contain 15 .asc files that are required for the PySEBS model calculation and one optional precipitation file (pre.asc). The required parameters include one DEM file (dem.asc), seven MOD09A1 files (band1.asc, band2.asc, band3.asc, band4.asc, band5.asc, band7.asc, and szen.asc), one MOD11A2 file (lst.asc), and six meteorological parameter files (tem.asc, ws.asc, pa.asc, roh.asc, rs.asc, and rnl.asc) ( Figure 6).

PySEBS Case Study-Study Area
After the development of PySEBS software, the ET a was calculated for spring maize in Jilin Province from 2001 to 2020, and the ET a was analyzed at both temporal and spatial scales during the growing season of spring maize in this study. Jilin Province is located in the central part of northeastern China (Figure 7). The terrain map shows obvious characteristics of high elevations in the southeast and low elevations in the northwest. The latitude varies from 40 • 50 to 46 • 19 N, and the longitude ranges from 121 • 38 to 131 • 19 E. Jilin Province is 769.6 km long from east to west and 606.6 km wide from north to south and has an area of 187,400 km 2 . As for climate conditions, it belongs to a temperate continental monsoon climate with four distinct seasons. Rain and heat are in the same season, such as a rainy summer with high temperatures and a cold and dry winter. The average winter temperature is below −11 • C, and the average summer temperature is above 23 • C. The annual average sunshine hours are 2259-3016 h. The annual average precipitation is 400-600 mm, and the annual frost-free period is 120-160 days. The climatic conditions of Jilin Province make it an important commodity grain production region in China. Due to the geographical advantage for agriculture, it has become one of the three famous Corn Belts in the world.

PySEBS Case Study-Study Area
After the development of PySEBS software, the ETa was calculated for spring maize in Jilin Province from 2001 to 2020, and the ETa was analyzed at both temporal and spatial scales during the growing season of spring maize in this study. Jilin Province is located in the central part of northeastern China (Figure 7). The terrain map shows obvious characteristics of high elevations in the southeast and low elevations in the northwest. The latitude varies from 40°50′ to 46°19′ N, and the longitude ranges from 121°38′ to 131°19′ E.

Calculation of ETa during Spring Maize Growing Season
The ETa of spring maize was calculated by PySEBS software. The data used were remote sensing data, meteorological data, and DEM data. The remote sensing data were obtained from NASA. The meteorological data are the daily data set of China's surface climate data (V3.0), which were accessed from the Chinese National Meteorological Information Center. The digital elevation data ASTER GDEM V2 were downloaded from the geospatial data cloud (http://www.gscloud.cn/ (accessed on 16 May 2018)) with a spatial resolution of 30 m × 30 m.
In the process of calculating ETa through the PySEBS model in this study, the weather data were synthesized for 8 days due to the remote sensing data of MOD09A1 and MOD11A2 used are 8-day synthesis products. Then, the weather data and remote sensing data with a unified time scale were used as the input parameters of the PySEBS model. Finally, the 8-day synthetic ETa was obtained. Then, according to the start and end time of the growth period of spring maize (1 May to 30 September), the ETa during the growth period was accumulated by adding up.

Extraction of Planting Area
The PySEBS software is a software for the calculation of the ETa on a spatial scale, and the results of the calculation are finally presented in the form of a raster dataset. Since the subject of research in this study is spring maize in Jilin Province, the planting area of spring maize in the study area needs to be determined first.
According to the distinct spectral characteristics of spring maize in different growth stages, the Normalized Difference Vegetation Index (NDVI) can effectively reflect the differences between various crops. Based on the curve of the MODIS-NDVI time series and the land use data of Jilin Province over 16 years, a decision tree classification model was constructed to achieve the extraction of the spring maize planting area. To obtain time series curves of NDVI of spring maize in Jilin Province for all years, the NDVI values of each sampling site were extracted as the main basis for extracting the spring maize

Calculation of ET a during Spring Maize Growing Season
The ET a of spring maize was calculated by PySEBS software. The data used were remote sensing data, meteorological data, and DEM data. The remote sensing data were obtained from NASA. The meteorological data are the daily data set of China's surface climate data (V3.0), which were accessed from the Chinese National Meteorological Information Center. The digital elevation data ASTER GDEM V2 were downloaded from the geospatial data cloud (http://www.gscloud.cn/ (accessed on 16 May 2018)) with a spatial resolution of 30 m × 30 m.
In the process of calculating ET a through the PySEBS model in this study, the weather data were synthesized for 8 days due to the remote sensing data of MOD09A1 and MOD11A2 used are 8-day synthesis products. Then, the weather data and remote sensing data with a unified time scale were used as the input parameters of the PySEBS model. Finally, the 8-day synthetic ET a was obtained. Then, according to the start and end time of the growth period of spring maize (1 May to 30 September), the ET a during the growth period was accumulated by adding up.

Extraction of Planting Area
The PySEBS software is a software for the calculation of the ET a on a spatial scale, and the results of the calculation are finally presented in the form of a raster dataset. Since the subject of research in this study is spring maize in Jilin Province, the planting area of spring maize in the study area needs to be determined first.
According to the distinct spectral characteristics of spring maize in different growth stages, the Normalized Difference Vegetation Index (NDVI) can effectively reflect the differences between various crops. Based on the curve of the MODIS-NDVI time series and the land use data of Jilin Province over 16 years, a decision tree classification model was constructed to achieve the extraction of the spring maize planting area. To obtain time series curves of NDVI of spring maize in Jilin Province for all years, the NDVI values of each sampling site were extracted as the main basis for extracting the spring maize planting area using the coordinate information of 38 sample sites with a long-term spring maize planting in this study (Figure 7).
After the sowing and emergence stages, the NDVI showed an increasing trend, which rose rapidly during the three-leaf, seven-leaf, and jointing stages of the spring maize vegetative growing period. When the spring maize entered the tasseling stage at the peak of vegetative growth, the NDVI value reached its peak. Then, spring maize entered the reproductive growth stages of the silking period, milk maturity period, and maturity period, during which the NDVI values gradually decreased. To eliminate the influence of noise on the NDVI values, the method of Harmonic Analysis of Time Series (HANTS) was used to smooth the time series of MODIS-NDVI. The comparison of before and after processing curves is shown in Figure 8.
Remote Sens. 2022, 14, 5629 12 planting area using the coordinate information of 38 sample sites with a long-term sp maize planting in this study (Figure 7). After the sowing and emergence stages, the NDVI showed an increasing trend, w rose rapidly during the three-leaf, seven-leaf, and jointing stages of the spring maize etative growing period. When the spring maize entered the tasseling stage at the pea vegetative growth, the NDVI value reached its peak. Then, spring maize entered th productive growth stages of the silking period, milk maturity period, and maturity riod, during which the NDVI values gradually decreased. To eliminate the influen noise on the NDVI values, the method of Harmonic Analysis of Time Series (HANTS) used to smooth the time series of MODIS-NDVI. The comparison of before and after cessing curves is shown in Figure 8. The construction of the decision tree for the extraction of the spring maize plan area in Jilin Province was based on the minimum and maximum time series curve NDVI for the period of 2001 to 2016. A total of 13 images of NDVI values reconstru by HANTS using a day-of-year (DOY) calendar from DOY 113 to 305 over 16 years w selected for this study. The pixels between the minimum and maximum NDVI valu the corresponding DOY were the candidate areas for the spring maize planting area eliminate the influence of other land uses, the agricultural land in the land cover provided by the remote sensing data of MCD12Q1 was extracted as a mask. Finally spring maize planting area was obtained. The extraction process is shown in Figure 9  The construction of the decision tree for the extraction of the spring maize planting area in Jilin Province was based on the minimum and maximum time series curves of NDVI for the period of 2001 to 2016. A total of 13 images of NDVI values reconstructed by HANTS using a day-of-year (DOY) calendar from DOY 113 to 305 over 16 years were selected for this study. The pixels between the minimum and maximum NDVI values of the corresponding DOY were the candidate areas for the spring maize planting area. To eliminate the influence of other land uses, the agricultural land in the land cover map provided by the remote sensing data of MCD12Q1 was extracted as a mask. Finally, the spring maize planting area was obtained. The extraction process is shown in Figure 9.
In this study, the identification of the spring maize planting area was based on the NDVI time series curve for the period of    The map of the multi-year average spring maize planting area in Jilin Province was obtained by accumulating the raster files of planting areas in previous years, and a raster was considered to belong to the multi-year average spring maize planting area when the number of years in which a raster appeared exceeded 7 years. This planting area map was used to clip the ET a raster data of Jilin Province calculated by PySEBS software to obtain the ET a of spring maize. The map of the multi-year average spring maize planting area in Jilin Province was obtained by accumulating the raster files of planting areas in previous years, and a raster was considered to belong to the multi-year average spring maize planting area when the number of years in which a raster appeared exceeded 7 years. This planting area map was used to clip the ETa raster data of Jilin Province calculated by PySEBS software to obtain the ETa of spring maize.

Temporal and Spatial Distribution Patterns of ETa
The temporal distribution pattern of ETa during the spring maize growing period (May to September) from 2001 to 2020 in Jilin Province is shown in Figure 11. The average cumulative ETa during the growing period of spring maize in Jilin Province ranged from 475 to 544 mm, with a multi-year average value of 517 mm and an average standard deviation value of 43 mm ( Figure 11A). The lowest and highest ETa happened in 2008 and 2020, respectively. The average annual cumulative ETa differed from 632 to 694 mm, with a mean value of 666 mm, and the standard deviation varied from 50 to 86 mm, with a mean value of 65 mm ( Figure 11B). The monthly distribution pattern of ETa showed a trend of increasing and then decreasing, with the monthly average ETa of 2, 9, 34, 58, 83,

Temporal and Spatial Distribution Patterns of ET a
The temporal distribution pattern of ET a during the spring maize growing period (May to September) from 2001 to 2020 in Jilin Province is shown in Figure 11. The average cumulative ET a during the growing period of spring maize in Jilin Province ranged from 475 to 544 mm, with a multi-year average value of 517 mm and an average standard deviation value of 43 mm ( Figure 11A). The lowest and highest ET a happened in 2008 and 2020, respectively. The average annual cumulative ET a differed from 632 to 694 mm, with a mean value of 666 mm, and the standard deviation varied from 50 to 86 mm, with a mean value of 65 mm ( Figure 11B). The monthly distribution pattern of ET a showed a trend of increasing and then decreasing, with the monthly average ET a of 2,9,34,58,83,110,132,114,78,35,9, and 2 mm from January to December, respectively, with the highest value of ET a in July. The cumulative ET a was 517 mm from May to September during the spring maize growing period ( Figure 11C). value of ETa in July. The cumulative ETa was 517 mm from May to September during the spring maize growing period ( Figure 11C). The spatial distribution of the multi-year average ETa and inter-annual coefficient of variation (CV) of spring maize in Jilin Province from 2001 to 2020 is shown in Figure 12. The multi-year average ETa was between 391 and 693 mm in Jilin Province, and the multiyear average ETa in Changchun City, Siping City, Jilin City, Liaoyuan City, and other regions was between 450 and 550 mm. Farmland near rivers, such as the Dongliao River in Siping City, the Songhua River across Songyuan City and Changchun City, and the Nen River across the border of Baicheng City, had the ETa above 600 mm. Areas with annual average ETa less than 450 mm were mainly distributed in the western region of Jilin Province, such as Baicheng City, Songyuan City, and the eastern part of Siping City. This might be because these areas had mostly sandy soil and insufficient rainfall. The CV of the multiyear average ETa of spring maize was mostly between 2% and 18% in Jilin Province. In a few areas, such as the western part of Baicheng City and Songyuan City, the inter-annual CV exceeded 10%, which might be caused by the large variation in inter-annual rainfall in these areas. Figure 13 shows the distribution of ETa during the growing season of spring The spatial distribution of the multi-year average ET a and inter-annual coefficient of variation (CV) of spring maize in Jilin Province from 2001 to 2020 is shown in Figure 12. The multi-year average ET a was between 391 and 693 mm in Jilin Province, and the multi-year average ET a in Changchun City, Siping City, Jilin City, Liaoyuan City, and other regions was between 450 and 550 mm. Farmland near rivers, such as the Dongliao River in Siping City, the Songhua River across Songyuan City and Changchun City, and the Nen River across the border of Baicheng City, had the ET a above 600 mm. Areas with annual average ET a less than 450 mm were mainly distributed in the western region of Jilin Province, such as Baicheng City, Songyuan City, and the eastern part of Siping City. This might be because these areas had mostly sandy soil and insufficient rainfall. The CV of the multi-year average ET a of spring maize was mostly between 2% and 18% in Jilin Province. In a few areas, such as the western part of Baicheng City and Songyuan City, the inter-annual CV exceeded 10%, which might be caused by the large variation in inter-annual rainfall in these areas.  The ETa values for each city in Jilin Province are shown in Figure 14. ETa for the maize growing season in each city basically remained between 450-550 mm from 2020, and the largest multi-year average ETa was in Tonghua, Jilin, and Liaoyuan w exceeded 530 mm. The smallest ETa was in Yanbian, Songyuan, and Changchun w ETa lower than 510 mm ( Figure 14A). The CV of ETa for the spring maize growing in each city ranged from 5% to 8%, with the smallest variability being 5% in Tongh the largest variability being 7.8% in Baicheng ( Figure 14B). The distribution of month in each city is shown in Figure 14C. ETa was mainly concentrated in May-S ber, i.e., during the growing period of spring maize, and reached its maximum with an average of approximately 150 mm.  The ET a values for each city in Jilin Province are shown in Figure 14. ET a for the spring maize growing season in each city basically remained between 450-550 mm from 2001 to 2020, and the largest multi-year average ET a was in Tonghua, Jilin, and Liaoyuan with ET a exceeded 530 mm. The smallest ET a was in Yanbian, Songyuan, and Changchun with an ET a lower than 510 mm ( Figure 14A). The CV of ET a for the spring maize growing season in each city ranged from 5% to 8%, with the smallest variability being 5% in Tonghua and the largest variability being 7.8% in Baicheng ( Figure 14B). The distribution of ET a by month in each city is shown in Figure 14C. ET a was mainly concentrated in May-September, i.e., during the growing period of spring maize, and reached its maximum in July with an average of approximately 150 mm.

Comparison of the Calculated ETa Results with Similar Studies
To verify the rationality of the ETa calculation results, this study has gathered and sorted out related studies for comparison ( Table 2). Guo et al. [47] estimated the ETa of spring maize in Yushu City, Jilin Province, in 2013, using a large weighing lysimeter on a farm scale, and found an average value of 362 mm for ETa. Qiu et al. [48] and Zhang et al

Comparison of the Calculated ET a Results with Similar Studies
To verify the rationality of the ET a calculation results, this study has gathered and sorted out related studies for comparison (Table 2). Guo et al. [47] estimated the ET a of spring maize in Yushu City, Jilin Province, in 2013, using a large weighing lysimeter on a farm scale, and found an average value of 362 mm for ET a . Qiu et al. [48] and Zhang et al. [49] used the Penman-Monteith and crop coefficient method to estimate the 3-year and 54-year ET a of spring maize in Jilin Province, and the average values were 523 mm and 538 mm, respectively, which were similar to the results of this study. Liu [50] used the SIMETAW (Simulation Evapotranspiration of Applied Water) model to simulate the evapotranspiration of the spring maize growing season under different hydrological year types in western Liaoning from 2007 to 2009, and the results were 464 mm, 485 mm, and 486 mm in abundant, flat, and dry years, respectively. In this study, the ET a values of spring maize in Jilin Province from 2001 to 2020 were between 391 mm and 693 mm, with a multi-year average of 517 mm. This was similar to the results of other related studies in Jilin Province. Jiang et al. [51] pointed out that the maximum error value of ET a estimated by using remote sensing data is reasonable within ±20%. However, Seguin et al. [52] suggested that the maximum error of ET a estimated by remote sensing changes with the actual applications, and the absolute value is usually reasonable around 15-30%. Since there is no relevant flux tower data in the study area that can be used directly for the validation of this model for ET a estimations, a relevant literature review was conducted to validate the estimation effectiveness of this PySEBS model on ET a in Jilin Province. The results of relevant literature studies on spring maize ET a in Jilin Province at the regional scale were 514 mm on average, while the average result of the PySEBS model was 517 mm, with a relative error of less than 1%, indicating a good applicability of the modified model in regional studies (Table 2).
PySEBS software has a high application value in water use evaluation and irrigation management. Based on this model, Ren et al. [30] studied the spatio-temporal patterns of water consumption and irrigation demand during winter wheat-summer maize agricultural production in the Huang-Huai-Hai Plain (3HP), and the scenario analysis of different cropping methods was evaluated. The annual average ET a of winter wheat-summer maize ranged from 700-900 mm with no significant change, while irrigation water demand showed a large interannual variation. In addition, four alternatives that may require less irrigation were tested, i.e., wheat-maize-wheat-fallow (WMWF), wheat-maize-fallow-maize (WMFM), wheat-maize-fallow-spring maize (WMFSM), and fallow-spring maize-fallowspring maize (FSMFSM), with the results showing that FSMFSM has the best water savings, reducing irrigation by 64%.

Advantages and Disadvantages of PySEBS Software
Before the development of PySEBS software, studies were conducted by calculating regional ET a using the SEBS module embedded in ILWIS software (https://52north.org/ software/software-projects/ilwis (accessed on 1 April 2019)), but the SEBS module has many shortcomings, such as a tedious calculation process, inefficient calculation, only one period of results can be calculated at one time, a lack of batch processing capability, and various coefficient transformations of data are required during processing to make the data into integer form for further operations, which is more challenging for beginners who need to calculate ET a products urgently, and unfavorable for ET a estimation at a large spatial and temporal scales.
To overcome the aforementioned disadvantages, this study developed a more userfriendly PySEBS, which has the following advantages compared with the conventional calculation method. The strengths are to (1) simplify the preprocessing process of remote sensing data and meteorological data, especially meteorological data can be directly interpolated from point-scale data to raster files; (2) develop a concise user interface that makes the input and output of data and modification of parameters more straightforward and convenient; and (3) reduce the time consumption greatly. Once the required parameters are set, the whole process will no longer require the users to manually input data year by year, and the ET a results can be calculated automatically and continuously for many years on a regional scale. In addition, the newly developed PySEBS model can be used to calculate crop water productivity [30], which is very useful for land degradation assessment [53]. Although the PySEBS model performed relatively well in this study, there were some limitations: since the cover type primarily targeted in this study was spring maize, the calculations of ET a in estimating other crops might be inadequate; another point is that the spatial and temporal resolutions of the remote sensing data used need to be further improved, and the spatial resolution of 500-1000 m is not sufficient to distinguish the evapotranspiration characteristics in areas with complex surface cover types, and the existence of mixed image elements may affect the estimation results to a certain extent. Finally, users can also use remote sensing data from other sources, such as Landsat and Sentinel satellites, rather than being limited to MODIS remote sensing data. However, some of the relevant equations in the Supplementary Materials, for example, Equation (3) for calculating surface albedo needs to be adjusted accordingly. Future research should pay more attention to the estimation of evapotranspiration on different surface cover types and the selection of remote sensing data with a higher resolution.

Conclusions
A scientific calculation software named PySEBS was developed to calculate ET a on a regional scale. Its most important feature is that it can preprocess remote sensing data and meteorological data quickly, efficiently, and continuously, and calculate the regional ET a on a daily basis. Detailed descriptions of the installation environment, input data requirements, selection method of meteorological data interpolation, formula sources and derivations, and operation procedures of the PySEBS software are presented in this study, which walk the users through the whole set-up and help them to understand the software principles. In addition, spring maize in Jilin Province was selected as the research object, and ET a was calculated during the growing season of spring maize from 2001 to 2020. Finally, the spatio-temporal distribution patterns of ET a during the growing season of spring maize in Jilin Province were mapped and analyzed in this study.
Results indicated that the multi-year average extraction accuracy of the spring maize planting area was 90.6% and that the overall R 2 between the statistical area of spring maize and the extracted area in all years was 0.81. The average cumulative ET a of spring maize in the growing season in Jilin Province from 2001 to 2020 ranged from 475 to 544 mm, with a multi-year mean of 517 mm and a standard deviation mean of 43 mm. The highest ET a was in July with a mean of 132 mm. On the spatial scale, the multi-year mean ET a ranged from 391 to 693 mm, and the CV was mostly between 2% and 18%. The cities with the largest multi-year average ET a values were Tonghua, Jilin, and Liaoyuan, and all of them had ET a over 530 mm. The multi-year average ET a of spring maize growing season in Jilin Province in this study was 517 mm, which was similar to the results of other related studies in Jilin Province, and the values were within a reasonable range, indicating that the PySEBS software could estimate the ET a values well. Therefore, PySEBS software can have useful applications in the spatio-temporal patterns of water consumption and irrigation demand in regional scales and associated water saving scenarios.