|
# http://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
String
help='Center latitude anchor.') help='center longitude anchor.') help='Timestamp of the master.') help='Timestamp of the slave.') help='Orbital inclination towards the equatorial plane [deg].') help='Apogee of the satellite in [m].') help='Swath width in [m].') help='Track length in [m].') help='Ground incident angle in [deg].') help='Resolution of raster in east x north [px].') ['Ascending', 'Descending'], help='Orbit heading.') default=True, help='Mask water bodies.')
gf.SatelliteTarget.__init__(self, *args, **kwargs)
self.scene_patch = scene_patch
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
Object.__init__(self, *args, **kwargs) self._mask_water = None
def width(self): track_shift = num.abs( num.cos(self.inclination*d2r) * self.track_length)
return self.swath_width + track_shift
return od.ne_to_latlon(self.lat_center, self.lon_center, -self.track_length/2, -self.width/2)
return od.ne_to_latlon(self.lat_center, self.lon_center, self.track_length/2, self.width/2)
return od.ne_to_latlon(self.lat_center, self.lon_center, self.track_length/2, -self.width/2)
return od.ne_to_latlon(self.lat_center, self.lon_center, -self.track_length/2, self.width/2)
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))
'''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)
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
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
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
# 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
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)
default=1., help='Amplitude of the atmospheric noise.')
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
# https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/acquisition-modes/interferometric-wide-swath default=DEFAULT_STORE_ID, help='Store ID for these stations.')
default=98.2, help='Inclination of the satellite orbit towards equatorial plane' ' in [deg]. Defaults to Sentinel-1 (98.1 deg).') default=693.*km, help='Apogee of the satellite in [m]. ' 'Defaults to Sentinel-1 (693 km).') default=250*km, help='Swath width in [m]. ' 'Defaults to Sentinel-1 Interfeometric Wide Swath Mode (IW).' ' (IW; 250 km).') default=250*km, help='Track length in [m]. Defaults to 200 km.') default=29.1, help='Near range incident angle in [deg]. Defaults to 29.1 deg;' ' Sentinel IW mode (29.1 - 46.0 deg).') default=(250, 250), help='Resolution of raster in east x north [px].') default=True, help='Mask out water bodies.') default=AtmosphericNoiseGenerator.D(), help='Add atmospheric noise model after Hansen, 2001.')
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
targets = [s.get_target() for s in self.get_scene_patches()]
for t in targets: t.store_id = self.store_id
return targets
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 = float(tmin) sc.meta.time_slave = 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
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)
logger.warning('InSAR mapping is not implemented!') return None |