Source code for pyrocko.scenario.targets.insar

# http://pyrocko.org - GPLv3
#
# The Pyrocko Developers, 21st Century
# ---|P------/S----------~Lg----------
from __future__ import absolute_import, division, print_function

import numpy as num
import logging
from os import path as op

from datetime import datetime
from pyrocko import gf, util
from pyrocko import orthodrome as od
from pyrocko.guts import Float, Timestamp, Tuple, StringChoice, Bool, Object,\
    String
from pyrocko.util import num_full

from .base import TargetGenerator, NoiseGenerator
from ..base import get_gsshg

DEFAULT_STORE_ID = 'ak135_static'

km = 1e3
d2r = num.pi/180.

logger = logging.getLogger('pyrocko.scenario.targets.insar')
guts_prefix = 'pf.scenario'


[docs]class ScenePatch(Object): lat_center = Float.T( help='Center latitude anchor.') lon_center = Float.T( help='center longitude anchor.') time_master = Timestamp.T( help='Timestamp of the master.') time_slave = Timestamp.T( help='Timestamp of the slave.') inclination = Float.T( help='Orbital inclination towards the equatorial plane [deg].') apogee = Float.T( help='Apogee of the satellite in [m].') swath_width = Float.T( help='Swath width in [m].') track_length = Float.T( help='Track length in [m].') incident_angle = Float.T( help='Ground incident angle in [deg].') resolution = Tuple.T( help='Resolution of raster in east x north [px].') orbital_node = StringChoice.T( ['Ascending', 'Descending'], help='Orbit heading.') mask_water = Bool.T( default=True, help='Mask water bodies.')
[docs] class SatelliteGeneratorTarget(gf.SatelliteTarget): def __init__(self, scene_patch, *args, **kwargs): gf.SatelliteTarget.__init__(self, *args, **kwargs) self.scene_patch = scene_patch def post_process(self, *args, **kwargs): resp = gf.SatelliteTarget.post_process(self, *args, **kwargs) from kite import Scene from kite.scene import SceneConfig, FrameConfig, Meta patch = self.scene_patch grid, _ = patch.get_grid() displacement = num.empty_like(grid) displacement.fill(num.nan) displacement[patch.get_mask()] = resp.result['displacement.los'] theta, phi = patch.get_incident_angles() llLat, llLon = patch.get_ll_anchor() urLat, urLon = patch.get_ur_anchor() dLon = num.abs(llLon - urLon) / patch.resolution[0] dLat = num.abs(llLat - urLat) / patch.resolution[1] scene_config = SceneConfig( meta=Meta( scene_title='Pyrocko Scenario Generator - {orbit} ({time})' .format(orbit=self.scene_patch.orbital_node, time=datetime.now()), orbital_node=patch.orbital_node, scene_id='pyrocko_scenario_%s' % self.scene_patch.orbital_node, satellite_name='Sentinel-1 (pyrocko-scenario)'), frame=FrameConfig( llLat=float(llLat), llLon=float(llLon), dN=float(dLat), dE=float(dLon), spacing='degree')) scene = Scene( displacement=displacement, theta=theta, phi=phi, config=scene_config) resp.scene = scene return resp
def __init__(self, *args, **kwargs): Object.__init__(self, *args, **kwargs) self._mask_water = None @property def width(self): track_shift = num.abs( num.cos(self.inclination*d2r) * self.track_length) return self.swath_width + track_shift def get_ll_anchor(self): return od.ne_to_latlon(self.lat_center, self.lon_center, -self.track_length/2, -self.width/2) def get_ur_anchor(self): return od.ne_to_latlon(self.lat_center, self.lon_center, self.track_length/2, self.width/2) def get_ul_anchor(self): return od.ne_to_latlon(self.lat_center, self.lon_center, self.track_length/2, -self.width/2) def get_lr_anchor(self): return od.ne_to_latlon(self.lat_center, self.lon_center, -self.track_length/2, self.width/2) def get_corner_coordinates(self): inc = self.inclination llLat, llLon = self.get_ll_anchor() urLat, urLon = self.get_ur_anchor() if self.orbital_node == 'Ascending': ulLat, ulLon = od.ne_to_latlon( self.lat_center, self.lon_center, self.track_length/2, -num.tan(inc*d2r) * self.width/2) lrLat, lrLon = od.ne_to_latlon( self.lat_center, self.lon_center, -self.track_length/2, num.tan(inc*d2r) * self.width/2) elif self.orbital_node == 'Descending': urLat, urLon = od.ne_to_latlon( self.lat_center, self.lon_center, self.track_length/2, num.tan(inc*d2r) * self.width/2) llLat, llLon = od.ne_to_latlon( self.lat_center, self.lon_center, -self.track_length/2, -num.tan(inc*d2r) * self.width/2) return ((llLat, llLon), (ulLat, ulLon), (urLat, urLon), (lrLat, lrLon))
[docs] def get_grid(self): '''Return relative positions of scatterer. :param track: Acquisition track, from `'asc'` or `'dsc'`. :type track: string ''' easts = num.linspace(0, self.width, self.resolution[0]) norths = num.linspace(0, self.track_length, self.resolution[1]) return num.meshgrid(easts, norths)
def get_mask_track(self): east_shifts, north_shifts = self.get_grid() norths = north_shifts[:, 0] track = num.abs(num.cos(self.inclination*d2r)) * norths track_mask = num.logical_and( east_shifts > track[:, num.newaxis], east_shifts < (track + self.swath_width)[:, num.newaxis]) if self.orbital_node == 'Ascending': track_mask = num.fliplr(track_mask) return track_mask def get_mask_water(self): if self._mask_water is None: east_shifts, north_shifts = self.get_grid() east_shifts -= east_shifts[0, -1]/2 north_shifts -= north_shifts[-1, -1]/2 latlon = od.ne_to_latlon(self.lat_center, self.lon_center, north_shifts.ravel(), east_shifts.ravel()) points = num.array(latlon).T self._mask_water = get_gsshg().get_land_mask(points)\ .reshape(*east_shifts.shape) return self._mask_water def get_mask(self): mask_track = self.get_mask_track() if self.mask_water: mask_water = self.get_mask_water() return num.logical_and(mask_track, mask_water) return mask_track def get_incident_angles(self): # theta: elevation angle towards satellite from horizon in radians. # phi: Horizontal angle towards satellite' :abbr:`line of sight (LOS)` # in [rad] from East. east_shifts, _ = self.get_grid() phi = num.empty_like(east_shifts) theta = num.empty_like(east_shifts) east_shifts += num.tan(self.incident_angle*d2r) * self.apogee theta = num.arctan(east_shifts/self.apogee) if self.orbital_node == 'Ascending': phi.fill(self.inclination*d2r + num.pi/2) elif self.orbital_node == 'Descending': phi.fill(2*num.pi-(self.inclination*d2r + 3/2*num.pi)) theta = num.fliplr(theta) else: raise AttributeError( 'Orbital node %s not defined!' % self.orbital_node) theta[~self.get_mask()] = num.nan phi[~self.get_mask()] = num.nan return theta, phi def get_target(self): gE, gN = self.get_grid() mask = self.get_mask() east_shifts = gE[mask].ravel() north_shifts = gN[mask].ravel() llLat, llLon = self.get_ll_anchor() ncoords = east_shifts.size theta, phi = self.get_incident_angles() theta = theta[mask].ravel() phi = phi[mask].ravel() if ncoords == 0: logger.warning('InSAR taget has no valid points,' ' maybe it\'s all water?') return self.SatelliteGeneratorTarget( scene_patch=self, lats=num_full(ncoords, fill_value=llLat), lons=num_full(ncoords, fill_value=llLon), east_shifts=east_shifts, north_shifts=north_shifts, theta=theta, phi=phi)
[docs]class AtmosphericNoiseGenerator(NoiseGenerator): amplitude = Float.T( default=1., help='Amplitude of the atmospheric noise.') beta = [5./3, 8./3, 2./3] regimes = [.15, .99, 1.] def get_noise(self, scene): nE = scene.frame.rows nN = scene.frame.cols if (nE+nN) % 2 != 0: raise ArithmeticError('Dimensions of synthetic scene must ' 'both be even!') dE = scene.frame.dE dN = scene.frame.dN rfield = num.random.rand(nE, nN) spec = num.fft.fft2(rfield) kE = num.fft.fftfreq(nE, dE) kN = num.fft.fftfreq(nN, dN) k_rad = num.sqrt(kN[:, num.newaxis]**2 + kE[num.newaxis, :]**2) regimes = num.array(self.regimes) k0 = 0. k1 = regimes[0] * k_rad.max() k2 = regimes[1] * k_rad.max() r0 = num.logical_and(k_rad > k0, k_rad < k1) r1 = num.logical_and(k_rad >= k1, k_rad < k2) r2 = k_rad >= k2 beta = num.array(self.beta) # From Hanssen (2001) # beta+1 is used as beta, since, the power exponent # is defined for a 1D slice of the 2D spectrum: # austin94: "Adler, 1981, shows that the surface profile # created by the intersection of a plane and a # 2-D fractal surface is itself fractal with # a fractal dimension equal to that of the 2D # surface decreased by one." beta += 1. # From Hanssen (2001) # The power beta/2 is used because the power spectral # density is proportional to the amplitude squared # Here we work with the amplitude, instead of the power # so we should take sqrt( k.^beta) = k.^(beta/2) RH # beta /= 2. amp = num.zeros_like(k_rad) amp[r0] = k_rad[r0] ** -beta[0] amp[r0] /= amp[r0].max() amp[r1] = k_rad[r1] ** -beta[1] amp[r1] /= amp[r1].max() / amp[r0].min() amp[r2] = k_rad[r2] ** -beta[2] amp[r2] /= amp[r2].max() / amp[r1].min() amp[k_rad == 0.] = amp.max() spec *= self.amplitude * num.sqrt(amp) noise = num.abs(num.fft.ifft2(spec)) noise -= num.mean(noise) return noise
[docs]class InSARGenerator(TargetGenerator): # https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath store_id = String.T( default=DEFAULT_STORE_ID, help='Store ID for these stations.') inclination = Float.T( default=98.2, help='Inclination of the satellite orbit towards equatorial plane' ' in [deg]. Defaults to Sentinel-1 (98.1 deg).') apogee = Float.T( default=693.*km, help='Apogee of the satellite in [m]. ' 'Defaults to Sentinel-1 (693 km).') swath_width = Float.T( default=250*km, help='Swath width in [m]. ' 'Defaults to Sentinel-1 Interfeometric Wide Swath Mode (IW).' ' (IW; 250 km).') track_length = Float.T( default=250*km, help='Track length in [m]. Defaults to 200 km.') incident_angle = Float.T( default=29.1, help='Near range incident angle in [deg]. Defaults to 29.1 deg;' ' Sentinel IW mode (29.1 - 46.0 deg).') resolution = Tuple.T( default=(250, 250), help='Resolution of raster in east x north [px].') mask_water = Bool.T( default=True, help='Mask out water bodies.') noise_generator = NoiseGenerator.T( default=AtmosphericNoiseGenerator.D(), help='Add atmospheric noise model after Hansen, 2001.') def get_scene_patches(self): lat_center, lon_center = self.get_center_latlon() scene_patches = [] for direction in ('Ascending', 'Descending'): patch = ScenePatch( lat_center=lat_center, lon_center=lon_center, time_master=0, time_slave=0, inclination=self.inclination, apogee=self.apogee, swath_width=self.swath_width, track_length=self.track_length, incident_angle=self.incident_angle, resolution=self.resolution, orbital_node=direction, mask_water=self.mask_water) scene_patches.append(patch) return scene_patches
[docs] def get_targets(self): targets = [s.get_target() for s in self.get_scene_patches()] for t in targets: t.store_id = self.store_id return targets
def get_insar_scenes(self, engine, sources, tmin=None, tmax=None): logger.debug('Forward modelling InSAR displacement...') scenario_tmin, scenario_tmax = self.get_time_range(sources) try: resp = engine.process( sources, self.get_targets(), nthreads=0) except gf.meta.OutOfBounds: logger.warning('Could not calculate InSAR displacements' ' - the GF store\'s extend is too small!') return [] scenes = [res.scene for res in resp.static_results()] tmin, tmax = self.get_time_range(sources) for sc in scenes: sc.meta.time_master = util.to_time_float(tmin) sc.meta.time_slave = util.to_time_float(tmax) scenes_asc = [sc for sc in scenes if sc.config.meta.orbital_node == 'Ascending'] scenes_dsc = [sc for sc in scenes if sc.config.meta.orbital_node == 'Descending'] def stack_scenes(scenes): base = scenes[0] for sc in scenes[1:]: base += sc return base scene_asc = stack_scenes(scenes_asc) scene_dsc = stack_scenes(scenes_dsc) if self.noise_generator: scene_asc.displacement += self.noise_generator.get_noise(scene_asc) scene_dsc.displacement += self.noise_generator.get_noise(scene_dsc) return scene_asc, scene_dsc def ensure_data(self, engine, sources, path, tmin=None, tmax=None): path_insar = op.join(path, 'insar') util.ensuredir(path_insar) tmin, tmax = self.get_time_range(sources) tts = util.time_to_str fn_tpl = op.join(path_insar, 'colosseo-scene-{orbital_node}_%s_%s' % (tts(tmin, '%Y-%m-%d'), tts(tmax, '%Y-%m-%d'))) def scene_fn(track): return fn_tpl.format(orbital_node=track.lower()) for track in ('ascending', 'descending'): fn = '%s.yml' % scene_fn(track) if op.exists(fn): logger.debug('Scene exists: %s' % fn) continue scenes = self.get_insar_scenes(engine, sources, tmin, tmax) for sc in scenes: fn = scene_fn(sc.config.meta.orbital_node) logger.debug('Writing %s' % fn) sc.save('%s.npz' % fn) def add_map_artists(self, engine, sources, automap): logger.warning('InSAR mapping is not implemented!') return None