|
#!/usr/bin/env python3
Dict, load
try: return Scene.load(filename) except (ImportError, UserIOWarning): pass try: scene = Scene() scene.import_data(filename) return scene except ImportError: pass raise ImportError('Could not read file %s' % filename)
else: raise TypeError('value must be of type numpy.ndarray')
"""Config object holding :class:`kite.scene.Scene` configuration """ default=0., help='Scene latitude of lower left corner') default=0., help='Scene longitude of lower left corner') default=25., help='Scene pixel spacing in north, give [m] or [deg]') default=25., help='Scene pixel spacing in east, give [m] or [deg]') choices=('degree', 'meter'), default='meter', help='Unit of pixel space')
'dE': 'dLon', 'dN': 'dLat' }
kwargs[new] = kwargs.pop(old) kwargs['spacing'] = 'degree' self.old_import = True
""" Frame holding geographical references for :class:`kite.scene.Scene`
The pixel spacing is given by ``dE`` and ``dN`` which can meters or degree. """
elif self.config != self._scene.config.frame: self.config = self._scene.config.frame else: return
self._log.warning('Importing an old kite format...') self._log.warning('Please check your pixel spacing - dE, dN!')
self.llNutm, self.utm_zone, self.utm_zone_letter) = utm.from_latlon(self.llLat, self.llLon)
def llLat(self):
def llLat(self, llLat): self.config.llLat = llLat self.updateExtent()
def llLon(self):
def llLon(self, llLon): self.config.llLon = llLon self.updateExtent()
def dN(self):
def dN(self, dN): self.config.dN = dN self.updateExtent()
def dE(self):
def dE(self, dE): self.config.dE = dE self.updateExtent()
def dEmeter(self):
self.llLat, self.llLon, self.llLat, self.llLon + self.dE * self.cols)
def dNmeter(self): self.llLat, self.llLon, self.llLat + self.dN * self.rows, self.llLon)
def dEdegree(self): if self.isDegree(): return self.dE
lat, lon = ne_to_latlon( self.llLat, self.llLon, 0., self.dE * self.cols) distLon = lon - self.llLon return distLon / self.cols
def dNdegree(self): if self.isDegree(): return self.dE
lat, lon = ne_to_latlon( self.llLat, self.llLon, self.dN * self.rows, 0.) distLat = lat - self.llLat return distLat / self.rows
def spacing(self): return self.config.spacing
def spacing(self, unit): self.config.spacing = unit
def E(self):
def Emeter(self): return num.arange(self.cols) * self.dEmeter
def N(self):
def lengthE(self): return self.cols * self.dE
def lengthN(self): return self.rows * self.dN
def Nmeter(self): return num.arange(self.rows) * self.dNmeter
def gridE(self): """ Grid holding local east coordinates of all pixels in ``NxM`` matrix of :attr:`~kite.Scene.displacement`.
:type: :class:`numpy.ndarray`, size ``NxM`` """ self.rows, axis=0)
def gridN(self): """ Grid holding local north coordinates of all pixels in ``NxM`` matrix of :attr:`~kite.Scene.displacement`.
:type: :class:`numpy.ndarray`, size ``NxM`` """ self.cols, axis=1)
raise ValueError('Frame is defined in meter! ' 'Use gridE and gridN for meter grids')
self.llLat, self.llLon, self.llLat + self.gridN.data.ravel(), self.llLon + self.gridE.data.ravel())
gridE.reshape(self.gridE.shape), valid_data, fill_value=num.nan) gridN.reshape(self.gridN.shape), valid_data, fill_value=num.nan)
def gridEmeter(self):
def gridNmeter(self):
def coordinates(self): """ Local east and north coordinates of all pixels in ``Nx2`` matrix.
:type: :class:`numpy.ndarray`, size ``Nx2`` """ self.rows, axis=0).flatten() self.cols, axis=1).flatten()
else:
def coordinatesMeter(self): """ Local east and north coordinates [m] of all pixels in ``NxM`` matrix.
:type: :class:`numpy.ndarray`, size ``NxM`` """ coords = num.empty((self.rows*self.cols, 2)) coords[:, 0] = num.repeat(self.Emeter[num.newaxis, :], self.rows, axis=0).flatten() coords[:, 1] = num.repeat(self.Nmeter[:, num.newaxis], self.cols, axis=1).flatten() return coords
""" Local map coordinates in east and north to matrix row and column
:param E: Easting in local coordinates :type E: float :param N: Northing in local coordinates :type N: float :returns: Row and column :rtype: tuple (int), (row, column) """
def shape(self):
def npixel(self): return self.cols * self.rows
self.llLon == other.llLon and\ self.dE == other.dE and\ self.dN == other.dN and\ self.rows == other.rows and\ self.cols == other.cols
""" Meta configuration for ``Scene``. """ default='Unnamed Scene', help='Scene title') default='None', help='Scene identification') default='Undefined Mission', help='Satellite mission name') optional=True, help='Wavelength in [m]') choices=['Ascending', 'Descending', 'Undefined'], default='Undefined', help='Orbital direction, ascending/descending') default=1481116161.930574, help='Timestamp for master acquisition') default=1482239325.482, help='Timestamp for slave acquisition') default={}, help='Extra header information') optional=True)
'orbit_direction': 'orbital_node' }
kwargs[new] = kwargs.pop(old, None) self.old_import = True
def time_separation(self): """ :getter: Absolute time difference between ``time_master`` and ``time_slave`` :type: timedelta """ return dt.fromtimestamp(self.time_slave) -\ dt.fromtimestamp(self.time_master)
""" Configuration object, gathering ``kite.Scene`` and sub-objects configuration. """ default=Meta.D(), help='Scene metainformation') default=FrameConfig.D(), help='Frame/reference configuration') default=QuadtreeConfig.D(), help='Quadtree parameters') default=CovarianceConfig.D(), help='Covariance parameters') default=APSConfig.D(), help='Empirical APS correction') default=GACOSConfig.D(), help='GACOS APS correction') default=PolygonMaskConfig.D(), help='Displacement mask polygon') default=DerampConfig.D(), help='Displacement deramp config')
def old_import(self):
"""Decorator for dynamic classmethod / instancemethod declaration """ if isinstance(args[0], Scene): return func(*args, **kwargs) else: return func(Scene(), *args, **kwargs)
frame_config.__setattr__(fattr, coord)
def displacement(self): """Surface displacement in meter on a regular grid.
:setter: Set the unwrapped InSAR displacement. :getter: Return the displacement matrix. :type: :class:`numpy.ndarray`, ``NxM`` """ return self._displacement
def displacement(self, value): _setDataNumpy(self, '_displacement', value) self.rows, self.cols = self._displacement.shape self.evChanged.notify()
def displacement_px_var(self): """ Variance of the surface displacement per pixel. Same dimension as displacement.
:setter: Set standard deviation of of the displacement. :getter: Return the standard deviation matrix. :type: :class:`numpy.ndarray`, ``NxM`` """
def displacement_px_var(self, value): self._displacement_px_var = value
def displacement_mask(self): """ Displacement :attr:`numpy.nan` mask
:type: :class:`numpy.ndarray`, dtype :class:`numpy.bool` """
def shape(self):
def phi(self): """ Horizontal angle towards satellite :abbr:`line of sight (LOS)` in radians counter-clockwise from East.
.. important ::
Kite's convention is:
* :math:`0` is **East** * :math:`\\frac{\\pi}{2}` is **North**!
:setter: Set the phi matrix for scene's displacement, can be ``int`` for static look vector. :type: :class:`numpy.ndarray`, size same as :attr:`~kite.Scene.displacement` or int """
def phi(self, value): self._phi = value else:
def theta(self): """ Theta is the look vector elevation angle towards satellite from the horizon in radians. Matrix of theta towards satellite's :abbr:`line of sight (LOS)`.
.. important ::
Kite convention!
* :math:`-\\frac{\\pi}{2}` is **Down** * :math:`\\frac{\\pi}{2}` is **Up**
:setter: Set the theta matrix for scene's displacement, can be ``int`` for static look vector. :type: :class:`numpy.ndarray`, size same as :attr:`~kite.Scene.displacement` or int """
def theta(self, value): self._theta = value else:
def thetaDeg(self): """ LOS elevation angle in degree, ``NxM`` matrix like :class:`kite.Scene.theta`
:type: :class:`numpy.ndarray` """ return num.rad2deg(self.theta)
def phiDeg(self): """ LOS horizontal orientation angle in degree, counter-clockwise from East,``NxM`` matrix like :class:`kite.Scene.phi`
:type: :class:`numpy.ndarray` """ return num.rad2deg(self.phi)
def los_rotation_factors(self): """ Trigonometric factors to rotate displacement matrices towards LOS
Rotation is as follows:
.. displacement_los =\ (los_rotation_factors[:, :, 0] * -down + los_rotation_factors[:, :, 1] * east + los_rotation_factors[:, :, 2] * north)
:returns: Factors for rotation :rtype: :class:`numpy.ndarray`, ``NxMx3`` :raises: AttributeError """ if (self.theta.size != self.phi.size): raise AttributeError('LOS angles inconsistent with provided' ' coordinate shape.') if self._los_factors is None: self._los_factors = num.empty((self.theta.shape[0], self.theta.shape[1], 3)) self._los_factors[:, :, 0] = num.sin(self.theta) self._los_factors[:, :, 1] = num.cos(self.theta)\ * num.cos(self.phi) self._los_factors[:, :, 2] = num.cos(self.theta)\ * num.sin(self.phi) return self._los_factors
# region = llLon, urLon, llLat, urLon
raise AssertionError( 'Region is outside of SRTMGL3 topo dataset')
raise AssertionError('Cannot get SRTMGL3 topo dataset')
elif interpolation == 'bivariate': interp = interpolate.RectBivariateSpline( tile.y(), tile.x(), tile.data) elevation = interp(lats, lons, grid=False)
ret = copy.deepcopy(self) ret.displacement *= -1 return ret
ret = copy.deepcopy(self) else:
raise ValueError('Scene frames do not align!')
if ret.meta.time_master < other.meta.time_master \ else other.meta.time_master
if ret.meta.time_slave > other.meta.time_slave \ else other.meta.time_slave
return self.__add__(-other)
return self.__add__(-scene, copy_obj=False)
"""Scene of unwrapped InSAR ground displacements measurements
:param config: Configuration object :type config: :class:`~kite.scene.SceneConfig`, optional
Optional parameters
:param displacement: Displacement in [m] :type displacement: :class:`numpy.ndarray`, NxM, optional :param theta: Theta look angle, see :attr:`BaseScene.theta` :type theta: :class:`numpy.ndarray`, NxM, optional :param phi: Phi look angle, see :attr:`BaseScene.phi` :type phi: :class:`numpy.ndarray`, NxM, optional
:param llLat: Lower left latitude in [deg] :type llLat: float, optional :param llLon: Lower left longitude in [deg] :type llLon: float, optional :param dLat: Pixel spacing in latitude [deg or m] :type dLat: float, optional :param dLon: Pixel spacing in longitude [deg or m] :type dLon: float, optional """
# wiring special methods
self.gacos: None, self.aps: None, self.polygon_mask: None, self.deramp: None }
def displacement(self):
plugin.__class__.__name__, time.time() - t)
def displacement(self, value):
def quadtree(self): """ Instantiates the scene's quadtree.
:type: :class:`kite.quadtree.Quadtree` """
def covariance(self): """ Instantiates the scene's covariance attribute.
:type: :class:`kite.covariance.Covariance` """
def plot(self): """ Shows a simple plot of the scene's displacement """ self._log.debug('Creating kite.ScenePlot instance') from kite.plot2d import ScenePlot return ScenePlot(self)
'processing states changed: %s', plugin.__class__.__name__)
sha = hashlib.sha1() for plugin, state in self.processing_states.items(): sha.update(plugin.get_state_hash().encode())
sha.update(self.covariance.get_state_hash().encode()) sha.update(self.quadtree.get_state_hash().encode()) return sha.hexdigest()
""" Start the spool user interface :class:`~kite.spool.Spool` to inspect the scene. """ if self.displacement is None: raise SceneError('Can not display an empty scene.')
from kite.spool import spool spool(scene=self)
except Exception as e: self._log.exception(e) raise ImportError('Something went wrong during import - ' 'see Exception!')
""" Save kite scene to kite file structure
Saves the current scene meta information and UTM frame to a YAML (``.yml``) file. Numerical data (:attr:`~kite.Scene.displacement`, :attr:`~kite.Scene.theta` and :attr:`~kite.Scene.phi`) are saved as binary files from :class:`numpy.ndarray`.
:param filename: Filenames to save scene to, defaults to ' :attr:`~kite.Scene.meta.scene_id` ``_`` :attr:`~kite.Scene.meta.scene_view` :type filename: str, optional """ self.meta.scene_view)
*[getattr(self, arr) for arr in components])
header='kite.Scene YAML Config')
def load(cls, filename): """ Load a kite scene from file ``filename.[npz,yml]`` structure.
:param filename: Filenames the scene data is saved under :type filename: str :returns: Scene object from data resources :rtype: :class:`~kite.Scene` """
except IOError: raise UserIOWarning('Could not load data from %s.npz' % basename)
except IOError: raise UserIOWarning('Could not load %s.yml' % basename)
displacement=displacement, theta=theta, phi=phi, config=config)
self._log.debug('Loading config from %s', filename) self.config = load(filename=filename) self.meta = self.config.meta
self.evConfigChanged.notify()
def _import_data(self, path, **kwargs): """ Import displacement data from foreign file format.
:param path: Filename of resource to import :type path: str :param kwargs: keyword arguments passed to import function :type kwargs: dict :returns: Scene from path :rtype: :class:`~kite.Scene` :raises: TypeError """ scene = self if not op.isfile(path) and not op.isdir(path): raise ImportError('File %s does not exist!' % path)
data = None
for mod_name in scene_io.__all__: cls = getattr( __import__('kite.scene_io', fromlist=mod_name), mod_name) module = cls() if module.validate(path, **kwargs): scene._log.debug('Importing %s using %s module' % (path, mod_name)) data = module.read(path, **kwargs) break if data is None: raise ImportError('Could not recognize format for %s' % path)
scene.meta.filename = op.abspath(path) return scene._import_from_dict(scene, data)
'\nSupported import for unwrapped InSAR data are:\n\n{}\n'\ .format('\n'.join(_class_list)) __import__('kite.scene_io', fromlist=mod_name), mod_name)
.format(name=mod_name, doc=cls.__doc__)
def _import_from_dict(scene, data): for sk in ['theta', 'phi', 'displacement']: setattr(scene, sk, data[sk])
for fk, fv in data['frame'].items(): setattr(scene.frame, fk, fv)
for mk, mv in data['meta'].items(): if mv is not None: setattr(scene.meta, mk, mv) scene.meta.extra.update(data['extra']) scene.frame.updateExtent()
scene._testImport() return scene
def __str__(self): return self.config.__str__()
""" Decompose line-of-sight (LOS) angles derived from :attr:`~kite.Scene.displacement` to unit vector. """
def unitE(self): """ Unit vector east component, ``NxM`` matrix like :attr:`~kite.Scene.displacement` :type: :class:`numpy.ndarray` """ return self._scene.los_rotation_factors[:, :, 1]
def unitN(self): """ Unit vector north component, ``NxM`` matrix like :attr:`~kite.Scene.displacement` :type: :class:`numpy.ndarray` """ return self._scene.los_rotation_factors[:, :, 2]
def unitU(self): """ Unit vector vertical (up) component, ``NxM`` matrix like :attr:`~kite.Scene.displacement` :type: :class:`numpy.ndarray` """ return self._scene.los_rotation_factors[:, :, 0]
"""Test scenes for synthetic displacement """
scene = cls() scene.meta.scene_title = 'Synthetic Displacement | Gaussian' scene = cls._prepareSceneTest(scene, nx, ny)
scene.displacement = scene._gaussAnomaly(scene.frame.E, scene.frame.N, **kwargs) if noise is not None: cls.addNoise(noise) return scene
scene = cls() scene.meta.title = 'Synthetic Displacement | Uniform Random' scene = cls._prepareSceneTest(scene, nx, ny)
rand_state = num.random.RandomState(seed=kwargs.pop('seed', None)) scene.displacement = (rand_state.rand(nx, ny)-.5)*2
return scene
noise=.5, **kwargs): scene = cls() scene.meta.title = 'Synthetic Displacement | Sine' scene = cls._prepareSceneTest(scene, nx, ny)
E, N = num.meshgrid(scene.frame.E, scene.frame.N) displ = num.zeros_like(E)
kE = num.random.rand(3) * kE kN = num.random.rand(3) * kN
for ke in kE: phase = num.random.randn(1)[0] displ += num.sin(ke * E + phase) for kn in kN: phase = num.random.randn(1)[0] displ += num.sin(kn * N + phase) displ -= num.mean(displ)
scene.displacement = displ * amplitude if noise is not None: scene.addNoise(noise) return scene
beta=[5./3, 8./3, 2./3], regime=[.15, .99, 1.], amplitude=1.): scene = cls() scene.meta.title =\ 'Synthetic Displacement | Fractal Noise (Hanssen, 2001)' scene = cls._prepareSceneTest(scene, nE, nN) if (nE+nN) % 2 != 0: raise ArithmeticError('Dimensions of synthetic scene must ' 'both be even!')
dE, dN = (scene.frame.dE, 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)
regime = num.array(regime) k0 = 0. k1 = regime[0] * k_rad.max() k2 = regime[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(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 *= amplitude * num.sqrt(amp) disp = num.abs(num.fft.ifft2(spec)) disp -= num.mean(disp)
scene.displacement = disp return scene
rand = num.random.RandomState(seed) noise = rand.randn(*self.displacement.shape) * noise_amplitude self.displacement += noise
scene.frame.llLat = 0. scene.frame.llLon = 0. scene.frame.dLat = 5e-4 scene.frame.dLon = 5e-4 # scene.frame.E = num.arange(nE) * 50. # scene.frame.N = num.arange(nN) * 50. scene.theta = num.repeat( num.linspace(0.8, 0.85, nE), nN).reshape((nE, nN)) scene.phi = num.rot90(scene.theta) scene.displacement = num.zeros((nE, nN)) return scene
amplitude=3., x0=None, y0=None): if x0 is None: x0 = x.min() + abs(x.max()-x.min())/2 if y0 is None: y0 = y.min() + abs(y.max()-y.min())/2 X, Y = num.meshgrid(x, y)
gauss_anomaly = amplitude * \ num.exp(-(((X-x0)**2/2*sigma_x**2)+(Y-y0)**2/2*sigma_y**2))
return gauss_anomaly
if __name__ == '__main__': testScene = TestScene.createGauss() |