test_geometry.py
import shutil
from typing import Tuple
import numpy as np
from os.path import join
import pytest
import pyrate.constants as C
from pyrate.core import ifgconstants as ifc
from pyrate.core.geometry import get_lonlat_coords, get_sat_positions, vincinv
from pyrate.core.refpixel import convert_pixel_value_to_geographic_coordinate
from tests import common
from pyrate.configuration import Configuration
from subprocess import run, PIPE
from pyrate import prepifg, correct
from pyrate.core.shared import Ifg, Geometry
def get_lonlat_coords_slow(ifg: Ifg) -> Tuple[np.ndarray, np.ndarray]:
"""
Function to get longitude and latitude coordinates for each pixel in the multi-looked.
interferogram dataset. Coordinates are identical for each interferogram in the stack.
:param ifg: pyrate.core.shared.Ifg Class object.
:return: lon: Longitude for each pixel (decimal degrees)
:return: lat: Latitude for each pixel (decimal degrees)
"""
# assume all interferograms have same projection and will share the same transform
transform = ifg.dataset.GetGeoTransform()
# number of rows and columns in dataset
nrows, ncols = ifg.shape
lon = np.zeros((nrows, ncols)) # pre-allocate 2D numpy array
lat = np.zeros((nrows, ncols)) # pre-allocate 2D numpy array
for i in range(0, nrows): # rows are y-direction
for j in range(0, ncols): # cols are x-direction
lon[i, j], lat[i, j] = convert_pixel_value_to_geographic_coordinate(j, i, transform)
return lon, lat
def test_get_lonlat_coords_vectorised(dem):
lon, lat = get_lonlat_coords_slow(dem)
lon_v, lat_v = get_lonlat_coords(dem)
np.testing.assert_array_almost_equal(lon, lon_v.data)
np.testing.assert_array_almost_equal(lat, lat_v.data)
@pytest.fixture(params=[(29, 50, -2.94634866714478), (94, 58, -2.94684600830078)])
def point_azimuth(request):
return request.param
@pytest.fixture(params=[(29, 50, 1.02217936515808), (94, 58, 1.0111095905304)])
def point_incidence(request):
return request.param
class TestPyRateAngleFiles:
@classmethod
def setup_class(cls):
cls.params = Configuration(common.MEXICO_CROPA_CONF).__dict__
# run prepifg
prepifg.main(cls.params)
# copy IFGs to temp folder
correct._copy_mlooked(cls.params)
@classmethod
def teardown_class(cls):
shutil.rmtree(cls.params[C.OUT_DIR], ignore_errors=True)
# the xy position in the original GAMMA geometry (before cropping and further multi-looking is required for
# comparison of PyRate with GAMMA angle values. To get this position for a particular coordinate, the following
# commands can be applied on Gadi:
# cd /g/data/dg9/INSAR_ANALYSIS/MEXICO/S1/UNIT-TEST-DATA/CROP-A/GEOTIFF
# gdallocationinfo -wgs84 cropA_20180106-20180130_VV_8rlks_eqa_unw.tif -99.06 19.37
# -> outputs the cropped pixel coordinates: (94P,58L) -> x0 = 94, y0 = 58
# cd /g/data/dg9/INSAR_ANALYSIS/MEXICO/S1/GAMMA/T005A/geotiff_files/unw_ifg_geotiffs/
# gdallocationinfo -wgs84 20180106-20180130_VV_8rlks_eqa_unw.tif -99.06 19.37
# -> outputs the original pixel coordinates: (940P,3271L)
# cd /g/data/dg9/INSAR_ANALYSIS/MEXICO/S1/GAMMA/T005A/DEM/
# gdallocationinfo 20180106_VV_8rlks_eqa_lv_phi.tif 940 3271
# -> outputs the GAMMA azimuth angle: -2.94684600830078
# gdallocationinfo 20180106_VV_8rlks_eqa_lv_theta.tif 940 3271
# > outputs the GAMMA incidence angle: 1.0111095905304
def get_pyrate_angle(self, x0, y0, tif_file):
"""
Get angle at particular pixel in the azimuth/incidence tif file
"""
# get azimuth angle value of PyRate file azimuth_angle.tif
temp = run(['gdallocationinfo', tif_file, str(x0), str(y0)], universal_newlines=True, stdout=PIPE)
out = temp.stdout
angle_value = float(out.split('Value:')[1].strip('\n'))
return angle_value
def test_pyrate_azimuth_matches_gamma_azimuth(self, point_azimuth):
x0, y0, exp = point_azimuth
# azimuth angle validation value extracted from GAMMA azimuth file for pixel location corresponding
# to (29,50) using gdallocationinfo 20180106_VV_8rlks_eqa_lv_phi.tif 613 3228
# exp = -2.94634866714478 # GAMMMA azimuth angle at pixel location
# azimuth angle validation value extracted from GAMMA azimuth file for pixel location corresponding
# to (94,58) using gdallocationinfo 20180106_VV_8rlks_eqa_lv_phi.tif 940 3271
# exp = -2.94684600830078
# GAMMA azimuth is defined towards satellite in an anti-clockwise angle system, with East being zero
# PyRate azimuth angle is defined towards the satellite in a clockwise angle system with North being zero
tif_file = join(self.params[C.GEOMETRY_DIR], 'azimuth_angle.tif')
azimuth_angle_pyrate = self.get_pyrate_angle(x0, y0, tif_file)
# convert PyRate azimuth into GAMMA azimuth
res = -(azimuth_angle_pyrate - np.pi / 2)
np.testing.assert_array_almost_equal(exp, res, decimal=2) # max difference < 0.01 rad
# screen output of difference if need be:
# print(exp - res)
def test_pyrate_incidence_matches_gamma_incidence(self, point_incidence):
x0, y0, exp = point_incidence
# incidence angle validation value extracted from GAMMA incidence file for pixel location corresponding
# to (29,50) using gdallocationinfo 20180106_VV_8rlks_eqa_lv_theta.tif 613 3228
# exp = 1.02217936515808
# incidence angle validation value extracted from GAMMA incidence file for pixel location corresponding
# to (94,58) using gdallocationinfo 20180106_VV_8rlks_eqa_lv_theta.tif 940 3271
# exp = 1.0111095905304
# GAMMA angle is defined from the horizontal plane with the zenith direction being pi / 2 radians (i.e.90 deg)
# PyRate angle is defined from the vertical axis with the zenith direction being zero
tif_file = join(self.params[C.GEOMETRY_DIR], 'incidence_angle.tif')
incidence_angle_pyrate = self.get_pyrate_angle(x0, y0, tif_file)
# convert PyRate incidence into GAMMA incidence
res = np.pi / 2 - incidence_angle_pyrate
# convert PyRate incidence into GAMMA incidence
np.testing.assert_array_almost_equal(exp, res, decimal=3) # max difference < 0.001 rad
# screen output of difference if need be:
# print(exp - res)
def test_azimuth_angle_calculation(self):
"""
Calculate local azimuth angle using a spherical model and compare to result using Vincenty's equations
"""
# get first IFG in stack to calculate lon/lat values
multi_paths = self.params[C.INTERFEROGRAM_FILES]
tmp_paths = [ifg_path.tmp_sampled_path for ifg_path in multi_paths]
# keep only ifg files in path list (i.e. remove coherence and dem files)
ifg_paths = [item for item in tmp_paths if 'ifg.tif' in item]
# read and open the first IFG in list
ifg0_path = ifg_paths[0]
ifg0 = Ifg(ifg0_path)
ifg0.open(readonly=True)
lon, lat = get_lonlat_coords(ifg0)
# read incidence and look angle files
tif_file = join(self.params[C.GEOMETRY_DIR], 'incidence_angle.tif')
geom = Geometry(tif_file)
incidence_angle = geom.data
tif_file = join(self.params[C.GEOMETRY_DIR], 'look_angle.tif')
geom = Geometry(tif_file)
look_angle = geom.data
# get metadata
a = float(ifg0.meta_data[ifc.PYRATE_SEMI_MAJOR_AXIS_METRES])
b = float(ifg0.meta_data[ifc.PYRATE_SEMI_MINOR_AXIS_METRES])
heading = float(ifg0.meta_data[ifc.PYRATE_HEADING_DEGREES])
azimuth = float(ifg0.meta_data[ifc.PYRATE_AZIMUTH_DEGREES])
# convert all angles from deg to radians
lon = np.radians(lon.data)
lat = np.radians(lat.data)
heading = np.radians(heading)
azimuth = np.radians(azimuth)
# calculate satellite positions
sat_lat, sat_lon = get_sat_positions(lat, lon, look_angle, incidence_angle, heading, azimuth)
# calculate azimuth angle from pixel to satellite using spherical trigonometry
# see Eq. 86 on page 4-11 in EARTH-REFERENCED AIRCRAFT NAVIGATION AND SURVEILLANCE ANALYSIS
# (https://ntlrepository.blob.core.windows.net/lib/59000/59300/59358/DOT-VNTSC-FAA-16-12.pdf)
azimuth_angle_spherical = np.arctan2(np.cos(sat_lat) * np.sin(sat_lon - lon), \
np.sin(sat_lat) * np.cos(lat) - np.cos(sat_lat) * np.sin(lat) * \
np.cos(sat_lon - lon))
# add 2 pi in case an angle is below zero
for azi in np.nditer(azimuth_angle_spherical, op_flags=['readwrite']):
if azi < 0:
azi[...] = azi + 2 * np.pi
azimuth_angle_elliposidal = vincinv(lat, lon, sat_lat, sat_lon, a, b)
# the difference between Vincenty's azimuth calculation and the spherical approximation is ~0.001 radians
azimuth_angle_diff = azimuth_angle_spherical - azimuth_angle_elliposidal
# max difference < 0.01 rad
np.testing.assert_array_almost_equal(azimuth_angle_spherical, azimuth_angle_elliposidal, decimal=2)
