Output resolution

[sixy6e]

The workflow is outputting ~7m metre pixels. The code is automatically determining an output resolution based on the input resolution.

https://github.com/GeoscienceAustralia/gamma_insar/blob/pygamma_workflow/insar/coregister_dem.py#L372

An enhancement to this could be to default the oversampling factor to 1, and still allow users to define their own oversampling factor.

[sixy6e]

What are your thoughts at @mcgarth @tfuhrmann ?

[mcgarth]

For Sentinel-1 ARD products (be it backscatter of InSAR) I don't think there is any sense in over-sampling beyond the resolution of the input DEM (i.e. SRTM 1-arcsecond). It's just increasing the disk size of the products unnecessarily in my view.

[tfuhrmann]

The DEM used in the GAMMA workflow along with its oversampling factor define the spatial resolution of the final ARD products. I think this resolution should be chosen according to the resolution of the input SAR data (and the chosen multi-looking factor for the SAR data). For example, if we process Sentinel-1 data at a spatial resolution of 30 m x 30 m, the SRTM 1-arcsecond would match without oversampling (since 1 arcsecond is approx 30 m x 30 m). If we decide to use full resolution Sentinel-1 data (square pixels) at a 15 m x 15 m resolution, the SRTM 1-arcsecond would need to be oversampled by a factor of 2.

[sixy6e]

Probably something that can be parameterised with a default of 1.

[tfuhrmann]

As a temporary fix (and to save disk space), I've hard-coded the DEM oversampling factor in the "develop" branch to 1. I'm currently testing the GAMMA workflow with that setting, and can confirm that the final output data sets then have a spatial resolution of around 30 m x 30 m.

