A Practical Remote Sensing Monitoring Framework for Late Frost Damage in Wine Grapes Using Multi-Source Satellite Data

: Late frost damage is one of the main meteorological disasters that affect the growth of wine grapes in spring, causing a decline in wine grapes quality and a reduction in yield in Northwest China. At present, remote sensing technology has been widely used in the ﬁeld of crop meteorological disasters monitoring and loss assessments, but little research has been carried out on late frost damage in wine grapes. To monitor the impact of late frost in wine grapes accurately and quickly, in this research, we selected the Ningxia planting area as the study area. A practical framework of late frost damage on wine grapes by integrating visible, near-infrared, and thermal infrared satellite data is proposed. This framework includes: (1) Wine grape planting area extraction using Gaofen-1 (GF-1), Landsat-8, and Sentinel-2 based on optimal feature selection and Random Forest (RF) algorithm; (2) retrieval of the land surface temperature (LST) using Landsat-8 thermal infrared data; (3) data fusion using Landsat-8 LST and MODIS LST for a high spatiotemporal resolution of LST with the Enhanced Spatial and Temporal Adaptive Reﬂectance Fusion Model (ESTARFM); (4) the estimation of daily minimum air temperature ( T min ) using downscaled LST and meteorological station data; (5) monitoring and evaluation of the degree of late frost damage in wine grapes in April 2020 by combining satellite-derived data and late frost indicators. The results show that the total area of wine grapes extracted in Ningxia was about 39,837 ha. The overall accuracy was 90.47%, the producer’s accuracy was 91.09%, and the user’s accuracy was 90.22%. The root mean square (RMSE) and the coefﬁcient of determination (R 2 ) of the T min estimation model were 1.67 °C and 0.91, respectively. About 41.12% of the vineyards suffered severe late frost damage, and the total affected area was about 16,381 ha during 20–25 April 2020. This suggests the satellite data can accurately monitor late frost damage in wine grapes by mapping the wine grape area and estimating T min . The results can help farmers to take remedial measures to reduce late frost damage in wine grapes, and provide an objective evaluation of late frost damage insurance claims for wine grapes. With the increasing weather extremes, this study has an important reference value for standardized global wine grape management and food security planning.


Introduction
Grape is one of the most productive fruits with a wide range of cultivation around the world, 80% of which are wine grapes. In evaluations of high-quality wine grape ecological regions, the Ningxia planting area has been classified as one of the best wine grape ecology districts in the world [1]. However, due to its location in the inland northwest, which has frequent cold air activity and temperature fluctuations in spring [2], this district is extremely vulnerable to low temperatures and late frost damages. Late frost damage seriously restricts the development of the local wine grape industry, yield, and quality stability [3]. Traditionally, monitoring and assessment of late frost damage in wine grapes have been based on meteorological stations and field investigations [4,5], but this progress is time-consuming, labor-intensive, and very expensive. Large-area and full-coverage remote sensing technology provide the promising potential for the monitoring and assessment of late frost damage as the remote sensing data can be used to extract the crop planting area [6-10], estimate the air and soil temperature [11,12], monitor the crop growth [13,14], assess the agro-meteorological disasters [15][16][17], and predict crop yields [18][19][20].
The wine grape planting area is the essential information for accurately monitoring and assessing late frost damage in wine grapes. Satellite data have been widely used to extract the planting area of crops, such as paddy rice [6, 10,21,22], corn [7], wheat [23], and soybean [24] with high to medium spatial resolution [25,26]. However, little research has been carried out on extracting the wine grape planting area using high and medium spatial resolution satellite data, and the few studies mainly focused on small-scale research, such as vineyards [27] or at county scales [28]. From using single images [29,30] to multi-source satellite images [22,31,32], the combination of multiple types of data and features, such as sensitive bands, texture features, and vegetation index, can further provide comprehensive information and accurate extraction [23,33]. Several machine learning methods have been used in planting area extraction, such as maximum likelihood [24], support vector machine [34], and Random Forest (RF) [35,36]. The RF model [37,38] has also been widely used to optimize features and remove redundant information. Ningxia was selected as the study area in this paper as it is one of the most important cultivation areas of wine grapes.
Land surface temperature (LST) is an important parameter for late frost monitoring and evaluation [39]. The representative algorithms for LST derived from satellite data mainly include the Jiménez-Muñoz single-channel algorithm [40], the Qin mono-window algorithm [41], the Rozenstein split-window algorithm [42], and the Sobrino radiation transfer equation method [43], etc. These algorithms have widespread applicability [44], and are relatively mature. The spatial resolution of Landsat-8 thermal infrared data is 100 m, but its revisit period is 16 d. MODIS products can provide daytime and nighttime LST from the Terra and Aqua satellite, but their spatial resolution is 1 km. As the area of most wine grape plantations in the study area is less than 1 km 2 , the higher the spatial and temporal resolution, the better for late frost damage monitoring and assessment. Data fusion was used to produce data with high spatial and temporal resolution. Scholars have proposed a variety of downscaling methods for LST, including the Disaggregation Procedure for Radiometric Surface Temperature (DisTrad) [45], the Algorithm for Sharpening Thermal Imagery (TsHARP) [46], the High resolution Urban Thermal Sharpener (HUTS) [47], etc. The spatio-temporal fusion method has been a research hotspot in recent years, and scholars have carried out extensive research on this issue. The Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) [48], the Spatial Temporal Adaptive Algorithm for Mapping Reflectance Change (STAARCH) [49], the Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM) [50], and other fusion methods [51][52][53] were proposed accordingly. In comparison, the ESTARFM takes the heterogeneity of the pixels into account by introducing the conversion coefficients between pure and mixed pixels, which has high fusion accuracy in areas of high heterogeneity [54].
Air temperature is a physical parameter that is often used as the distinguishing factor for low-temperature freezing damage in crops. Since the daily minimum air temperature (T min ) is one of the important indicators of agro-meteorological disasters, LST needs to be converted to T min . Obtaining continuous temporal and spatial daily T min images is the research focus of damage monitoring. Although the observational air temperature of traditional ground meteorological stations usually has high accuracy, due to the limited quantity and uneven distribution of meteorological stations, they cannot meet the demand for spatial continuity. Common remote sensing estimation methods for air temperature include statistical methods, the Temperature-Vegetation Index (TVX) [55], the atmospheric temperature profile method, and the energy balance method, etc. In recent years, based on the correlations [11,56] between air temperature and LST, the method of estimating air temperature by using full coverage LST has been recognized widely.
As current wine grape late frost damage monitoring is relatively sparse distributed, we proposed a practical remote sensing monitoring framework to monitor late frost damage in wine grapes and estimate the damaged area at a large scale. Specifically, this framework includes extraction of the wine grape planting area using Gaofen-1 (GF-1), Landsat-8, and Sentinel-2 satellite data; retrieval of LST using the Landsat-8 satellite data; fusion of Landsat LST and MODIS LST data; and estimation of daily T min . Ultimately, the spatial distribution and the degree of late frost damage in wine grapes were realized according to the disaster indicators, the results of extracting grape planting area, and the daily T min estimation. (Tmin) is one of the important indicators of agro-meteorological disasters, LST needs to be converted to Tmin. Obtaining continuous temporal and spatial daily Tmin images is the research focus of damage monitoring. Although the observational air temperature of traditional ground meteorological stations usually has high accuracy, due to the limited quantity and uneven distribution of meteorological stations, they cannot meet the demand for spatial continuity. Common remote sensing estimation methods for air temperature include statistical methods, the Temperature-Vegetation Index (TVX) [55], the atmospheric temperature profile method, and the energy balance method, etc. In recent years, based on the correlations [11,56] between air temperature and LST, the method of estimating air temperature by using full coverage LST has been recognized widely. As current wine grape late frost damage monitoring is relatively sparse distributed, we proposed a practical remote sensing monitoring framework to monitor late frost damage in wine grapes and estimate the damaged area at a large scale. Specifically, this framework includes extraction of the wine grape planting area using Gaofen-1 (GF-1), Landsat-8, and Sentinel-2 satellite data; retrieval of LST using the Landsat-8 satellite data; fusion of Landsat LST and MODIS LST data; and estimation of daily Tmin. Ultimately, the spatial distribution and the degree of late frost damage in wine grapes were realized according to the disaster indicators, the results of extracting grape planting area, and the daily Tmin estimation.

