A NOVEL PROCEDURE FOR GENERATION OF SAR-DERIVED ZTD MAPS FOR WEATHER PREDICTION: APPLICATION TO SOUTH AFRICA USE CASE

The knowledge of tropospheric water vapor distribution can significantly improve the accuracy of Numerical Weather Prediction (NWP) models. The present work proposes an automatic and fast procedure for generating reliable water vapor products from the synergic use of Sentinel-1 Synthetic Aperture Radar (SAR) imagery and Global Navigation Satellite System (GNSS) observations. Moreover, a compression method able to drastically reduce, without significant accuracy loss, the water vapor dataset dimension has been implemented to facilitate the sharing through cloud services. The activities have been carried in the EU H2020 TWIGA project framework, aimed at providing water vapor maps at Technology Readiness Level 7. .


INTRODUCTION
Numerical Weather Prediction models (NWPMs) currently represent the most important tool for weather forecasting and provide crucial information for lives and property protection and conscious human activities managing. Enhancing the performance of these models, especially for the prediction of heavy convective precipitation events, is still a challenging task.
It has been demonstrated in the literature that the knowledge of tropospheric water vapor distribution can significantly improve the accuracy of NWPMs forecasts. The assimilation in NWPMs of water vapor observations derived by meteorological and Global Navigation Satellite System (GNSS) stations represents a well-established and effective technique (Guerova et al., 2016), and several studies showed its positive impact on the prediction of precipitation events (Kuo et al., 1993;Nakamura et al., 2004;Marcus et al., 2007). Although characterized by high temporal resolution, the GNSS observations are usually provided as sparse datasets due to the low density of available measurements. This means that they cannot comprehensively describe the spatial variations of water vapor, which in storm phenomena can occur within tens of meters (Stensrud, 2007).
Satellite observations can undoubtedly give a valuable and complementary contribution in this sense, as they can provide additional water vapor observations over broad areas. Of particular interest is the Synthetic Aperture Radar (SAR), which can exploit the SAR Interferometry (InSAR) technique to obtain water vapor distribution products characterized by high spatial resolution and high accuracy, about 2 mm (Tang et al., 2016). Experiments on the ingestion in NWPMs of water vapor data derived by ENVISAT-ASAR satellite are reported by Pichelli et al. (2015) and Mateus et al. (2016). The authors observed an improvement in forecasting weak to moderate precipitation. In recent years, researches are increasing thanks to the availability of imagery from Sentinel-1 satellites of the European * Corresponding author Space Agency (ESA). In orbit since 2014-2016, the satellites have the advantage to guarantee the free availability of time series of SAR data with an unprecedented time coverage, i.e. up to 1-3 days. The temporal high-resolution of water vapor variability information is crucial in obtaining benefits from the assimilation process. Several investigations have been performed taking advantage of Sentinel-1 data (Mateus et al., 2018;Miranda et al., 2019;Lagasio et al., 2019;Mateus et al., 2020), confirming the InSAR water vapor products potential in enhancing NWP models accuracy.
The present work fits within this context and has been carried out in the framework of TWIGA, the EU Horizon2020 project aimed at providing currently unavailable geo-information on weather, water, and climate for sub-Saharan Africa. Here, an automatic and fast procedure integrating Sentinel-1 SAR imagery data with in situ sensors observations has been designed and tuned to obtain highly accurate, dense, and wide water vapor products, also called Zenith Total Delay (ZTD). These products will be supplied to local meteorological services to be ingested by NWP models and improve the prediction of heavy rainfall.
The paper is structured as follows. In Section 2, the adopted methodology for water vapor product generation is illustrated. Section 3 provides information about the use case area, the processed SAR imagery, and the obtained results. In the last section, conclusions are reported.

METHOD
The main steps for generating absolute ZTD maps from a stack of Single Look Complex (SLC) Sentinel-1 SAR images are outlined in Figure 1. The procedure combines into a unique processing chain both standard techniques and some novel contributions to enhance the quality of the estimated water vapor maps.
The proposed methodology exploits the availability of a stack of complex SAR images acquired over the same area at different epochs to estimate the difference in the optical path between radar and ground targets (Ferretti et al., 2007).

