A NOVEL PROCEDURE FOR GENERATION OF SAR-DERIVED ZTD MAPS FOR WEATHER PREDICTION: APPLICATION TO SOUTH AFRICA USE CASE
A NOVEL PROCEDURE FOR GENERATION OF SAR-DERIVED ZTD MAPS FOR WEATHER PREDICTION: APPLICATION TO SOUTH AFRICA USE XXXX
X. X. Xxxxxxxx 0 *, X. Xxxxxxx 1, N. Petrushevsky1, A.M. Xxxxxxxxx 1, X. Xxxxxx 2, A. N. Meroni 2, X.Xxxxxxxxxx 2, X. Xxxxxx 3 1 Politecnico di Milano - Dipartimento di Elettronica, Informazione e Bioingegneria (DEIB), Milan, Italy (moniaelisa.xxxxxxxx,
marco.manzoni, xxxxx.xxxxxxxxxxxx, xxxxxx.xxxxxxxxxxxxxx)@xxxxxx.xx
2 Politecnico di Milano - Dipartimento di Ingegneria Civile e Ambientale (DICA), Milan, Italy (xxxxxxxx.xxxxxx, agostinoniyonkuru.xxxxxx, xxxxxxxxxx.xxxxxxxxxx)@xxxxxx.xx
3 Centro Internazionale in Monitoraggio Ambientale Research Foundation, Savona, Italy - xxxxxxx.xxxxxx@xxxxxxxxxxxxxx.xxx
KEY WORDS: SAR, GNSS, ZTD, water vapor, weather prediction, TWIGA
ABSTRACT:
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.
.
1. 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 (Xxx et al., 1993; Xxxxxxxx et al., 2004; Xxxxxx 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 (Xxxxxxxx, 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 (Xxxx et al., 2016).
Experiments on the ingestion in NWPMs of water vapor data derived by ENVISAT-ASAR satellite are reported by Xxxxxxxx et al. (2015) and Xxxxxx 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
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 (Xxxxxx et al., 2018; Xxxxxxx et al., 2019; Xxxxxxx et al., 2019; Xxxxxx 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.
2. 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.
SLC Sentinel-1
stack
Precise DEM orbit file
Ionospheric
maps
Ionospheric map estimation and
compensation
Differential ZTD estimation Unwrapping
GNSS time
series
Precise orbit correction
Phase offset compensation Geocoding
Mosaicking
GACOS
Absolute ZTD
maps (GeoTIFF)
Mask
Nodata filling
Gaussian smoothing
worldfile
Absolute ZTD
maps (JPEG)
Auxiliary
file
Mask
Absolute ZTD
maps (GeoTIFF)
Masking
Rescaling
JPEG to GeoTIFF conversion
GeoTIFF to JPEG conversion
Scaling
Master reconstruction
Coregistration Debursting
Topographic phase compensation
ZTD MAP ESTIMATION
PRE-PROCESSING
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).
COMPRESSION/DECOMPRESSION
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.
2.1 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, Xxxxx Xxxxxxxxx, 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 (Xxxxx et al., 2016). The estimation is performed jointly over the whole stack of images, reaching in this way a better accuracy.
2.2 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 (Xxxxxxxxx, Xxxxxxxxx, 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 kth target, the following cost function:
ποΏ½π(ππ ) = arg min π π(βπ βπ βπ πΌ(ππβπ,πβπ) β πππ(ποΏ½πβποΏ½π)), | (1) |
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
pseudo-range measured at the time of the SAR acquisitions (Xxxxxxx 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, Xxxxxxxx, 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).
2.3 Compression/Decompression
The final aim of the procedure is to generate absolute ZTD products to be supplied to African stakeholders through cloud- based 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 pre- processing 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 back- conversion 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.
3. RESULTS AND DISCUSSION
3.1 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 |
Relative Orbit | 43 |
Frames | 1082, 1087, 1092, 1097, 1102 |
Pass | Ascending |
Number of images | 6 |
First image acquisition date | 2018-02-01 |
Last image acquisition date | 2018-04-02 |
Revisiting time | 12 days |
Table 1. Information related to Synthetic Aperture Radar (SAR) imagery collected for South Africa use case.
.
As highlighted in Figure 2, an acceptable number of GNSS South African Network stations is available within the area of interest outlined by SAR frames.
SAR frames GNSS
Figure 2. South Africa use case: SAR frames and GNSS stations.
3.2 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 km2.
Figure 3. Correlation between SAR and differential ZTD reconstructed in 17 GNSS stations before and after residual orbit compensation.
3.3 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.
Absolute ZTD [cm]
264
160
Figure 4. InSAR absolute ZTD map (09 March 2018 β 21 march 2018) obtained by merging five frames.
Min: -23.6
Max: 20.9
Mean: -0.00009
Std Dev: 1.1
3.0
Number of pixels (x 106)
2.5
2.0
1.5
1.0
0.5
0
-25 -20 -15 -10 -5 0 5 10 15 20
Pixel values (mm)
Figure 5. Histogram of the differences between the original GeoTIFF and the smoothed one obtained during compression.
Min: -24.0
Max: 22.4
Mean: -0.006
Std Dev: 1.3
2.5
Number of pixels (x 106)
2.0
1.5
1.0
0.5
0
-25 -20 -15 -10 -5 0 5 10 15 20
Pixel values (mm)
Figure 6. Histogram of the differences between the smoothed GeoTIFF and the recovered one.
4. 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).
Funding: TWIGA has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No 776691
REFERENCES
Xx Xxx, F., Xxxxx Xxxxxxxxx, A., 2006. TOPSAR: Terrain Observation by Progressive Scans. IEEE Transactions on Geoscience and Remote Sensing, 44 (9), 2352-2360.
Xxxxxxxx, A., Monti-Xxxxxxxxx, A., Xxxxx, C., Xxxxx, F., Xxxxxxxxx, D., 2007. InSAR Principles: Guidelines for SAR Interferometry Processing and Interpretation. TM-19. ESA Publications, The Netherlands.
Xxxxxxx, D.C., Xxxxx, M.D., 1998: Two-dimensional phase unwrapping: Theory, algorithms, and software. New York: Wiley.
Xxxxxxxxx, X., Xxxxxxxxx, S., 2008. On the Exploitation of Target Statistics for SAR Interferometry Applications. IEEE Trans. Geosci. Remote Sensing, 46, 3436β3443.
Xxxxxxx, X., Xxxxx, X., Xxxxx, X., Xxxx, G., xx Xxxx, S., Xxxxxxxx, E., Xxxx, X., Xxxxxxx, R., Xxxxxxx, G., Vedel, X., Xxxxxx, M., 2016. Review of the state of the art and future prospects of the ground-based GNSS meteorology in Europe. Atmos. Meas. Tech., 9, 5385β5406.
Xxx, X.-X., Xxx, X.-X., Xxxxxxxxx, X.X., 0000. Assimilation of precipitable water measurements into a mesoscale numerical model. Mon. Wea. Rev., 121, 1215β1238.
Xxxxxxx, M., Xxxxxx, X., Pulvirenti, L., Xxxxxx, A.N., Xxxx, G., Xxxxxxxxx, N., Xxxxxxx, F.S., Luini, L., Xxxxxx, G., Xxxxxxx, E., Xxxxx, A, Xxxxxxxxxxx, G., Xxxxxxxxxx, S., Xxxxx Xxxxxxxxx,X., Xxxx, K, Terzo, O, Xxxxx, A., Xxxxxxx, X, Xxxxxxxxxxxxx, X., Xxxxxx, X., 0000. A Synergistic Use of a High-Resolution Numerical Weather Prediction Model and High-Resolution Earth Observation Products to Improve Precipitation Forecast. Remote Sens., 11 (20), 2387.
Xxxxxxx, M., Monti-Xxxxxxxxx, A.V., Xxxxxxx, E., Xxxxxx, G., 2020. Joint Exploitation of SAR and GNSS for Atmospheric Phase Screens Retrieval Aimed at Numerical Weather Prediction Model Ingestion. Remote Sens., 12, 654.
Xxxxxx, S., Xxx, X., Xxxx, T., Xxxxxxxxx D., Xxxxx, X., 2007. Influence of GPS Precipitable Water Vapor Retrievals on Quantitative Precipitation Forecasting in Southern California. Journal of Applied Meteorology and Climatology, 46(11), 1828- 1839.
Xxxxxx X., Xxxx, X., Xxxx, G., Member S., Xxxxxxx, J., 2016. Three-Dimensional Variational Assimilation of InSAR PWV Using the WRFDA Model. IEEE Trans. Geosci. Remote Sens., 54, 7323β7330.
Xxxxxx, P., Xxxxxxx, P. M. A., Xxxx, G., CatalΓ£o, X., Xxxxx, P., Xxxx, R., 2018. Assimilating InSAR maps of water vapor to improve heavy rainfall forecasts: A case study with two successive storms. Journal of Geophysical Research: Atmospheres, 123(7), 3341-3355.
Xxxxxx, P., Xxxxxxx, P. M. A., Xxxx, G., CatalΓ£o, J., 2020. Continuous Multitrack Assimilation of Sentinelβ1 Precipitable Water Vapor Maps for Numerical Weather Prediction: How Far Can We Go With Current InSAR Data?. Journal of Geophysical Research: Atmospheres, 126(3), e2020JD034171.
Xxxxxx, A.N., Xxxxxxxxx, M., Xxxxxx, X., et al., 2020. On the Definition of the Strategy to Obtain Absolute InSAR Zenith Total Delay Maps for Meteorological Applications. Frontiers in Earth Science, 8, 359.
Xxxxxxx, P.M.A., Xxxxxx, P., Xxxx, G., CatalΓ£o, J., Xxxx, R., Xxxxxxxx, M., 2019. InSAR meteorology: Highβresolution geodetic data can increase atmospheric predictability. Geophysical Research Letters, 46(5), 2949β2955.
Xxxxxxxx, X., Xxxxxxx, X., Xxxxxxx, X., 0000. Data assimilation of GPS precipitable water vapor into the JMA mesoscale numerical weather prediction model and its impact on rainfall forecasts. J. Meteorol. Soc. Jpn. Ser., 82(1B), 441β452.
Xxxxxxxx, E., Xxxxxxxx, R., Xxxxxx, D., Xxxxxxxxxx, G., Xxxxxxxx, D., Xxxxxxxxx, N., et al. 2015. InSAR water vapor data assimilation into mesoscale model MM5: Technique and pilot study. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 8(8), 3859β 3875.
Sansosti, X., Xxxxxxxxx, X., Xxxxxxx, X., Xxxxxxxx, X., Xxxxxxx, X., 0000. Geometrical SAR image registration. IEEE Trans. Geosci. Remote Sens., 44, 2861β2870.
Xxxxx, D., Xxxxxxxx, A., 2019. Guide to Sentinel-1 Geocoding. University of Zurich, Report UZH-S1-GC-AD.
Xxxxxxxx, X., 2007. Parameterization schemes: Keys to understanding numerical weather prediction models. Cambridge, UK: Cambridge University Press.
Xxxx, X., Xxxx, M., Xxxxx, X., Xx, X., Xx, W., 2016. High- spatial-resolution mapping of precipitable water vapour using SAR interferograms, GPS observations and ERA-Interim reanalysis. Atmos. Meas. Tech., 9, 4487β4501.