Sentinel 1 interface

kafka/__init__.py

@@ -1,5 +1,7 @@
-__all__ = ['inference','linear_kf.LinearKalman','input_output']
+__all__ = ['inference','linear_kf.LinearKalman','input_output',
+           'observation_operators']
 # deprecated to keep older scripts who import this from breaking
 from .inference import *
 from .input_output import *
+from .observation_operators import *
 from linear_kf import LinearKalman

kafka/inference/kf_tools.py

@@ -10,6 +10,11 @@
 
 from utils import block_diag
 
+class NoHessianMethod(Exception):
+    """An exception triggered when the forward model isn't able to provide an
+    estimation of the Hessian"""
+    def __init__(self, message):
+        self.message = message
 
 def band_selecta(band):
     if band == 0:
@@ -33,6 +38,10 @@ def hessian_correction(gp, x0, R_mat, innovation, mask, state_mask, band,
                        nparams):
     """Calculates higher order Hessian correction for the likelihood term.
     Needs the GP, the Observational uncertainty, the mask...."""
+    if not hasattr(gp, "hessian"):
+        # The observation operator does not provide a Hessian method. We just
+        # return 0, meaning no Hessian correction.
+        return 0.
     C_obs_inv = R_mat.diagonal()[state_mask.flatten()]
     mask = mask[state_mask].flatten()
     little_hess = []