Figure 1.
Implemented procedure to generate ZTD maps from a stack of SAR data.
The computed optical path variation can be due to a few contributions: i) the possible displacement of the target; ii) the different satellite view angle between the acquisitions, which generates local topographic effects; and iii) the changes in atmospheric condition (mainly the troposphere and the ionosphere).
As we are interested in estimating the tropospheric effect only, the other contributions, if not negligible, must be removed. This step is performed in the pre-processing phase, which prepares the SAR image stack for the interferometric process.

Pre-processing
The pre-processing is mainly performed by exploiting existing tools in the free and open-source software SNAP (Sentinel Application Platform). It is based on well-known techniques such as co-registration (Sansosti et al., 2006) and debursting (De Zan, Monti Guarnieri, 2006). The topographic contribution to the interferogram phases is compensated by using precise orbit information and a DEM of the scene. The ionospheric contribution is estimated and then compensated by a novel method based on the split spectrum technique (Gomba et al., 2016). The estimation is performed jointly over the whole stack of images, reaching in this way a better accuracy.

ZTD map estimation
The processing phase includes both the interferometric technique and all the steps required to obtain absolute ZTD maps. The estimation of differential ZTD maps is usually performed using techniques that exploit the Permanent Scatterers (PS). The drawback of these approaches is that PS density is usually low in rural areas. Thus, the final water vapor products cannot be as dense as required. Our procedure implements the phase linking method (Guarnieri, Tebaldini, 2008) to extract differential ZTD dense maps by exploiting both PS and Distributed Scatterers (DS).
Given a stack of N coregistered SAR images, the phase linking technique generates all the possible N(N-1)/2 interferograms. The procedure iterates by minimizing, in each small window centered on the k th target, the following cost function: where � ( ) = estimated phase I = interferogram = k-th target Notice that the cost function involves only phase differences. Therefore, a constant can be added to all the terms � ( ) without changing the solution. In other words, for each pixel Pk, all the phases can be estimated but for a constant phase offset.
After applying the unwrapping technique (Ghiglia, Pritt, 1998), the estimated differential ZTD products need to be further enhanced. Due to the lack of millimetric accuracy in the satellite orbit determination, the pre-processing phase cannot completely remove the topographic contribution. The joint exploitation of SAR data and GNSS atmospheric products allows for the correction of this error. It is worth noting that the synergic use of satellites (SAR) and in-situ (GNSS) observations is one of the significant innovations in TWIGA paradigm.
The orbital parameters are estimated by matching the SAR phases evaluated in correspondence of the GNSS and the GNSS The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition) pseudo-range measured at the time of the SAR acquisitions (Manzoni et al., 2020). The matching is performed in a robust L1 norm to provide the needed orbital error parameters.
The subsequent step of the procedure still involves GNSS measures. Since InSAR is a differential technique both in space and in time, a constant phase term is missing in all the images. Thus, GNSS observations are used to estimate the missing constants required to adjust the interferometric phase. Once calibrated, the differential ZTD maps are then geocoded (Small, Schubert, 2019) and, if adjacent ZTD images exist, they are merged into a unique wide map using a mosaicking tool.
Finally, the differential ZTD maps are computed with respect to an initial map, the so-called master. To obtain absolute ZTD products the master map has to be estimated. For this purpose, the map provided by the Generic Atmospheric Correction Online Service (GACOS) is exploited following Meroni et al. (2020).

Compression/Decompression
The final aim of the procedure is to generate absolute ZTD products to be supplied to African stakeholders through cloudbased services, such as File Transfer Protocol (FTP). The main drawback is certainly related to the size of the GeoTIFF ZTD maps, which can reach the order of hundreds of MegaBytes per map. This problem has been overcome by implementing a compression method that converts GeoTIFFs to georeferenced JPEGs.
Since JPEG is a lossy image compression technique, some issues arise in correspondence of image discontinuities, e.g. at the edges between data and no data or when a sharp difference in adjacent pixels values exist. Hence, before the conversion, some preprocessing is carried out to face these artifacts. Firstly, a data interpolation is performed for filling regions with no data values. Then, a multidimensional Gaussian filter is applied to smooth sharp differences in pixel values. After that, the GeoTIFF pixel values are rescaled in the range 0-255, supported by JPEG format, and the image is converted to JPEG format. A worldfile containing information for georeferencing is generated. The information for backconversion and rescaling is provided in an XML auxiliary file. Moreover, a GIF mask file contains the no data values original distribution. The JPEG compressed image can be easily stored, transmitted to the cloud, and delivered to users.