For the bash workflow, I would like to propose to add a couple of lines of code to "calc_multi-look_values.bash" to read the spatial resolution of the input DEM (given in /gamma_dem/*.dem.par) and calculate an oversampling factor that fits to the multi-looked spatial resolution of the radar data. The ground pixel resolution of the radar is already calculated in the .bash script (parameters: "azspmean" and "grrgspmean". Any thoughts?

[sixy6e]

A function of some sort would be useful, but anything developed for your BASH workflow would need to be transcribed into a Python alternative. At least doing anything numerical would be easier in Python than BASH.

[mcgarth]

At least doing anything numerical would be easier in Python than BASH.

Damn straight!

[sixy6e]

@tfuhrmann can you provide some formula details on how you would calculate a nominal oversampling factor that fits the multi-looked spatial resolution? My assumption is that an oversampling factor can be a float (not just multiples of 2).

In the short term I've set a default of 1, but I've assumed integers based on the original code only specifying 2, 4 or 8.

[tfuhrmann]

I've looked up the info on the GAMMA function "gc_map" which is used for the oversampling (see below). This suggests the oversampling factor (lat_ovr and lon_ovr) can be float.

I can generate some code to calculate an oversampling factor that suits the spatial resolution of the multi-looked SAR data. Are we all happy with using this approach?

` gc_map Calculate terrain-geocoding lookup table and DEM derived data products Copyright 2019, Gamma Remote Sensing, v1.1 7-Nov-2019 cm/clw/uw

usage: gc_map [lat_ovr] [lon_ovr] [sim_sar] [u] [v] [inc] [psi] [pix] [ls_map] [frame] [ls_mode] [r_ovr]

input parameters: MLI_par (input) ISP MLI or SLC image parameter file (slant range geometry) OFF_par (input) ISP offset/interferogram parameter file (enter - if geocoding SLC or MLI data) DEM_par (input) DEM/MAP parameter file DEM (input) DEM data file (or constant height value) DEM_seg_par (input/output) DEM/MAP segment parameter file used for output products

NOTE: If already exists, then the output DEM parameters will be read from this file, otherwise they are estimated from the image data. If and are the same file, is ignored.

DEM_seg (output) DEM segment used for output products, interpolated if lat_ovr > 1.0 or lon_ovr > 1.0 lookup_table (output) geocoding lookup table (fcomplex) lat_ovr latitude or northing output DEM oversampling factor (enter - for default: 1.0) lon_ovr longitude or easting output DEM oversampling factor (enter - for default: 1.0) sim_sar (output) simulated SAR backscatter image in DEM geometry (enter - for none) u (output) zenith angle of surface normal vector n (angle between z and n, enter - for none) v (output) orientation angle of n (between x and projection of n in xy plane, enter - for none) inc (output) local incidence angle (between surface normal and look vector, enter - for none) psi (output) projection angle (between surface normal and image plane normal, enter - for none) pix (output) pixel area normalization factor (enter - for none) ls_map (output) layover and shadow map (in map projection, enter - for none) frame number of DEM pixels to add around area covered by SAR image (enter - for default = 8) ls_mode output lookup table values in regions of layover, shadow, or DEM gaps (enter - for default) 0: set to (0.,0.) 1: linear interpolation across these regions 2: actual value (default) 3: nn-thinned r_ovr range over-sampling factor for nn-thinned layover/shadow mode (enter - for default: 2.0)

---

NOTE: This script accepts the syntax of the initial version of gc_map (now named gc_map1) while calling gc_map2 for performing the calculations. gc_map2 produces more accurate layover and shadow maps, as well as incidence angle maps and its derived products such as the simulated SAR image. We recommend directly calling gc_map2. When prefered, the initial version of gc_map can still be used by calling gc_map1. *****`

[sixy6e]

I'm not quite following what the gc_map approach is, but the code for determining the lat and lon overview is taken directly from dem_ovr.

https://github.com/sixy6e/gamma_insar/blob/issue-46/insar/coregister_dem.py#L508

Although, we could've missed something in the BASH -> Python translation.

[tfuhrmann]

gc_map is a GAMMA function. I've just copied in the GAMMA docu in my above comment for your reference on where the dem_ovr factor is used.

[sixy6e]

Ah, ok.

[tfuhrmann]

I've added some lines of code to the existing bash script calc_multi-look_values.bash to calculate a DEM oversampling factor (see below, not pushed to the repo yet). For Sentinel-1 data in the Newcastle area the output of the script saved in ./results/T147D_multi-look_results would be as follows:

MULTI-LOOKING VALUE RESULTS

Input multi-look value: 2
Mean azimuth spacing [m]: 14.069718
Mean range spacing (ground) [m]: 3.702077
Azimuth multi-looking factor to obtain square pixels: -
Range multi-looking factor to obtain square pixels: 4

MLI range and azimuth looks: 8 2

DEM spacing latitude [m]: -30.8792853858
DEM spacing longitude [m]: 25.9156445572
Average DEM spacing [m]: 28.3974649715
Average spacing of multi-looked SAR images [m]: 28.878026
DEM oversampling factor suiting resolution of multi-looked SAR data: .983358

In this example, the oversampling factor is very close to 1 (which we are currently using as a hard-coded value for multi-looked Sentinel-1 data). The output of the script could either be rounded to the closest integer for the DEM oversampling or the float value can be directly used (to be tested if the DEM coregistration runs through without issues when using a float value).

Code added to calc_multi-look_values.bash (for discussion!):

[...]  
### calculate DEM oversampling factor that fits the resolution of mulit-looked SAR data

  # load DEM file names from gamma_functions
  dem_file_names
  # grep the posting in lat and lon, half posting to be added to offset (pixel centre)
  postlat=`grep post_lat: $dem_par | awk '{print $2}'`
  postlon=`grep post_lon: $dem_par | awk '{print $2}'`
  post_lat=`echo ${postlat} | sed -e 's/[eE]+*/\\*10\\^/'`
  post_lon=`echo ${postlon} | sed -e 's/[eE]+*/\\*10\\^/'`
  # also get the centre latitude
  corner_lat=`grep corner_lat: $dem_par | awk '{print $2}'`
  nlines=`grep nlines: $dem_par | awk '{print $2}'`
  centre_lat=`echo "scale=10 ; $corner_lat+($nlines/2*$post_lat)" | bc -l`

  # convert lat/lon posting to metric
  # ellipsoid constants (as given in DEM.par file)
  a=`grep ellipsoid_ra: $dem_par | awk '{print $2}'`
  rep_f=`grep ellipsoid_reciprocal_flattening: $dem_par | awk '{print $2}'`
  f=`echo "scale=10 ; 1/$rep_f" | bc -l`
  # calculate radii of curvature for degree -> metric conversion
  e2=`echo "scale=10 ; 1/(1-$f)^2-1" | bc -l`
  c=`echo "scale=10 ; $a*sqrt(1+$e2)" | bc -l`
  V=`echo "scale=10 ; sqrt((1+($e2*c($centre_lat/180*$pi)^2)))" | bc -l`
  N=`echo "scale=10 ; $c/$V" | bc -l`
  M=`echo "scale=10 ; $c/($V^3)" | bc -l`
  post_lat_m=`echo "scale=10 ; $post_lat/180*$pi*(sqrt($M*$N))" | bc -l`
  post_lon_m=`echo "scale=10 ; $post_lon/180*$pi*c($centre_lat/180*$pi)*sqrt($M*$N)" | bc -l`
  dem_post_m=`echo "scale=10 ; (${post_lat_m#-}+${post_lon_m#-})/2" | bc -l`
  echo "DEM spacing latitude [m]: "$post_lat_m >> $multi_results   
  echo "DEM spacing longitude [m]: "$post_lon_m >> $multi_results
  echo "Average DEM spacing [m]: "$dem_post_m >> $multi_results

  # calculate DEM oversampling factor from average spacing of SAR data and DEM
  ml_sar_sp_mean=`echo "scale=6 ; (($alks*$azspmean)+($rlks*$grrgspmean))/2" | bc -l`
  dem_ovr=`echo "scale=6 ; $dem_post_m/$ml_sar_sp_mean" | bc -l` 
  echo "Average spacing of multi-looked SAR images [m]: "$ml_sar_sp_mean >> $multi_results
  echo "DEM oversampling factor suiting resolution of multi-looked SAR data: "$dem_ovr >> $multi_results 

[sixy6e]

A quick Python equivalent, but missing the incident angle:

import math
import osr
import rasterio

DEM_FNAME = "/g/data/dg9/MASTER_DEM/GAMMA_DEM_SRTM_1as_mosaic.img" 

def ovr_factor(DEM_FNAME):
    # initialise a crs
    crs = osr.SpatialReference()

    with rasterio.open(DEM_FNAME) as src:
        lat_step = src.transform.e
        lon_step = src.transform.a

        # inverting to xy order
        centre_xy = [i // 2 for i in list(src.shape)[::-1]]
        centre_lon, centre_lat = centre_xy * src.transform

        # spheroidal parameters
        crs.ImportFromWkt(src.crs.wkt)
        a = crs.GetSemiMajor()
        rep_f = crs.GetInvFlattening()
        f = 1 / rep_f
        e2 = 1 / (1 - f)**2 - 1
        c = a * math.sqrt(1 + e2)
        V = math.sqrt((1 + (e2 * math.cos(math.radians(centre_lat))**2)))
        N = c / V
        M = c / V**3

        # y and x step size in metres
        step_lat_m = math.radians(lat_step) * (math.sqrt(M * N))
        step_lon_m = math.radians(lon_step) * math.cos(math.radians(centre_lat)) * math.sqrt(M * N)

        # from an SLC.par for Sentinel-1
        azsp = 14.066250
        rgsp = 9.318248

In this example, I used the DEM covering Australia, as it doesn't make sense to calculate a new ovr for every frame. It should be consistent. Also, the incidence angle, rather than finding an average for every SLC file, is there simply a nominal value that can be used. My reasoning is that there appears to be reverse logic here, in that we create the SLC files to then determine the nominal ovr factor. Doesn't the workflow require the ovr beforehand?

Optionally, this is a separate utility that is used to find an ovr to then be used in production and set in stone as part of some configuration file.

[tfuhrmann]

The spatial resolution of SLC data changes slightly with latitude. That's why we had placed the calculation of multi-looking values in azimuth and range direction after SLC creation. However, the effect is small enough that we can still go with one consist DEM oversampling factor for the whole of Australia, I think (see below example). Note that the range spacing given in GAMMA *.par files is the slant range (along the radar line-of-sight). The ground range spacing can be derived by dividing by sin(inc.).

Neglecting the latitude change in the SLC resolution for different frames, the DEM oversampling factor purely depends on the multi-looking factor. If the decision on the multi-looking factor is made at the very start of the processing, then the DEM ovr factor can also be calculated before starting the SLC creation, for example for ml=1 -> DEM_ovr=2, for ml=2 -> DEM_ovr=1 (for S1 data). In that case, the whole script calc_multi-look_values.bash is redundant and we can set the ml/ovr factors in stone for the S1 pipeline (e.g. ml_az=2, ml_rg=8, DEM_ovr=1).

comparison of multi-looked spatial resolution of S1 data for the northern NT and ACT region (we would need to have a look at Tassie data to check out the effect in the very South!) NT:

multi-looked azimuth spacing [m]: 28.12
multi-looked range spacing (ground) [m]: 29.33

ACT:

multi-looked azimuth spacing [m]: 28.14
multi-looked range spacing (ground) [m]: 29.72

[sixy6e]

I realise that projected resolutions differ with latitude, I was thinking that a consistent resolution is more ideal (at least it is in the optical RS world) for a continental product. Probably worth a quick telecon.

[tfuhrmann]

Ah yeah, that's what I wanted to raise as well (but forgot). The metric ground resolution of SRTM-1 data is latitude-dependent as well, since the posting is given in degrees (i.e. 1 arc-second). That's why my code example calculates the centre latitude of the DEM used for a particular frame. For northern Australia (-12 deg) the SRTM-1 ground resolution would be 29.91 m x 30.15 m, for Tassie (-44 deg) it would be 30.92 m x 22.24 m.

[mcgarth]

Decision: for NSW datasets multilook = 2 and dem_oversampling = 1