Unverified Commit 87d32c56 authored by Jesse Mapel's avatar Jesse Mapel Committed by GitHub
Browse files

Added ISD generation notebooks (#282)

* Added notebooks

* Updated SAR notebook
parent 70eda0af
Loading
Loading
Loading
Loading
+97 −0
Original line number Original line Diff line number Diff line
%% Cell type:code id: tags:

``` python
import numpy as np
np.set_printoptions(suppress=True)
import pyproj
```

%% Output

    ---------------------------------------------------------------------------
    ModuleNotFoundError                       Traceback (most recent call last)
    <ipython-input-1-7e71c9f92a96> in <module>
          1 import numpy as np
          2 np.set_printoptions(suppress=True)
    ----> 3 import pyproj

    ModuleNotFoundError: No module named 'pyproj'

%% Cell type:code id: tags:

``` python
def reproject(record, semi_major, semi_minor, source_proj, dest_proj, **kwargs):
    """
    Thin wrapper around PyProj's Transform() function to transform 1 or more three-dimensional
    point from one coordinate system to another. If converting between Cartesian
    body-centered body-fixed (BCBF) coordinates and Longitude/Latitude/Altitude coordinates,
    the values input for semi-major and semi-minor axes determine whether latitudes are
    planetographic or planetocentric and determine the shape of the datum for altitudes.
    If semi_major == semi_minor, then latitudes are interpreted/created as planetocentric
    and altitudes are interpreted/created as referenced to a spherical datum.
    If semi_major != semi_minor, then latitudes are interpreted/created as planetographic
    and altitudes are interpreted/created as referenced to an ellipsoidal datum.
    Parameters
    ----------
    record : object
          Pandas series object
    semi_major : float
              Radius from the center of the body to the equater
    semi_minor : float
              Radius from the pole to the center of mass
    source_proj : str
                      Pyproj string that defines a projection space ie. 'geocent'
    dest_proj : str
                   Pyproj string that defines a project space ie. 'latlon'
    Returns
    -------
    : list
    Transformed coordinates as y, x, z
    """
    source_pyproj = pyproj.Proj(proj = source_proj, a = semi_major, b = semi_minor)
    dest_pyproj = pyproj.Proj(proj = dest_proj, a = semi_major, b = semi_minor)

    y, x, z = pyproj.transform(source_pyproj, dest_pyproj, record[0], record[1], record[2], **kwargs)

    return y, x, z
```

%% Cell type:markdown id: tags:

## Input Data

%% Cell type:code id: tags:

``` python
semi_major = 1100
semi_minor = 1000
altitude = 50
lon = 35
lat = 30
radius = 0
```

%% Cell type:markdown id: tags:

## Compute the quaternion (as x, y, z, w)

%% Cell type:code id: tags:

``` python
ground_pt = np.asarray(reproject([lon, lat, radius], semi_major, semi_minor, 'latlong', 'geocent'))
sensor_position = np.asarray(reproject([lon, lat, radius + altitude], semi_major, semi_minor, 'latlong', 'geocent'))
look_vec = ground_pt - sensor_position
look_vec = look_vec / np.linalg.norm(look_vec)
quat = np.zeros(4)
quat[:3] = np.cross(np.array([0, 0, 1]), look_vec)
quat[3] = np.dot(np.array([0, 0, 1]), look_vec)
```

%% Output

    [ 0.49673176 -0.70940648  0.         -0.5       ]

%% Cell type:code id: tags:

``` python
```
+96 −0
Original line number Original line Diff line number Diff line
%% Cell type:code id: tags:

``` python
import numpy as np
np.set_printoptions(suppress=True)
```

%% Cell type:code id: tags:

``` python
radius = 10 # km
altitude = 9990 # km
detector_size = 0.1 # mm
focal_length = 50 # mm
lines = 16
exposure_duration = 0.1 # seconds
```

%% Cell type:code id: tags:

``` python
positions = []
velocities = []
quats = []
for i in np.arange(0,16):
    angle = i*detector_size/(radius+altitude)
    position = (radius+altitude) * np.array([np.cos(angle), 0, -np.sin(angle)])
    positions.append(position)
    velocity = (radius+altitude) * np.array([-np.sin(angle), 0, -np.cos(angle)])
    velocities.append(velocity)
    camera_angle = -np.pi/2 + angle
    quat = np.array([0, np.sin(camera_angle/2), 0, np.cos(camera_angle/2)])
    quats.append(quat)
for pos in positions:
    print(pos)
for vel in velocities:
    print(vel)
for quat in quats:
    print(quat)
```

%% Output

    [1050.    0.   -0.]
    [1049.99999524    0.           -0.1       ]
    [1049.99998095    0.           -0.2       ]
    [1049.99995714    0.           -0.3       ]
    [1049.99992381    0.           -0.39999999]
    [1049.99988095    0.           -0.49999998]
    [1049.99982857    0.           -0.59999997]
    [1049.99976667    0.           -0.69999995]
    [1049.99969524    0.           -0.79999992]
    [1049.99961429    0.           -0.89999989]
    [1049.99952381    0.           -0.99999985]
    [1049.99942381    0.           -1.0999998 ]
    [1049.99931429    0.           -1.19999974]
    [1049.99919524    0.           -1.29999967]
    [1049.99906667    0.           -1.39999959]
    [1049.99892857    0.           -1.49999949]
    [   -0.     0. -1050.]
    [   -0.1            0.         -1049.99999524]
    [   -0.2            0.         -1049.99998095]
    [   -0.3            0.         -1049.99995714]
    [   -0.39999999     0.         -1049.99992381]
    [   -0.49999998     0.         -1049.99988095]
    [   -0.59999997     0.         -1049.99982857]
    [   -0.69999995     0.         -1049.99976667]
    [   -0.79999992     0.         -1049.99969524]
    [   -0.89999989     0.         -1049.99961429]
    [   -0.99999985     0.         -1049.99952381]
    [   -1.0999998      0.         -1049.99942381]
    [   -1.19999974     0.         -1049.99931429]
    [   -1.29999967     0.         -1049.99919524]
    [   -1.39999959     0.         -1049.99906667]
    [   -1.49999949     0.         -1049.99892857]
    [ 0.         -0.70710678  0.          0.70710678]
    [ 0.         -0.70707311  0.          0.70714045]
    [ 0.         -0.70703943  0.          0.70717412]
    [ 0.         -0.70700576  0.          0.70720779]
    [ 0.         -0.70697208  0.          0.70724146]
    [ 0.         -0.7069384   0.          0.70727512]
    [ 0.         -0.70690472  0.          0.70730878]
    [ 0.         -0.70687104  0.          0.70734244]
    [ 0.         -0.70683736  0.          0.7073761 ]
    [ 0.         -0.70680367  0.          0.70740976]
    [ 0.         -0.70676998  0.          0.70744342]
    [ 0.         -0.70673629  0.          0.70747707]
    [ 0.         -0.7067026   0.          0.70751073]
    [ 0.         -0.70666891  0.          0.70754438]
    [ 0.         -0.70663522  0.          0.70757803]
    [ 0.         -0.70660152  0.          0.70761168]

%% Cell type:code id: tags:

``` python
```
+175 −0
Original line number Original line Diff line number Diff line
%% Cell type:code id: tags:

``` python
import numpy as np
from itertools import product
```

%% Cell type:code id: tags:

``` python
radius   = 1737400
alt      = 2000000
ground   = 7.5
exposure = 0.005
samples  = 1000
lines    = 1000
```

%% Cell type:code id: tags:

``` python
sensor_rad = radius + alt
angle_per_line = ground / radius
angle_per_samp = angle_per_line
angle_per_Second = angle_per_line / exposure
```

%% Cell type:code id: tags:

``` python
line_vec = np.arange(0, lines+0.00000001)
sample_vec = np.arange(0, samples+0.00000001)
# From here on, matrix indexing is [line, sample, (xyz)]
sample_mat, line_mat = np.meshgrid(line_vec, sample_vec)
```

%% Cell type:code id: tags:

``` python
positions = sensor_rad * np.vstack([np.cos(-angle_per_line * line_vec), np.zeros(line_vec.shape), np.sin(-angle_per_line * line_vec)]).T
# Note: The chain rule results in an extra negative on the velocity calculations
velocities = sensor_rad * np.vstack([np.sin(-angle_per_line * line_vec), np.zeros(line_vec.shape), -np.cos(-angle_per_line * line_vec)]).T
print('Positions')
for pos in positions[::int(0.25/exposure)]:
    print(list(pos))
print('Velocities')
for vel in velocities[::int(0.25/exposure)]:
    print(list(vel))
```

%% Output

    Positions
    [3737400.0, 0.0, -0.0]
    [3737399.912943243, 0.0, -806.6795148600813]
    [3737399.651772976, 0.0, -1613.3589921395437]
    [3737399.216489211, 0.0, -2420.03839425777]
    [3737398.607091969, 0.0, -3226.7176836341473]
    [3737397.823581278, 0.0, -4033.396822688066]
    [3737396.8659571735, 0.0, -4840.075773838925]
    [3737395.734219702, 0.0, -5646.754499506133]
    [3737394.4283689144, 0.0, -6453.432962109106]
    [3737392.9484048723, 0.0, -7260.111124067275]
    [3737391.2943276456, 0.0, -8066.788947800085]
    [3737389.4661373096, 0.0, -8873.466395726995]
    [3737387.4638339505, 0.0, -9680.14343026748]
    [3737385.287417662, 0.0, -10486.820013841041]
    [3737382.936888544, 0.0, -11293.496108867193]
    [3737380.4122467074, 0.0, -12100.17167776548]
    [3737377.713492269, 0.0, -12906.846682955462]
    [3737374.840625355, 0.0, -13713.521086856732]
    [3737371.7936460995, 0.0, -14520.194851888911]
    [3737368.5725546437, 0.0, -15326.867940471646]
    [3737365.177351138, 0.0, -16133.540315024618]
    Velocities
    [-0.0, 0.0, -3737400.0]
    [-806.6795148600813, 0.0, -3737399.912943243]
    [-1613.3589921395437, 0.0, -3737399.651772976]
    [-2420.03839425777, 0.0, -3737399.216489211]
    [-3226.7176836341473, 0.0, -3737398.607091969]
    [-4033.396822688066, 0.0, -3737397.823581278]
    [-4840.075773838925, 0.0, -3737396.8659571735]
    [-5646.754499506133, 0.0, -3737395.734219702]
    [-6453.432962109106, 0.0, -3737394.4283689144]
    [-7260.111124067275, 0.0, -3737392.9484048723]
    [-8066.788947800085, 0.0, -3737391.2943276456]
    [-8873.466395726995, 0.0, -3737389.4661373096]
    [-9680.14343026748, 0.0, -3737387.4638339505]
    [-10486.820013841041, 0.0, -3737385.287417662]
    [-11293.496108867193, 0.0, -3737382.936888544]
    [-12100.17167776548, 0.0, -3737380.4122467074]
    [-12906.846682955462, 0.0, -3737377.713492269]
    [-13713.521086856732, 0.0, -3737374.840625355]
    [-14520.194851888911, 0.0, -3737371.7936460995]
    [-15326.867940471646, 0.0, -3737368.5725546437]
    [-16133.540315024618, 0.0, -3737365.177351138]

%% Cell type:code id: tags:

``` python
lat = -angle_per_line * line_mat
# Image is a right look, so the longitude goes negative
lon = -angle_per_samp * sample_mat
ground_points = radius * np.stack([np.multiply(np.cos(lat), np.cos(lon)), np.multiply(np.cos(lat), np.sin(lon)), np.sin(lat)], axis=-1)
print("Ground point at line: 500, sample: 500")
print(ground_points[500, 500])
```

%% Output

    Ground point at line: 500, sample: 500
    [1737391.90602155   -3749.98835331   -3749.99708833]

%% Cell type:code id: tags:

``` python
slant_range = np.array([[np.linalg.norm(point) for point in row] for row in ground_points - positions[:, None, :]])
ground_range = radius * np.abs(lon)
```

%% Cell type:code id: tags:

``` python
range_poly = np.polynomial.polynomial.polyfit(ground_range.flatten(), slant_range.flatten(), 3)
print("Ground range to slant range polynomial coefficients")
print(list(range_poly))
```

%% Output

    Ground range to slant range polynomial coefficients
    [2000000.000003915, -1.0488420462327845e-08, 5.377893056639776e-07, -1.3072058387193456e-15]

%% Cell type:code id: tags:

``` python
sc_pos = positions[500]
sc_vel = velocities[500]
off_ground_pt = ground_points[500, 500] - np.array([100, 100, 100])
look_vec = off_ground_pt - positions[500]
zero_doppler_look_vec = look_vec - np.dot(look_vec, sc_vel) / np.dot(sc_vel, sc_vel) * sc_vel
locus_point = sc_pos + slant_range[500, 500] / np.linalg.norm(zero_doppler_look_vec) * zero_doppler_look_vec
# Image is a right look, so do look X velocity
locus_direction = np.cross(zero_doppler_look_vec, sc_vel)
locus_direction = 1.0 / np.linalg.norm(locus_direction) * locus_direction
print("Locus point:", list(locus_point))
print("Locus direction:", list(locus_direction))
```

%% Output

    Locus point: [1737392.0956685692, -3849.7959147875467, -3749.9887626446034]
    Locus direction: [0.0019248861951120758, -0.9999981473962212, -4.154676206387554e-06]

%% Cell type:code id: tags:

``` python
remote_look_vec = -slant_range[500, 500] / sensor_rad * sc_pos
remote_zero_doppler_look_vec = remote_look_vec - np.dot(remote_look_vec, sc_vel) / np.dot(sc_vel, sc_vel) * sc_vel
remote_locus_point = sc_pos + remote_zero_doppler_look_vec
remote_locus_direction = np.cross(remote_zero_doppler_look_vec, sc_vel)
remote_locus_direction = 1.0 / np.linalg.norm(remote_locus_direction) * remote_locus_direction
print("Remote locus point:", list(remote_locus_point))
print("Remote locus direction:", list(remote_locus_direction))
```

%% Output

    Remote locus point: [1737388.3904556318, 0.0, -3749.980765309453]
    Remote locus direction: [-0.0, -1.0, 0.0]

%% Cell type:code id: tags:

``` python
```