Use case
The proposed procedure has been exploited to generate absolute ZTD maps for the whole South Africa country. Figure 2 shows the footprints of the Sentinel-1A frames (blue boxes) which ensure the use case area coverage: five frames belonging to ascending relative orbit 43 were selected. Table 1 provides the main characteristics of the SAR dataset used for the study.
The temporal period of images has been chosen based on a severe storm event that hit the South-Africa, particularly the Johannesburg and Pretoria cities, between 22 and 24 March 2018. For each Sentinel-1 frame, a stack of 6 images at 12 days revisit frequency related to timespan February-April 2018 has been considered. Of particular interest is the imagery of 21 March 2018, which can provide useful information for NWP models since it has been acquired just a few hours before the storm event.
It is worth noting that the short temporal extent of the stack is required to avoid the phase noise induced by temporal decorrelation of the distributed target as well as the impact of possible subsidences.

Dataset characteristics South Africa use case
Sensor S1A  Figure 2, an acceptable number of GNSS South African Network stations is available within the area of interest outlined by SAR frames.  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition)

Differential and absolute ZTD maps
The SAR-derived differential ZTD maps have been computed by considering an estimation window of 400 m x 400 m, which is enough to cope with decorrelation noise. The precise orbit correction and the phase offset compensation take advantage of the South African Network GNSS observations, which have been processed to have ionosphere-free ZTD measurements. Figure 3 proposes a comparison between SAR-derived and GNSS differential ZTD values computed for the interval time 25-02-2018 and 09-03-2018. Results prove the procedure's capability to obtain reliable ZTD measurements. In fact, a correlation value of 0.984 is observed even before applying precise calibration. Once this correction is performed, the correlation index increases to 0.993. Figure 4 shows an example of absolute ZTD map obtained after geocoding, mosaicking, and master reconstruction steps. The long strip covers an area of about 850x250 km 2 .

Figure 3.
Correlation between SAR and differential ZTD reconstructed in 17 GNSS stations before and after residual orbit compensation.

Compression/Decompression method: first experiments
A single ZTD map of South Africa use case in Packbits GeoTIFF format is 170 MB in size. After the compression method, the JPEG image is about 0.5 MB.
The decompression method recovers the original GeoTIFF, with some errors due to the Gaussian smoothing and the JPEG compression. Figure 5 shows the histogram of differences between the original and the smoothed GeoTIFF, with pixel values converted in millimeters. The mean error is zero, while the standard deviation is about one millimeter, which is below the ZTD accuracy discussed in the introduction. There are peaks in changes of about 20 mm that can be attributed to the smoothing effect in correspondence of discontinuities. Figure 6 reports the histogram of differences between the smoothed GeoTIFF and the recovered one. It is worth noting that the final map can be affected by additional minor errors introduced by the rescaling operation during compression. The rescaling implies to round float value to an integer, and this operation cannot be reversed.  The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition) Figure 6. Histogram of the differences between the smoothed GeoTIFF and the recovered one.

CONCLUSIONS
The presented work proposes a novel procedure developed for the generation of ZTD products from the synergic use of SAR and GNSS data. The procedure has been implemented in the framework of the TWIGA project with the aim to provide to African stakeholders a fast, automatic, and easy-to-replicate method to retrieve useful geo-information, currently unavailable, for weather forecasting and better management of water resources. ZTD products have been derived for the use case of South Africa by taking advantage of the ubiquitous and freely available time series of Sentinel-1 SAR data. A compression method to drastically reduce ZTD file size has been developed to allow data delivery and storage. Further experiments need to be performed to evaluate the effects of compression on product quality. From the computational side, a ZTD product covering an area of 850 km x 250 km was generated in a reasonably fast time, i.e. about half of an hour. The STAMPS software, which works only with PS, obtains water-vapor maps in several hours.
As a final remark, the procedure has been implemented by making the maximum use of open-free tools, mainly the ESA SNAP tool and several Python libraries (numpy, scipy, rasterio, Python-gdal). The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2021 XXIV ISPRS Congress (2021 edition)