Gridding of L2a and L2b EGMS data to minimise the overlap artefacts

By Alexis Hrysiewicz, Jan. 2025

For your information

This documentation is an early version.

This example offers the way in order to download EGMS data in L2a (or L2b) level and to use interpolations for minimising of burst-overlap artefacts. This example therefore proposes an advanced use of EGMS toolkit.

First, we import the EGMS toolkit package and the required packages.

from EGMStoolkit.classes import EGMSS1burstIDapi
from EGMStoolkit.classes import EGMSS1ROIapi
from EGMStoolkit.classes import EGMSdownloaderapi
from EGMStoolkit.functions import egmsdatatools
import fiona, os
from shapely.geometry import Polygon, mapping, shape, LineString, Point
import pyproj
import numpy as np
import pandas as pd
from osgeo import gdal, osr
import glob

Then we detect the available data based on the following parameters: * L2a level; * Region of Interest over Dublin, Ireland * 2018-2020 release * Track number 1 * Ascending data

info = EGMSS1burstIDapi.S1burstIDmap(verbose=True,log=None)
info.downloadfile(verbose=True)

ROIpara = EGMSS1ROIapi.S1ROIparameter(verbose=False)

ROIpara.egmslevel = 'L2a'
ROIpara.bbox = [-6.427059639290446,53.2606655698541,-6.0952332730202095,53.41811986118854]
ROIpara.release = '2018_2022'

ROIpara.createROI(verbose=True)
ROIpara.detectfromIDmap(info,Track_user=1,Pass_user='Ascending',verbose=True)

ROIpara.displaymap(output='fig_search_adexample1.jpg')

The search parameters are now ready. We can send the request to retreive the files.

downloadpara = EGMSdownloaderapi.egmsdownloader(verbose=False)
downloadpara.updatelist(infoS1ROIparameter=ROIpara)
downloadpara.printlist(verbose=True)

So there are 4 files to be downloaded.

downloadpara.token = 'xxxxxxxxxxxxxxxxxxxx'
downloadpara.download(outputdir='./Output')
downloadpara.unzipfile(outputdir='./Output')

The files therefore are downloaded and unzipped. For this example, only two parameters will be manipulated: (1) mean_velocity, (2) mean_velocity_std.

Note

If all parameters need to be stored, a parameter variable can be generated from any file header of EGMS dataset. For example,

with open('Output/L2a/2018_2022/EGMS_L2a_001_0314_IW2_VV_2018_2022_1/EGMS_L2a_001_0314_IW2_VV_2018_2022_1.csv','r') as fi:
    first_line = fi.readline()
paralist = first_line.split()[0].split(',')

Here, all parameters (i.e., mean_velocity) and displacement for each date will be manipulated and stored. There are 294 parameters, including 269 observation dates.

The following lines will run the merging. The duplicated-point removal must be deactived as the interpolation will be average the observation points.

paralist = ['mean_velocity','mean_velocity_std']

egmsdatatools.datamergingcsv(infoEGMSdownloader=downloadpara,inputdir='./Output',outputdir='./Output',mode='onlist',verbose=True,paratosave=paralist,__removeduplicate__=False)

Now, we therefore have a .csv file with the EGMS observations. The idea is to grid these observations into the regular grid (no clip or cropped will be required). Based on our Region of Interest shapefile, we can extract the final image extent.

schema = {
    'geometry': 'Polygon',
    'properties' : {'id':'int'}
    }

latlon_to_meter = pyproj.Transformer.from_crs('epsg:4326', 'epsg:3035')

with fiona.open('bbox.shp','r','ESRI Shapefile', schema) as shpfile:
    coordinates = []
    Xcoord = []
    Ycoord = []
    for feature in shpfile:
        line = shape(feature['geometry'])
        if isinstance(line, LineString):
            for index, point in enumerate(line.coords):
                if index == 0:
                    first_pt = point
                a, b = latlon_to_meter.transform(point[1],point[0])
                Xcoord.append(b)
                Ycoord.append(a)