Study Area
Located in the Ningxia Hui Autonomous Region of Northwest China, the study area along the Helan Mountain extends from 37.47°N to 39.42°N and 105.62°E to 106.72°E (Figure 1). The Ningxia wine grape planting area includes four growing districts: Shizuishan, Yinchuan, Qingtongxia, and Hongsipu. It lies the temperate continental climate zone, with an average annual temperature range of 8-10 °C and an average annual precipitation range of approximately 130-250 mm. The annual sunshine hours are up to 3000 h, and the frost-free period is about 150 days.  With moderate heat and abundant sunshine, the study area has a typical continental climate, which is very conducive to the accumulation of sugar and the transformation of pigment, making it one of the major wine grape planting areas of China and worldrecognized "golden area" and "Chinese Bordeaux" for wine grape cultivation. Through

Satellite Data and Preprocessing
The satellite data used in this paper can be divided into two parts: (1) GF-1, Landsat-8, and Sentinel-2 data for extracting wine grape planting area; and (2) Landsat-8 thermal infrared data and the MODIS product set (LST and EVI) for LST retrieval and data fusion.
Sentinel-2 images are provided freely by the European Space Agency (ESA) scientific data hub portal (https://scihub.copernicus.eu/, accessed on 31 May 2020). Landsat-8 images are freely available for users to download from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, accessed on 1 June 2020). GF-1 satellite images are provided freely by the China Centre for Resources Satellite Data and Application (CRESDA) (http://www.cresda.com, accessed on 12 April 2020). GF-1 is the first satellite of a series of high-resolution earth observation systems from China, launched in April 2013. GF-1 Wide Field of View (WFV) images have an overlapping swath of 830 km and contain four bands (blue, green, red, and near-infrared) with a spatial resolution of 16 m. MODIS products are freely available to users via the National Aeronautics and Space Administration (NASA) satellite data website and can be downloaded from https://ladsweb.modaps.eosdis.nasa.gov/search/ (accessed on 3 July 2020).
To obtain enough images to cover the main phenological stages of wine grape in the study area, 10 cloud-free multi-source satellite images from April 2019 to April 2020 were used to extract the planting area, and they included two GF-1 WFV images with a spatial resolution of 16 m, four Landsat-8 Operational Land Imager (OLI) images with a spatial resolution of 30 m and four Sentinel-2 Multispectral Instrument (MSI) images with a spatial resolution of 10 m (Table 1).

Figure 2.
The seasonal pattern of the Normalized Difference Vegetation Index (NDVI) corresponding to different wine grape phenology.

Satellite Data and Preprocessing
The satellite data used in this paper can be divided into two parts: (1) GF-1, Landsat-8, and Sentinel-2 data for extracting wine grape planting area; and (2) Landsat-8 thermal infrared data and the MODIS product set (LST and EVI) for LST retrieval and data fusion.
Sentinel-2 images are provided freely by the European Space Agency (ESA) scientific data hub portal (https://scihub.copernicus.eu/, accessed on 31 May 2020). Landsat-8 images are freely available for users to download from the United States Geological Survey (USGS) website (https://earthexplorer.usgs.gov/, accessed on 1 June 2020). GF-1 satellite images are provided freely by the China Centre for Resources Satellite Data and Application (CRESDA) (http://www.cresda.com, accessed on 12 April 2020). GF-1 is the first satellite of a series of high-resolution earth observation systems from China, launched in April 2013. GF-1 Wide Field of View (WFV) images have an overlapping swath of 830 km and contain four bands (blue, green, red, and near-infrared) with a spatial resolution of 16 m. MODIS products are freely available to users via the National Aeronautics and Space Administration (NASA) satellite data website and can be downloaded from https://ladsweb.modaps.eosdis.nasa.gov/search/ (accessed on 3 July 2020).
To obtain enough images to cover the main phenological stages of wine grape in the study area, 10 cloud-free multi-source satellite images from April 2019 to April 2020 were used to extract the planting area, and they included two GF-1 WFV images with a spatial resolution of 16 m, four Landsat-8 Operational Land Imager (OLI) images with a spatial resolution of 30 m and four Sentinel-2 Multispectral Instrument (MSI) images with a spatial resolution of 10 m (Table 1). Three cloud-free Landsat-8 Thermal Infrared Sensor (TIRS) images were acquired on 22 March, 23 April, and 9 May 2020 at a spatial resolution of 100 m. Thirty Aqua MODIS images of nighttime LST (MYD11A1) from 1 April to 30 April 2020; three images of daytime LST (MOD11A1) on 22 March, 23 April, and 9 May 2020; and seven MODIS images of EVI products (MOD13A2.A2020081/A2020097/A2020113/A2020129, MYD13A2.A2020089/A2020105/A2020121) were acquired with a spatial resolution of 1 km.
The FLAASH module [57] of ENVI was used for atmospheric correction of the GF-1 and Landsat-8 images, and the sentinel-2 images were preprocessed by Sen2Cor (http://step.esa. int/main/third-party-plugins-2/sen2cor/, accessed on 5 June 2020). and used as reference images to geometrically correct the GF-1 images from the same day. The GF-1 and Sentinel-2 images were all resampled to 30 m spatial resolution. The digital elevation model (DEM) used in preprocessing was Shuttle Radar Topography Mission (SRTM) DEM. All MODIS products downloaded used the MODIS Reprojection Tool (MRT) provided by NASA for reprojection, data layer extraction, and format conversion (HDF to TIF).

Meteorological Data
In this study, daily T min data from 1 April to 30 April 2020 were collected from 180 meteorological stations (Figure 1) in the study area. The LST values used as validation data on April 2020 were acquired from the Real-Time Product Dataset of the China Meteorological Administration Land Data Assimilation System (CLDAS-V2.0), which were downloaded from China meteorological data service center (CMDC) (http://data.cma.cn/, accessed on 15 October 2020).

Field Survey and Sample Set
The field survey was carried out from July 30 to August 3 in 2019. We visited 10 vineyards in Yinchuan, Qingtongxia, and Hongsipu to collect information about the wine grape distribution, grape species, and cultivation measures in the study area. Based on the field survey and visual interpretations of Google Earth, more than 2000 samples (consisting of over 150,000 pixels) of wine grapes, farmland, woodland, meadow, desert steppe, desert, building, and water were selected as training and validating samples ( Figure 3). In this study, 60% were used as training samples; the remaining 40% were used to validate and evaluate the extraction results. Remote Sens. 2021, 13, 3231 6 of 23

Methods
The framework proposed in this study is presented in Figure 4. Firstly, the wine grape planting area was extracted by optimal feature selection and the RF method. Secondly, Landsat-8 LST (using two retrieval algorithms) and MODIS LST (after cloudy gapfilling) were fused to obtain a seamless daily high spatial resolution of LST (100 m). Thirdly, fully covered Tmin was estimated based on the correlations between downscaled LST and meteorological station Tmin data. Finally, by using daily Tmin images of the study area, wine grape planting area image, and the evaluation indicator of wine grape late frost damage, we obtained a map of the late frost damaged area of wine grapes in the study area in April 2020.

Extraction of the Wine Grape Planting Area Using RF
RF is a machine learning algorithm that uses multiple decision trees to train and predict samples. It was first proposed by Breiman [58]. This method can effectively solve the over-fitting phenomenon by double random sampling of the training samples and selected feature variables and constructing multiple decision trees through the idea of integrated learning. Each decision tree grows to the maximum without any pruning. Through internal evaluation, the unbiased error is estimated in the process of RF generation, and the final prediction and classification results depend on the mean prediction value and the voting results of multiple decision trees.
In this study, the RF algorithm and Gini index were used to evaluate the importance of features and determine the optimal feature subset according to the voting mechanism [59]. Gini index evaluates the classification contributions of the feature variables in each tree of RF, normalizes all the feature importance scores, and compares the feature importance weights. The Python language and the Scikit-learn library were used to implement the feature optimization progress [60].

Methods
The framework proposed in this study is presented in Figure 4. Firstly, the wine grape planting area was extracted by optimal feature selection and the RF method. Secondly, Landsat-8 LST (using two retrieval algorithms) and MODIS LST (after cloudy gap-filling) were fused to obtain a seamless daily high spatial resolution of LST (100 m). Thirdly, fully covered T min was estimated based on the correlations between downscaled LST and meteorological station T min data. Finally, by using daily T min images of the study area, wine grape planting area image, and the evaluation indicator of wine grape late frost damage, we obtained a map of the late frost damaged area of wine grapes in the study area in April 2020.

Extraction of the Wine Grape Planting Area Using RF
RF is a machine learning algorithm that uses multiple decision trees to train and predict samples. It was first proposed by Breiman [58]. This method can effectively solve the over-fitting phenomenon by double random sampling of the training samples and selected feature variables and constructing multiple decision trees through the idea of integrated learning. Each decision tree grows to the maximum without any pruning. Through internal evaluation, the unbiased error is estimated in the process of RF generation, and the final prediction and classification results depend on the mean prediction value and the voting results of multiple decision trees.
In this study, the RF algorithm and Gini index were used to evaluate the importance of features and determine the optimal feature subset according to the voting mechanism [59]. Gini index evaluates the classification contributions of the feature variables in each tree of RF, normalizes all the feature importance scores, and compares the feature importance weights. The Python language and the Scikit-learn library were used to implement the feature optimization progress [60]. Different crops show different spectral characteristics, due to their different pigment and moisture content and different canopy structures [61]. Therefore, spectral features can be used for crop recognition. In total, 52 spectral features were obtained from 10 multisource satellite images, including eight bands (Blue, Green, Red, Near-infrared) of two GF-1 images, 16 bands (blue, green, red, and near-infrared) of four Landsat-8 images, and 28 bands (blue, green, red, near-infrared, and three red-edge bands) of four Sentinel-2 images.
The vegetation indices were calculated from different spectral reflectance for vegetation discrimination and quantitative retrieval of parameters. NDVI is one of the most important indices for deriving crop phenological parameters and distinguishes different types of ground objects, which can reflect the difference between the vegetation and soil and is suitable for dynamic monitoring during early growth stages. As a common indicator parameter of green plants, the Normalized Difference Red Edge Index (NDRE) has a high sensitivity for vegetation recognition and classification, as it can take advantage of the red-edge bands in the field of agricultural remote sensing monitoring [62]; it is widely used in agriculture and forest monitoring. In this study, NDVI values were calculated from the 10 multi-source satellite images. Sentinel-2 contains three red-edge bands, so 12 NDRE features were calculated from the four Sentinel-2 images. NDRE features were not calculated from the GF-1 and Landsat-8 images, due to their lack of red-edge bands.
Texture features measures the surface roughness and uniformity of image pixel by the similarity of spatial structure, reflecting the close relationship between texture information and gray-scale [63], which has become one of the major factors in image classification [64]. The Principal Component Analysis (PCA) transformation was implemented on the spectral bands of 10 images, and only the first PCA, which had the highest amount of information, was selected to calculate five texture features based on the Gray Level Cooccurrence Matrix (GLCM). The five texture features were variance, contrast, correlation, angular second moment, and entropy.
In total, 124 features were obtained from the NDVI, NDRE, spectra, and texture data, from which the feature set was selected, as listed in Table 2. Different crops show different spectral characteristics, due to their different pigment and moisture content and different canopy structures [61]. Therefore, spectral features can be used for crop recognition. In total, 52 spectral features were obtained from 10 multisource satellite images, including eight bands (Blue, Green, Red, Near-infrared) of two GF-1 images, 16 bands (blue, green, red, and near-infrared) of four Landsat-8 images, and 28 bands (blue, green, red, near-infrared, and three red-edge bands) of four Sentinel-2 images.
The vegetation indices were calculated from different spectral reflectance for vegetation discrimination and quantitative retrieval of parameters. NDVI is one of the most important indices for deriving crop phenological parameters and distinguishes different types of ground objects, which can reflect the difference between the vegetation and soil and is suitable for dynamic monitoring during early growth stages. As a common indicator parameter of green plants, the Normalized Difference Red Edge Index (NDRE) has a high sensitivity for vegetation recognition and classification, as it can take advantage of the red-edge bands in the field of agricultural remote sensing monitoring [62]; it is widely used in agriculture and forest monitoring. In this study, NDVI values were calculated from the 10 multi-source satellite images. Sentinel-2 contains three red-edge bands, so 12 NDRE features were calculated from the four Sentinel-2 images. NDRE features were not calculated from the GF-1 and Landsat-8 images, due to their lack of red-edge bands.
Texture features measures the surface roughness and uniformity of image pixel by the similarity of spatial structure, reflecting the close relationship between texture information and gray-scale [63], which has become one of the major factors in image classification [64]. The Principal Component Analysis (PCA) transformation was implemented on the spectral bands of 10 images, and only the first PCA, which had the highest amount of information, was selected to calculate five texture features based on the Gray Level Co-occurrence Matrix (GLCM). The five texture features were variance, contrast, correlation, angular second moment, and entropy.

Data Fusion of LST
The mono-window algorithm (LST_MW) [65], and the nonlinear split-window algorithm [66] were utilized to calculate LST using Landsat-8 TIRS on 22 March, 23 April, and 9 May 2020, respectively. The LST_MW constructs the calculation model of LST, atmospheric transmittance, brightness temperature, and average atmospheric temperature based on the thermal radiation transfer equation and Landsat-8 TIRS 10. The split-window algorithm (LST_SW) used the atmospheric absorption difference between the two thermal infrared channels to correct the atmospheric impact, establishing the relationship between the LST and the brightness temperature.
The MODIS LST products with a temporal resolution of 1 day provided an excellent source for the temporal continuity analysis of LST. Since wine grape late frost damage usually occurs during 03:00-04:00, which is close to the observation time of the Aqua MODIS nighttime LST (MYD11A1), MYD11A1 products were selected as the data source for LST with a high temporal resolution of 1 day and a low spatial resolution of 1 km. However, the quality of pixels largely depends on the cloud coverage, so cloudy weather may lead to poor pixel quality. Since the Helan Mountain lies to the west of the study area, which is affected by snow accumulation in early spring, the mountainous area often has many low-quality pixels. In this study, the spatiotemporal linear interpolation method was used to find similar properties between the effective pixels and invalid pixels (the pixels with poor quality were removed) in the adjacent time and space; then, the missing pixels were interpolated by calculating the conversion relationship [67].
Assuming that the LST of similar pixels can be known at time t 0 and t 1 , then the linear relationship between two LST images can be defined by statistical linear regression.
where LST t 0 is the LST at t 0 , and LST t 1 is the LST at t 1 . If the LST of the pixel at t 0 is unknown, but that at t 1 is known, then the LST of the pixel at t 0 can be obtained through the linear relationship above.
Given the cloud cover in the study area, LST images with cloud-free pixel coverage of more than 70% were selected as the reference image, and the image with cloud-free pixel coverage less than 70% as interpolated images. We then established the linear model linking the cloud-free pixels of the interpolated images, cloud-free pixels of the reference images, DEM and EVI. Due to the different measurements of LST, DEM, and EVI, the three values were normalized before establishing the linear relationship.
where LST int is the array of interpolated images, LST re f is the cloud-free pixel array of reference image, and a, b, c, and d are regression coefficients. With this equation, the LST pixels of cloud cover can be calculated using the reference images, DEM and EVI. For the missing pixels after filling, the CLDAS-V2.0 0 cm LST data set was used for further filling. The downloaded CLDAS-V2.0 0 cm LST can be used as supplementary data after format conversion, coordinate system conversion, resampling, clipping, and other processing.
To obtain LST with a high spatiotemporal resolution, ESTARFM was used to fuse the Landsat-8 LST and MODIS LST. ESTARFM is a data fusion algorithm based on a moving window [54,68]. It is based on the premise that the image registration of fine-resolution and coarse-resolution images is accurate, and it is assumed that the reflectivity difference between the images is only caused by system deviation, and there is no significant difference between the images from two periods [50,69]. The satellite images used for downscaling model are listed in Table 3. Assuming that the land cover types and the sensor calibration of the two images are consistent, the relationship between fine-resolution (Landsat-8 LST) and coarse-resolution (MODIS LST) images can be described by a linear model.  (3) and (4): Because of the complexity of the land surface, most of the pixels in the images are mixed pixels, that is, one pixel contains different land cover types. It is difficult to ensure the accuracy of the prediction results by using only a single pixel. Therefore, assuming that the reflectivity of the mixed pixels is a linear combination of that in different land cover types, and that the proportion of land cover types in the two images remains basically unchanged, the reflectivity of the mixed pixel can be expressed as: L(x, y, t p , B) = L(x, y, t l , B) + v(x, y) · M(x, y, t p , B) − M(x, y, t l , B) where v(x, y) is the conversion coefficient of mixed pixels.
According to the principle that the similar neighboring pixels have similar reflectivity, a moving window is used to obtain the information of the neighboring pixels, and the information of the similar pixels is then integrated into the calculation of fine-resolution reflectivity. From the similarity between neighboring pixels, the central pixel values of the predicted images are calculated as follows: where (x w/2 , y w/2 ) is the location of the central pixel, w is the size of the moving window, N is the number of similar pixels, including the central predicted pixel, W i is the weight of the ith similar pixel, V i is the conversion coefficient of the ith similar pixel, and (x i , y i ) is the location of ith similar pixel.
Two images from different dates were used to calculate the value of the center pixel. Landsat-8 LST and MODIS LST images at t m and t n were used to predict the fine-resolution reflectance of date t p , and the weighted combination of the two prediction results was calculated to obtain a more accurate value. The temporal weight was calculated by the magnitude of change in the MODIS values between t l (l = m, n) and t p , and Equation (8) is as follows: Then the final predicted fine-resolution reflectance of the center pixel at t p is as follows: L(x w/2 , y w/2 , t p , B) = T m · L m (x w/2 , y w/2 , t m , B) + T n · L n (x w/2 , y w/2 , t n , B) where L m (x w/2 , y w/2 , t m , B) and L n (x w/2 , y w/2 , t n , B) are fine-resolution pixel value, respectively, at t m and t n , and T m and T n are the temporal weight of t m and t n .

Estimation Model of Daily Minimum Air Temperature
The daily T min was estimated by the 30 images of downscaled LST (LST E ) in April 2020 with the spatial resolution of 100 m. The observation time for MODIS LST was 1:30 a.m. every day, and the late frost of wine grapes often occurs during 03:00-04:00, and the daily T min in the study area occurs during approximately 05:00-07:00 [4]. The research showed that temperature at these key times had a strong linear correlation, so we can monitor the late frost damage by building a linear estimation model. Because LST E was calculated from MODIS LST, LST E is considered to have a similar property with MODIS LST and also has a good linear relationship with the daily T min . Therefore, the study constructed a linear model of the daily T min and LST E to estimate the daily T min using Equation (10). The daily T min in April 2020, as input variable in the model, was derived from meteorological stations.
T min = a · LST E + b where a and b are the regression coefficients of ordinary least squares regression.

Mapping the Late Frost Damaged Area
The evaluation indicators of wine grape late frost damage from the Yinchuan Meteorological Bureau (Figure 4), were proved to be reliable through experiments and verification. The occurrence of wine grape late frost damage is related to air temperature, precipitation, and soil conditions, but the low air temperature is the most important factor. In this study, we classified the degree of late frost damage according to the evaluation indicator and the daily T min of each wine grape pixel.
As described in Sections 2.5.1 and 2.5.3, daily minimum air temperature images (DMTI), wine grape planting area image (WGPAI), and the evaluation index of wine grape late frost damage (LFDEI) were obtained and integrated to derive the late frost damage images (LFDI) [17]. Equation (11) is as follows: The spatial distribution information of late frost damage was then further combined with the estimated area of wine grapes to estimate the late frost damaged area in the Ningxia planting area in April 2020.

Extraction of the Wine Grape Planting Area
In total, 124 features, including spectral reflectance, vegetation indices, and texture features, were calculated. Rankings of the feature importance and performance are shown in Figure 5, and the changes in the overall accuracy (OA) and kappa coefficient with the number of feature variables are shown in Figure 6.

Mapping the Late Frost Damaged Area
The evaluation indicators of wine grape late frost damage from the Yinchuan Meteorological Bureau (Figure 4), were proved to be reliable through experiments and verification. The occurrence of wine grape late frost damage is related to air temperature, precipitation, and soil conditions, but the low air temperature is the most important factor. In this study, we classified the degree of late frost damage according to the evaluation indicator and the daily Tmin of each wine grape pixel.
As described in Sections 2.5.1 and 2.5.3, daily minimum air temperature images (DMTI), wine grape planting area image (WGPAI), and the evaluation index of wine grape late frost damage (LFDEI) were obtained and integrated to derive the late frost damage images (LFDI) [17]. Equation (11) is as follows: The spatial distribution information of late frost damage was then further combined with the estimated area of wine grapes to estimate the late frost damaged area in the Ningxia planting area in April 2020.

Extraction of the Wine Grape Planting Area
In total, 124 features, including spectral reflectance, vegetation indices, and texture features, were calculated. Rankings of the feature importance and performance are shown in Figure 5, and the changes in the overall accuracy (OA) and kappa coefficient with the number of feature variables are shown in Figure 6.   The OA and the kappa coefficient first increased rapidly with the number of feature variables ( Figure 5), then, increased slowly and fluctuated up to 13 features. When the number of features reached 53, the total accuracy and kappa coefficient value reached the maximum. Therefore, the top 53 features were used as the optimal feature subset to extract wine grape planting area. Among the top 53 features (Figure 6), which included 21 spectral features, 20 texture features, and 12 vegetation index features (6 NDVI features and 6 NDRE features), NDVI 3 had the highest feature importance score, and NIR 10 had the lowest.
Based on the testing samples, the extraction results were assessed using the confusion matrix with OA, producer accuracy (PA), user accuracy (UA), and the kappa coefficient of the optimal feature subset in Table 4. The OA and kappa coefficient were 90.47% and 0.89, respectively. The wine grape PA and UA were 91.09% and 90.22%, respectively. The omission and commission error may have been caused by mixed pixels around the vineyard border, due to the limitation of the spatial resolution of the satellite images. On the other hand, wine grapes in new plantations can be classified as desert, due to their low vegetation coverage and similarity to bare soil. The extracted area of wine grape planting regions (Figure 7) was calculated to be about 39,837 ha. The OA and the kappa coefficient first increased rapidly with the number of feature variables ( Figure 5), then, increased slowly and fluctuated up to 13 features. When the number of features reached 53, the total accuracy and kappa coefficient value reached the maximum. Therefore, the top 53 features were used as the optimal feature subset to extract wine grape planting area. Among the top 53 features (Figure 6), which included 21 spectral features, 20 texture features, and 12 vegetation index features (6 NDVI features and 6 NDRE features), NDVI 3 had the highest feature importance score, and NIR 10 had the lowest.
Based on the testing samples, the extraction results were assessed using the confusion matrix with OA, producer accuracy (PA), user accuracy (UA), and the kappa coefficient of the optimal feature subset in Table 4. The OA and kappa coefficient were 90.47% and 0.89, respectively. The wine grape PA and UA were 91.09% and 90.22%, respectively. The omission and commission error may have been caused by mixed pixels around the vineyard border, due to the limitation of the spatial resolution of the satellite images. On the other hand, wine grapes in new plantations can be classified as desert, due to their low vegetation coverage and similarity to bare soil. The extracted area of wine grape planting regions (Figure 7) was calculated to be about 39,837 ha.

LST Retrieval Using Landsat-8 Thermal Infrared Data
The LST_MW and the LST_SW were used to retrieve the LST from Landsat-8 TIRS (TIRS 10: 10.60~11.19 μm; TIRS 11: 11.50~12.51 μm) at a resolution of 100 m. One hundred pixels were randomly selected from every LST_MW and LST_SW image. The verification pixels of three dates obtained by the same algorithm were combined into a group, with a total of 300 pixels in each group. Through a comparison with the corresponding pixel values of CLDAS-V2.0 at 11:00 (the observation time of Landsat-8 is about 11:30 a.m. Beijing time), the accuracy of the two LST retrieval algorithms was analyzed. Figure 8 shows the scatter plot of the LST calculated by the two algorithms versus CLDAS-V2.0 LST.

Wine grape Others
Boundary of wine grape growing district Figure 7. Wine grape map of the study area.

LST Retrieval Using Landsat-8 Thermal Infrared Data
The LST_MW and the LST_SW were used to retrieve the LST from Landsat-8 TIRS (TIRS 10: 10.60~11.19 µm; TIRS 11: 11.50~12.51 µm) at a resolution of 100 m. One hundred pixels were randomly selected from every LST_MW and LST_SW image. The verification pixels of three dates obtained by the same algorithm were combined into a group, with a total of 300 pixels in each group. Through a comparison with the corresponding pixel values of CLDAS-V2.0 at 11:00 (the observation time of Landsat-8 is about 11:30 a.m. Beijing time), the accuracy of the two LST retrieval algorithms was analyzed. Figure 8

LST Retrieval Using Landsat-8 Thermal Infrared Data
The LST_MW and the LST_SW were used to retrieve the LST from Landsat-8 TIRS (TIRS 10: 10.60~11.19 μm; TIRS 11: 11.50~12.51 μm) at a resolution of 100 m. One hundred pixels were randomly selected from every LST_MW and LST_SW image. The verification pixels of three dates obtained by the same algorithm were combined into a group, with a total of 300 pixels in each group. Through a comparison with the corresponding pixel values of CLDAS-V2.0 at 11:00 (the observation time of Landsat-8 is about 11:30 a.m. Beijing time), the accuracy of the two LST retrieval algorithms was analyzed. Figure 8 shows the scatter plot of the LST calculated by the two algorithms versus CLDAS-V2.0 LST.  We revealed that in the scatter plot of the two groups of 300 verification pixels; the root mean square (RMSE) of LST_MW and LST_SW were 2.15 • C and 2.51 • C, respectively.
In addition, the coefficient of determination (R 2 ) of LST_ MW was 0.78, and the R 2 of LST_ SW was 0.76. Based on the above analysis, although there were some differences between the two algorithms in the process of LST retrieval, the two algorithms both produced ideal retrieval results, and the accuracy of the LST_MW algorithm was relatively higher than that of the LST_SW algorithm in this study. Therefore, the LST calculated by the LST_MW algorithm was used for subsequent research.

Cloud Gap-Filling of MODIS LST Data
During April, the valid pixel rates of the original daily MODIS LST images in the study area were 8.77% to 82.95%. The cloud-free pixel coverage on 4, 11, 16, 26, 27, 28, and 30 April were 71.57, 79.46, 79.80, 82.95, 79.04, 71.65, and 80.85, respectively. Therefore, these LST images were selected as the reference image, and the others were interpolated images. The multivariate regression model was built by using the cloud-free pixels of the interpolated images as the dependent variables and the cloud-free pixels of the reference image, DEM and EVI, as the independent variables. The model increased the valid pixel rates of all interpolated images. The effect of pixel interpolation in daily MODIS LST images is shown in Figure 9. After interpolation, the valid pixel rates of all interpolated images were 71.57% to 86.40%.
Remote Sens. 2021, 13, 3231 14 of 23 We revealed that in the scatter plot of the two groups of 300 verification pixels; the root mean square (RMSE) of LST_MW and LST_SW were 2.15 °C and 2.51 °C, respectively. In addition, the coefficient of determination (R 2 ) of LST_ MW was 0.78, and the R 2 of LST_ SW was 0.76. Based on the above analysis, although there were some differences between the two algorithms in the process of LST retrieval, the two algorithms both produced ideal retrieval results, and the accuracy of the LST_MW algorithm was relatively higher than that of the LST_SW algorithm in this study. Therefore, the LST calculated by the LST_MW algorithm was used for subsequent research.

Cloud Gap-Filling of MODIS LST Data
During April, the valid pixel rates of the original daily MODIS LST images in the study area were 8.77% to 82.95%. The cloud-free pixel coverage on 4, 11, 16, 26, 27, 28, and 30 April were 71.57, 79.46, 79.80, 82.95, 79.04, 71.65, and 80.85, respectively. Therefore, these LST images were selected as the reference image, and the others were interpolated images. The multivariate regression model was built by using the cloud-free pixels of the interpolated images as the dependent variables and the cloud-free pixels of the reference image, DEM and EVI, as the independent variables. The model increased the valid pixel rates of all interpolated images. The effect of pixel interpolation in daily MODIS LST images is shown in Figure 9. After interpolation, the valid pixel rates of all interpolated images were 71.57% to 86.40%. The MODIS LST images in April were compared with the pixels of CLDAS-V2.0 LST as follows: All 30 images were taken as a group, randomly choose 50 interpolated pixels in each image; that is, 1500 random pixels were selected for the validation. We validated the cloud gap-filling methods ( Figure 10); overall, the R 2 between the interpolated MODIS LST and CLDAS-V2.0 LST was 0.87 with RMSE of 1.79 °C, which indicates a strong correlation between MODIS LST and CLDAS-V2.0 LST. Furthermore, for the pixels still vacant after interpolation, CLDAS-V2.0 0 cm LST real-time products were used to fill them indirectly, then the valid pixel rates of all images were 100%. The final nighttime LST images with full spatial coverage at 1 km resolution in the study area were obtained and used as the data source for the data fusion stage of the research. The MODIS LST images in April were compared with the pixels of CLDAS-V2.0 LST as follows: All 30 images were taken as a group, randomly choose 50 interpolated pixels in each image; that is, 1500 random pixels were selected for the validation. We validated the cloud gap-filling methods ( Figure 10); overall, the R 2 between the interpolated MODIS LST and CLDAS-V2.0 LST was 0.87 with RMSE of 1.79 • C, which indicates a strong correlation between MODIS LST and CLDAS-V2.0 LST. Furthermore, for the pixels still vacant after interpolation, CLDAS-V2.0 0 cm LST real-time products were used to fill them indirectly, then the valid pixel rates of all images were 100%. The final nighttime LST images with full spatial coverage at 1 km resolution in the study area were obtained and used as the data source for the data fusion stage of the research.

Data Fusion of Landsat-8 LST and MODIS LST Using the ESTARFM Method
Temporal resolution changes before and after data fusion are shown in Figure 11a. Taking four 13 km × 13 km local areas in Figure 11b as an example, the spatial resolution results of data fusion can be visualized by comparing the differences in the LST image pixels before and after data fusion for 24 April. Figure 11 shows the data fusion result of LST both in temporal and spatial resolution. In Figure 11a, only three images of Landsat-8 LST were used before data fusion. After fusion by ESTARFM, a complete data set containing 30 images of LST in April was obtained. In Figure 11b, compared with the LST images with a spatial resolution of 1 km, the downscaled images have more detailed information and a clear description of the LST pixels. With regard to the phenomenon of local blurring and mutation at the boundary of different types of ground objects, this study analyzed the reasons, which may be the degradation of the downscaling quality caused by the high degree of fragmentation of mixed pixels. Overall, downscaled LST had a high correlation with MODIS LST and could be used to construct the model of the following Tmin estimation. (a)

Data Fusion of Landsat-8 LST and MODIS LST Using the ESTARFM Method
Temporal resolution changes before and after data fusion are shown in Figure 11a. Taking four 13 km × 13 km local areas in Figure 11b as an example, the spatial resolution results of data fusion can be visualized by comparing the differences in the LST image pixels before and after data fusion for 24 April. Figure 11 shows the data fusion result of LST both in temporal and spatial resolution. In Figure 11a, only three images of Landsat-8 LST were used before data fusion. After fusion by ESTARFM, a complete data set containing 30 images of LST in April was obtained. In Figure 11b, compared with the LST images with a spatial resolution of 1 km, the downscaled images have more detailed information and a clear description of the LST pixels. With regard to the phenomenon of local blurring and mutation at the boundary of different types of ground objects, this study analyzed the reasons, which may be the degradation of the downscaling quality caused by the high degree of fragmentation of mixed pixels. Overall, downscaled LST had a high correlation with MODIS LST and could be used to construct the model of the following T min estimation.

Data Fusion of Landsat-8 LST and MODIS LST Using the ESTARFM Method
Temporal resolution changes before and after data fusion are shown in Figure 11a. Taking four 13 km × 13 km local areas in Figure 11b as an example, the spatial resolution results of data fusion can be visualized by comparing the differences in the LST image pixels before and after data fusion for 24 April. Figure 11 shows the data fusion result of LST both in temporal and spatial resolution. In Figure 11a, only three images of Landsat-8 LST were used before data fusion. After fusion by ESTARFM, a complete data set containing 30 images of LST in April was obtained. In Figure 11b, compared with the LST images with a spatial resolution of 1 km, the downscaled images have more detailed information and a clear description of the LST pixels. With regard to the phenomenon of local blurring and mutation at the boundary of different types of ground objects, this study analyzed the reasons, which may be the degradation of the downscaling quality caused by the high degree of fragmentation of mixed pixels. Overall, downscaled LST had a high correlation with MODIS LST and could be used to construct the model of the following Tmin estimation.

Calibration and Validation of Daily Minimum Air Temperature Estimation Using the Downscaled LST Data
A linear regression model was constructed between Tmin measured from 180 meteorological stations and the fused LSTE derived from the Aqua nighttime MODIS LST and Landsat LST. Two-thirds of the meteorological data were used for modeling and one-third for validation.
The regression model between Tmin and LSTE is formatted in Equation (12).
RMSE was 3.07 °C; the R 2 and Bias for the estimation model of daily minimum air temperature were 0.85 and 0.008, respectively.
Fully covered Tmin of the study area can be calculated by Equation (12). In Figure 12b, the scattered points are distributed near the 1:1 line. The RMSE and R 2 of the estimated Tmin and measured Tmin were 1.67 °C and 0.91, respectively. The low error value and good correlation indicated that the estimation model had high applicability for the estimation of Tmin in the study area, and the estimated results can be used for remote sensing monitoring and grade evaluation of wine grape late frost damage.

Calibration and Validation of Daily Minimum Air Temperature Estimation Using the Downscaled LST Data
A linear regression model was constructed between T min measured from 180 meteorological stations and the fused LST E derived from the Aqua nighttime MODIS LST and Landsat LST. Two-thirds of the meteorological data were used for modeling and one-third for validation.
The regression model between T min and LST E is formatted in Equation (12).
RMSE was 3.07 • C; the R 2 and Bias for the estimation model of daily minimum air temperature were 0.85 and 0.008, respectively.
Fully covered T min of the study area can be calculated by Equation (12). In Figure 12b, the scattered points are distributed near the 1:1 line. The RMSE and R 2 of the estimated T min and measured T min were 1.67 • C and 0.91, respectively. The low error value and good correlation indicated that the estimation model had high applicability for the estimation of T min in the study area, and the estimated results can be used for remote sensing monitoring and grade evaluation of wine grape late frost damage.

Calibration and Validation of Daily Minimum Air Temperature Estimation Using the Downscaled LST Data
A linear regression model was constructed between Tmin measured from 180 meteorological stations and the fused LSTE derived from the Aqua nighttime MODIS LST and Landsat LST. Two-thirds of the meteorological data were used for modeling and one-third for validation.
The regression model between Tmin and LSTE is formatted in Equation (12).
RMSE was 3.07 °C; the R 2 and Bias for the estimation model of daily minimum air temperature were 0.85 and 0.008, respectively.
Fully covered Tmin of the study area can be calculated by Equation (12). In Figure 12b, the scattered points are distributed near the 1:1 line. The RMSE and R 2 of the estimated Tmin and measured Tmin were 1.67 °C and 0.91, respectively. The low error value and good correlation indicated that the estimation model had high applicability for the estimation of Tmin in the study area, and the estimated results can be used for remote sensing monitoring and grade evaluation of wine grape late frost damage.

Mapping of the Late Frost Damaged Area
The distribution map of the wine grape planting area with a spatial resolution of 30m was resampled to 100 m, which was consistent with the spatial resolution of T min . The spatial distribution of wine grapes and full spatial coverage of T min at 100 m resolution in each planting region in April 2020 were combined to obtain the daily late frost damaged area. The results showed that in April 2020, affected by the very cold air, days with low air temperature appeared frequently and intensively. There were about 14 days when the daily T min was lower than 0 • C in the planting area; that is, there were 14 late frost days, which were 2-4 April, 11-13 April, 16-17 April, and 20-25 April. Among these days, extremely low air temperature appeared on the morning of 24 April, when the T min suddenly dropped to about −6 • C; the T min on the other days of the fourth late frost period were all below 0 • C. Although the T min was lower than 0 • C, most of the wine grapes were not frozen before 20 April because they had just been unearthed from the soil, but had not yet sprouted. This period had strong resistance to low air temperature, which had little effect on the growth of wine grape buds and leaves. Therefore, we monitored the late frost situation during 20-25 April. Figure 13 shows the remote sensing monitoring of late frost in the wine grape planting area in April 2020. The severe damage area is mainly concentrated in the central part of the study area, with less frost damaged in the northern part. Combined with the local field survey, we found that the degree of late frost damage was relatively severe in areas with low or flat terrain. Table 5 shows the affected area during the late frost damage in April 2020. We found that a total of 77.48% of the vineyards suffered from late frost damage, and the area of wine grapes affected by severe late frost damage reached 16,381 ha, accounting for 41.12% of the total planting area. Among the areas, Yinchuan was the most severely damaged area, with about 8501 ha, and the total damaged area of Qingtongxia, and Hongsipu was about 7393 ha. The planting scale of wine grapes in the Shizuishan growing district was relatively small, so the damaged area was significantly smaller than other growing districts. According to the local planting area statistics, the mortality rate of wine grape buds in some wine chateaus with severe freezing damage was 70%; only about 30% of them survived, and vineyards were faced with great loss of production [70]. The results of late frost damage obtained by estimating the T min was in good agreement with the statistics from the agricultural meteorological disaster department.

Discussion
In this paper, a practical remote sensing monitoring framework of wine grape late frost damage was proposed using multi-source satellite data and auxiliary ground meteorological data. This framework includes extraction of the wine grape planting area and Figure 13. Late frost damaged area of wine grapes in the study area in April 2020.

Discussion
In this paper, a practical remote sensing monitoring framework of wine grape late frost damage was proposed using multi-source satellite data and auxiliary ground meteorological data. This framework includes extraction of the wine grape planting area and estimating the daily minimum air temperature using satellite data.
In April 2020, wine grapes in the study area suffered severe late frost damage. According to the estimation results of T min in this paper, 14 frost days occurred in April, which may have caused the leaves and buds of wine grapes to be frozen. During 20-25 April, several extreme low air temperature events occurred, particularly the sudden drop in air temperature in the early morning of 24 April, which turned the wine chateaus with slight or moderate late frost damage in the early stage into moderate or severe late frost damage finally. In addition, because most of the wine grapes had begun to sprout, this growth stage was greatly affected by low air temperature; the whole region suffered severe late frost damage with long duration, low air temperature, and wide affected area, especially in low-lying areas. Because cold air is heavier than warm air, and there are more water molecules in cold air. When the cold air sinks and collects in valleys or low-lying areas, frost forms. Through field investigations and meteorological monitoring, it is found that this late frost was the mixed frost of radiation and advection [70], causing different degrees of freezing damage to the buds, leaves, and tender shoots of wine grapes in the germination and leaf-expansion stages. At the same time, because the occurrence of late frost was mostly around the germination and leaf-expansion stages of wine grapes, the sensitivity and tolerance of different wine grape varieties to low air temperature are different.
This late frost had a great impact on the mid-late maturing red wine grape varieties, such as Mathelan, Merlot, and Shiraz, and white wine grape varieties, such as Riesling and Chardonnay. The effect on Cabernet Sauvignon, a late maturing variety, was relatively slight. The influence of cold air can be reduced by smoking, fan, drip irrigation, and shallow tillage to prevent late frost damage. In view of the late frost damages that happened, we can help promote bud germination and inflorescence growth to reduce the yield loss and increase the fruit setting rate by providing a small amount of irrigation, increasing the nitrogen fertilizer, deep plowing the soil, and covering grapevines with plastic film.
In previous studies, the agricultural research on wine grapes mostly focused on agronomic mechanisms and cultivation, while research on late frost damage is mainly based on field experiments and meteorological stations [71][72][73], which have certain spatiotemporal limitations. It is rare to use remote sensing to monitor late frost damage in wine grapes and to estimate the damaged area at a large scale. In view of the above problems, this paper proposed a monitoring model of wine grape late frost damage based on remote sensing, which made full use of multi-source data, such as long-term series of remote sensing satellite data, ground meteorological station data, terrain data, and phenological data of wine grape, to monitor late frost damage in wine grapes and estimate the damaged area at a large scale. Previous researches on wine grapes focused on late frost duration patterns [4], the impact of climate change [74][75][76], risk assessment [77,78], and regional planning [2]. Compared with the daily minimum temperature estimation and frost monitoring models used in these studies, the multi-source data used in our model ensured the integrity of the model framework and the high accuracy of the results. In addition, the advantage of our model was that we used fewer modeling indices, but obtained higher accuracy. Although the complete research framework could apply to wine grape remote sensing monitoring, more work is needed to explore the optimal scale for wine grape monitoring. If the LST is downscaled to a higher spatial resolution, more potentially available information may be added. In addition, we should further implement a corresponding evaluation index of late frost damage for different kinds of wine grapes to comprehensively evaluate the late frost damage under cumulative effects, such as temperature, precipitation, and soil condition. Meanwhile, the application of the remote sensing monitoring framework for other crops in other regions will be explored from the perspective of practical use.

Conclusions
This study presented a practical remote sensing monitoring framework for late frost damage in wine grapes based on in-situ measurements and multi-source satellite data. The in situ data included meteorological data and field survey data. These data obtained the accurate wine grape planting area, daily T min data, and late frost damage in wine grapes, but were limited by their uneven and sparse distribution. In this study, we extracted the wine grape planting area by selecting the optimal feature subset, estimating the daily minimum air temperature (T min ) with the spatial resolution of 100 m, and mapped the wine grape late frost damage in April 2020. Spatially, about 41.12% of the vineyards suffered severe frost damage, and the total affected area was about 16,381 ha.
The results of late frost damage obtained by estimating the T min were in good agreement with the statistics of the agricultural meteorological disaster department. The proposed framework innovatively solved the problems of a lack of meteorological data or coarse spatial resolution, and realized accurate monitoring of full coverage in the study area. It can also provide a reference and technical support for estimations of the wine grape area, large scale dynamic monitoring of late frost damage, and agro-meteorological services in the whole planting area in the future.