@@ -59,12 +68,13 @@ def tip_prior():
     Returns
     -------
     The mean prior vector, covariance and inverse covariance matrices."""
-    sigma = np.array([0.12, 0.7, 0.0959, 0.15, 1.5, 0.2, 0.5]) # broadly TLAI 0->7 for 1sigma
+    # broadly TLAI 0->7 for 1sigma
+    sigma = np.array([0.12, 0.7, 0.0959, 0.15, 1.5, 0.2, 0.5])
     x0 = np.array([0.17, 1.0, 0.1, 0.7, 2.0, 0.18, np.exp(-0.5*1.5)])
     # The individual covariance matrix
-    little_p = np.diag ( sigma**2).astype(np.float32)
-    little_p[5,2] = 0.8862*0.0959*0.2
-    little_p[2,5] = 0.8862*0.0959*0.2
+    little_p = np.diag(sigma**2).astype(np.float32)
+    little_p[5, 2] = 0.8862*0.0959*0.2
+    little_p[2, 5] = 0.8862*0.0959*0.2
 
     inv_p = np.linalg.inv(little_p)
     return x0, little_p, inv_p

kafka/input_output/Sentinel1_Observations.py

@@ -0,0 +1,193 @@
+#!/usr/bin/env python
+"""
+
+"""
+import datetime
+import glob
+import os
+from collections import namedtuple
+
+import gdal
+
+import numpy as np
+import pdb
+import osr
+
+import scipy.sparse as sp
+
+WRONG_VALUE = -999.0  # TODO tentative missing value
+
+SARdata = namedtuple('SARdata',
+                     'observations uncertainty mask metadata emulator')
+
+
+def reproject_image(source_img, target_img, dstSRSs=None):
+    """Reprojects/Warps an image to fit exactly another image.
+    Additionally, you can set the destination SRS if you want
+    to or if it isn't defined in the source image."""
+    g = gdal.Open(target_img)
+    geo_t = g.GetGeoTransform()
+    x_size, y_size = g.RasterXSize, g.RasterYSize
+    xmin = min(geo_t[0], geo_t[0] + x_size * geo_t[1])
+    xmax = max(geo_t[0], geo_t[0] + x_size * geo_t[1])
+    ymin = min(geo_t[3], geo_t[3] + y_size * geo_t[5])
+    ymax = max(geo_t[3], geo_t[3] + y_size * geo_t[5])
+    xRes, yRes = abs(geo_t[1]), abs(geo_t[5])
+    if dstSRSs is None:
+        dstSRS = osr.SpatialReference()
+        raster_wkt = g.GetProjection()
+        dstSRS.ImportFromWkt(raster_wkt)
+    else:
+        dstSRS = dstSRSs
+    g = gdal.Warp('', source_img, format='MEM',
+                  outputBounds=[xmin, ymin, xmax, ymax], xRes=xRes, yRes=yRes,
+                  dstSRS=dstSRS)
+    if g is None:
+        raise ValueError("Something failed with GDAL!")
+    return g
+
+
+class S1Observations(object):
+    """
+    """
+
+    def __init__(self, data_folder, state_mask,
+                 emulators={'VV': 'SOmething', 'VH': 'Other'}):
+
+        """
+
+        """
+        # 1. Find the files
+        files = glob.glob(os.path.join(data_folder, '*.nc'))
+        files.sort()
+        self.state_mask = state_mask
+        self.dates = []
+        self.date_data = {}
+        for fich in files:
+            fname = os.path.basename(fich)
+            splitter = fname.split('_')
+            this_date = datetime.datetime.strptime(splitter[5],
+                                                   '%Y%m%dT%H%M%S')
+            self.dates.append(this_date)
+            self.date_data[this_date] = fich
+        # 2. Store the emulator(s)
+        self.emulators = emulators
+
+    def _read_backscatter(self, obs_ptr):
+        """
+        Read backscatter from NetCDF4 File
+        Should return a 2D array that should overlap with the "state grid"
+
+        Input
+        ------
+        fname (string, filename of the NetCDF4 file, path included)
+        polarisation (string, used polarisation)
+
+        Output
+        ------
+        backscatter values stored within NetCDF4 file for given polarisation
+
+        """
+        backscatter = obs_ptr.ReadAsArray()
+        return backscatter
+
+    def _calculate_uncertainty(self, backscatter):
+        """
+        Calculation of the uncertainty of Sentinel-1 input data
+
+        Radiometric uncertainty of Sentinel-1 Sensors are within 1 and 0.5 dB
+
+        Calculate Equivalent Number of Looks (ENL) of input dataset leads to
+        uncertainty of scene caused by speckle-filtering/multi-looking
+
+        Input
+        ------
+        backscatter (backscatter values)
+
+        Output
+        ------
+        unc (uncertainty in dB)
+
+        """
+
+        # first approximation of uncertainty (1 dB)
+        unc = backscatter*0.05
+
+
+        # need to find a good way to calculate ENL
+        # self.ENL = (backscatter.mean()**2) / (backscatter.std()**2)
+
+        return unc
+
+    def _get_mask(self, backscatter):
+        """
+        Mask for selection of pixels
+
+        Get a True/False array with the selected/unselected pixels
+
+
+        Input
+        ------
+        this_file ?????
+
+        Output
+        ------
+
+        """
+
+        mask = np.ones_like(backscatter, dtype=np.bool)
+        mask[backscatter == WRONG_VALUE] = False
+        return mask
+
+    def get_band_data(self, timestep, band):
+        """
+        get all relevant S1 data information for one timestep to get processing
+        done
+
+
+        Input
+        ------
+        timestep
+        band
+
+        Output
+        ------
+        sardata (namedtuple with information on observations, uncertainty,
+                mask, metadata, emulator/used model)
+
+        """
+
+        if band == 0:
+            polarisation = 'VV'
+        elif band == 1:
+            polarisation = 'VH'
+        this_file = self.date_data[timestep]
+        fname = 'NETCDF:"{:s}":sigma0_{:s}'.format(this_file, polarisation)
+        obs_ptr = reproject_image(fname, self.state_mask)
+        observations = self._read_backscatter(obs_ptr)
+        uncertainty = self._calculate_uncertainty(observations)
+        mask = self._get_mask(observations)
+        R_mat = np.zeros_like(observations)
+        R_mat = uncertainty
+        R_mat[np.logical_not(mask)] = 0.
+        N = mask.ravel().shape[0]
+        R_mat_sp = sp.lil_matrix((N, N))
+        R_mat_sp.setdiag(1./(R_mat.ravel())**2)
+        R_mat_sp = R_mat_sp.tocsr()
+
+        emulator = self.emulators[polarisation]
+        # TODO read in angle of incidence from netcdf file
+        # metadata['incidence_angle_deg'] =
+        fname = 'NETCDF:"{:s}":{:s}'.format(this_file, "theta")
+        obs_ptr = reproject_image(fname, self.state_mask)
+        relativeorbit = obs_ptr.GetMetadata()['NC_GLOBAL#relativeorbit']
+        orbitdirection = obs_ptr.GetMetadata()['NC_GLOBAL#orbitdirection']
+        satellite = obs_ptr.GetMetadata()['NC_GLOBAL#satellite']
+        metadata = {'incidence_angle': obs_ptr.ReadAsArray(), 'relativeorbit': relativeorbit, 'orbitdirection': orbitdirection, 'satellite': satellite}
+        sardata = SARdata(observations, R_mat_sp, mask, metadata, emulator)
+        return sardata
+
+
+if __name__ == "__main__":
+    data_folder = "/media/ucfajlg/WERKLY/jose/new/"
+    sentinel1 = S1Observations(data_folder)

kafka/input_output/__init__.py

@@ -1,3 +1,4 @@
-__all__ = ["observations"]
+__all__ = ["observations", "Sentinel1_Observations"]
 
 from .observations import *
+from .Sentinel1_Observations import *