And we can create the interpolation parameters. The grid is aligned to the coordinates used by EGMS for L3 datasets.

paragrid = dict()
paragrid['Xmin'] = np.fix(np.min(Xcoord)/100)*100
paragrid['Ymin'] = np.fix(np.min(Ycoord)/100)*100
paragrid['Xmax'] = np.fix((np.max(Xcoord)/100)+1)*100
paragrid['Ymax'] = np.fix((np.max(Ycoord)/100)+1)*100
paragrid['xres'] = 100
paragrid['yres'] = 100
paragrid['algo'] = 'average:radius1=100:radius2=100:angle=0.0:nodata=-9999' # Alfgorithm used and options
paragrid['variable'] = ','.join(paralist)

Finally, we can run the interpolation via the EGMS toolkit tool.

egmsdatatools.datagridding(inputdir='./Output',outputdir='./Geotiff',namefile='all',verbose=True,paragrid=paragrid,AREA_OR_POINT='Area')

---

This following sections are optional but improve the result quality.

Even if the rasters are correct, the GDAL-based interpolation uses a circular kernel. It means that we extrapolated some points. However the goal is to avoid any extrapolation, and to ideally interpolate only pixels with EGMS measurement points. The solution is to create a mask by counting numbers of measurement points.

data = pd.read_csv('Output/EGMS_L2a_001_VV_2018_2022_1.csv',delimiter=';')

centersX=np.arange(paragrid['Xmin'],paragrid['Xmax']+paragrid['xres'],paragrid['xres']) + paragrid['xres']/2
centersY=np.arange(paragrid['Ymin'],paragrid['Ymax']+paragrid['yres'],paragrid['yres']) - paragrid['yres']/2

H, *_ = np.histogram2d(data['easting'],data['northing'],bins=[centersX, centersY])
H = np.flipud(H.T)

The variable H contains the numbers of measurement points for each pixel.

To be more conservative, we can define a threshold to keep pixels containing at least x points. Here only the pixels with at least 1 measurement points will be kept.

H[H<1] = 0

Then we save a GTiff image of this mask using GDAL:

image_size = H.shape
nx = image_size[0]
ny = image_size[1]

dst_ds = gdal.GetDriverByName('GTiff').Create('Geotiff/mask.tif', ny, nx, 1)
geotransform = (paragrid['Xmin'], paragrid['xres'], 0, paragrid['Ymax'], 0, -paragrid['yres'])

dst_ds.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(3035)
dst_ds.SetProjection(srs.ExportToWkt())
dst_ds.GetRasterBand(1).WriteArray(H)
dst_ds.GetRasterBand(1).SetNoDataValue(np.nan)
dst_ds.FlushCache()
dst_ds = None

The last code section is for the application of the computed mask.

# Read the mask
src_ds = gdal.Open('Geotiff/mask.tif')
mask = src_ds.GetRasterBand(1).ReadAsArray()
src_ds = None

# Loop for each file
for fi in glob.glob('Geotiff/EGMS_*.tif'):

    src_ds = gdal.Open(fi)
    band = src_ds.GetRasterBand(1).ReadAsArray()
    band[mask==0] = -9999

    dst_ds = gdal.GetDriverByName("GTiff").Create(fi.replace('.tif','_masked.tif'), src_ds.RasterXSize, src_ds.RasterYSize, 1, gdal.GDT_Float64)
    dst_ds.SetGeoTransform(src_ds.GetGeoTransform())
    dst_ds.SetProjection(src_ds.GetProjection())
    dst_ds.GetRasterBand(1).WriteArray(band)
    dst_ds.GetRasterBand(1).SetNoDataValue(-9999)
    dst_ds.FlushCache()

    src_ds = None
    dst_ds = None

Each file ending by _masked.tif - stored in the Geotiff directory - will be the EGMS rasterised datasets.