1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186 | 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)
|