From fedf315a76ba4cca50a352b157b06b530406660b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20G=C3=B3mez-Dans?= Date: Tue, 9 Jan 2018 17:48:58 +0000 Subject: [PATCH 01/18] Added a BHR test class The new observation class allows one to hopefully read in a single "pixel" of BHR data over time to do some quick tests. In essence, it just needs a set of dates, albedos and a gp (this bit is missing) It then needs to pack the data into a BHR structure. --- kafka/input_output/observations.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/kafka/input_output/observations.py b/kafka/input_output/observations.py index 9ad6fd4..82c1b2a 100644 --- a/kafka/input_output/observations.py +++ b/kafka/input_output/observations.py @@ -294,6 +294,25 @@ def get_band_data(self, the_date, band_no): bhr_data = BHR_data(bhr, mask, R_mat_sp, None, self.emulator) return bhr_data + +class BHRObservationsTest(object): + """A class to test BHR data "one pixel at a time". In essence, one only needs + to define a self.dates dictionary (keys are datetime objects), and a 2 element + list or array with the VIS/NIR albedo. then we need the get_band_method...""" + def __init__(self, dates, vis_albedo, nir_albedo): + assert (len(dates) == len(vis_albedo)) + assert (len(dates) == len(nir_albedo)) + self.dates = {} + for ii, the_date in enumerate(dates): + self.dates[the_date] = [vis_albedo[ii], nir_albedo[ii]] + + def get_band_data(self, the_date, band_no): + bhr = self.dates[the_date][band_no] + mask = np.array(1, dtype=np.bool) + R_mat = 1./(np.maximum(2.5e-3, bhr * 0.05))**2 + + + class KafkaOutput(object): """A very simple class to output the state.""" From 5b9490b78a8b7c6b468b2c37cdd4e37f4bdfa6cd Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 10:05:19 +0000 Subject: [PATCH 02/18] Added SAR observation operator Added SAR forward model, as well as code that crates the linear approximation to the SAR model. The SAR observations have yet to be included in the library. Obviously, untested code. --- .../sar_forward_model.py | 149 ++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 kafka/observation_operators/sar_forward_model.py diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py new file mode 100644 index 0000000..2a8d846 --- /dev/null +++ b/kafka/observation_operators/sar_forward_model.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python +""" + +""" + +import numpy as np +import pdb + +def sar_observation_operator(x, polarisation): + + """ + For the sar_observation_operator a simple Water Cloud Model (WCM) is used + We assume that the WCM is given by + + ----------------------------------------------------- + tau = exp(-2*B*V/cos(theta)) + sigma_veg = A*V**E*cos(theta)*(1-tau) + sigma_soil = 10**((C+D*SM)/10) + sigma_0 = sigma_veg + tau*sigma_soil + + A, B, C, D, E are parameters determined by fitting of the model + V is a vegetation descriptor e.g. LAI, LWAI, VWM + SM is the soil moisture [m^3/m^3] + sigma_veg is the volume component of the backscatter [m^3/m^3] + sigma_soil is the surface component of the backscatter [m^3/m^3] + tau is the two way attenuation through the canopy (unitless) + sigma_0 is the total backscatter in [m^3/m^3] + ---------------------------------------------------- + + Input + ----- + polarisation: considered polarisation as string + x: 2D array where every row is the set of parameters for one pixel + + Output + ------ + sigma_0: predicted backscatter for each individual parameter set + grad: gradient for each individual parameter set and each parameter determined by 2D input array x. The gradient should thus have the same size and shape as x + sigma_veg: predicted volume component of sigma_0 + sigma_surf: predicted surface component of sigma_0 + tau: predicted two-way attenuation through the canopy + """ + + # x 2D array where every row is the set of parameters for one pixel + x = np.atleast_2d(x) + + # conversion of incidence angle to radiant + # the incidence angle itself should probably implemented in x) + theta = np.deg2rad(23.) + theta = np.deg2rad(np.arange(1.,80.)) + + # Simpler definition of cosine of theta + mu = np.cos(theta) + + # the model parameters (A, B, C, D, E) for different polarisations + parameters = {'VV': [0.0846, 0.0615, -14.8465, 15.907, 1.], + 'VH': [0.0795, 0.1464, -14.8332, 15.907, 0.]} + + # Select model parameters + try: + A, B, C, D, E = parameters[polarisation] + except KeyError: + raise ValueError('Only VV and VH polarisations available!') + + # Calculate Model + tau = np.exp(-2 * B / mu * x[0, :]) + sigma_veg = A * (x[0, :] ** E) * mu * (1 - tau) + sigma_surf = 10 ** ((C + D * x[1, :]) / 10.) + + sigma_0 = sigma_veg + tau * sigma_surf + #pdb.set_trace() + + # Calculate Gradient (grad has same dimension as x) + grad = x*0 + n_elems = x.shape[1] + for i in range(n_elems): + tau_value = np.exp(-2 * B / mu * x[0, i]) + grad[0, i] = A * E * mu * (x[0, i] ** (E - 1)) * (1 - tau_value) + \ + 2 * A * B * (x[0, i] ** E) * tau_value - (\ + (2 ** (1/10. * (C + D * x[1, i]) + 1)) * \ + (5 ** (1/10. * (C + D * x[1, i])) * B * tau_value) \ + ) / mu + grad[1, i] = D * np.log(10) * tau_value * 10 ** (1/10. * (C + D * x[1, i]) - 1) + + + # returned values are linear scaled not dB!!! + # return sigma_0, grad, sigma_veg, sigma_surf, tau + return sigma_0, grad + + +def create_nonlinear_observation_operator(n_params, emulator, metadata, + mask, state_mask, x_forecast, band): + """Using an emulator of the nonlinear model around `x_forecast`. + This case is quite special, as I'm focusing on a BHR SAIL + version (or the JRC TIP), which have spectral parameters + (e.g. leaf single scattering albedo in two bands, etc.). This + is achieved by using the `state_mapper` to select which bits + of the state vector (and model Jacobian) are used.""" + LOG.info("Creating the ObsOp for band %d" % band) + n_times = x_forecast.shape[0] / n_params + good_obs = mask.sum() + H_matrix = sp.lil_matrix((n_times, n_params * n_times), + dtype=np.float32) + H0 = np.zeros(n_times, dtype=np.float32) + + + # So the model has spectral components. + if band == 0: + # VV + polarisation = "VV" + elif band == 1: + # VH + polarisation = "VH" + + + # This loop can be JIT'ed + x0 = np.zeros((n_times, n_params)) + for i, m in enumerate(mask[state_mask].flatten()): + if m: + x0[i, :] = x_forecast[(n_params * i): (n_params*(i+1))] + LOG.info("Running SAR forward model") + # Calls the run_emulator method that only does different vectors + # It might be here that we do some sort of clustering + + H0_, dH = sar_observation_operator(x0[mask[state_mask]], polarisation) + + LOG.info("Storing emulators in H matrix") + # This loop can be JIT'ed too + n = 0 + for i, m in enumerate(mask[state_mask].flatten()): + if m: + H_matrix[i, (n_params * i): (n_params*(i+1))] = dH[n] + H0[i] = H0_[n] + n += 1 + LOG.info("\tDone!") + + return (H0, H_matrix.tocsr()) + + # # Calculate Gradient without conversion of sigma_soil from dB to linear + # grad = x*0 + # n_elems = x.shape[0] + # for i in range(n_elems): + # tau = np.exp(-2 * B / mu * x[i, 0]) + # grad[i, 0] = 2 * A * B * x[i, 0] * tau - \ + # A * mu * tau + \ + # A * mu - \ + # 2 * B * C * tau / mu - \ + # 2 * B * D * x[i, 1] * tau / mu + # grad[i, 1] = D * tau From 920502f7a76cdefb134dbae3e3f2c9cf9dbb2352 Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 10:11:34 +0000 Subject: [PATCH 03/18] Added documentation --- .../sar_forward_model.py | 33 +++++++++++++++---- 1 file changed, 26 insertions(+), 7 deletions(-) diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index 2a8d846..0f5df5b 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -88,14 +88,33 @@ def sar_observation_operator(x, polarisation): return sigma_0, grad -def create_nonlinear_observation_operator(n_params, emulator, metadata, +def create_sar_observation_operator(n_params, forward_model, metadata, mask, state_mask, x_forecast, band): - """Using an emulator of the nonlinear model around `x_forecast`. - This case is quite special, as I'm focusing on a BHR SAIL - version (or the JRC TIP), which have spectral parameters - (e.g. leaf single scattering albedo in two bands, etc.). This - is achieved by using the `state_mapper` to select which bits - of the state vector (and model Jacobian) are used.""" + """Creates the SAR observation operator using the Water Cloud SAR forward + model (defined above). + + Parameters + ----------- + n_params: int + Number of parameters in the state vector per pixel + emulator: function + The function to call the forward model. Defined above + metadata: list + Not used + mask: array + A 2D mask with the observational valid pixels + state_mask: array + A 2D mask with the pixels that will be used for state inference + x_forecast: array + The state vector around which the linearisation will be done. Order is + parameters per pixel (e.g. LAI1, SM1, LAI2, SM2, ..., LAIN, SMN) + band: array + The band number. Assumed 0 is VV and 1 is VH. + + Returns + -------- + H0, dH + """ LOG.info("Creating the ObsOp for band %d" % band) n_times = x_forecast.shape[0] / n_params good_obs = mask.sum() From 930d22ab9e1768eb06feaf344e1c748be634b7bf Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 10:13:05 +0000 Subject: [PATCH 04/18] Added sparse imports --- kafka/observation_operators/sar_forward_model.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index 0f5df5b..79adee6 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -4,7 +4,9 @@ """ import numpy as np -import pdb + +import scipy.sparse as sp + def sar_observation_operator(x, polarisation): From 0f69bc8bd3818ab4d1e65b10b0ba617d5fd3536c Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 10:21:50 +0000 Subject: [PATCH 05/18] Fixed forward model definition --- kafka/observation_operators/sar_forward_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index 79adee6..47b6346 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -99,7 +99,7 @@ def create_sar_observation_operator(n_params, forward_model, metadata, ----------- n_params: int Number of parameters in the state vector per pixel - emulator: function + forward_model: function The function to call the forward model. Defined above metadata: list Not used @@ -143,7 +143,7 @@ def create_sar_observation_operator(n_params, forward_model, metadata, # Calls the run_emulator method that only does different vectors # It might be here that we do some sort of clustering - H0_, dH = sar_observation_operator(x0[mask[state_mask]], polarisation) + H0_, dH = forward_model(x0[mask[state_mask]], polarisation) LOG.info("Storing emulators in H matrix") # This loop can be JIT'ed too From 8c4fedf4151aacbe7bb09a34c469720da7a065bb Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 11:53:42 +0000 Subject: [PATCH 06/18] SAR observations class added Currently using NC, but might need changing to GeoTIFF. Might also need reprojecting/resampling to fit state grid etc. Fill value is currently set to 0, which isn't great as values in file span from -ve to +ve. Also data in dB probably needs to be converted to linear scale to work with RT model. --- kafka/input_output/Sentinel1_Observations.py | 138 +++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 kafka/input_output/Sentinel1_Observations.py diff --git a/kafka/input_output/Sentinel1_Observations.py b/kafka/input_output/Sentinel1_Observations.py new file mode 100644 index 0000000..9de4a7e --- /dev/null +++ b/kafka/input_output/Sentinel1_Observations.py @@ -0,0 +1,138 @@ +#!/usr/bin/env python +""" + +""" +import datetime +import glob +import os +from collections import namedtuple + +import gdal + +import numpy as np + +WRONG_VALUE = -999.0 # TODO tentative missing value + +SARdata = namedtuple('SARdata', + 'observations uncertainty mask metadata emulator') + + +class S1Observations(object): + """ + """ + + def __init__(self, data_folder, + emulators={'VV': 'SOmething', 'VH': 'Other'}): + + """ + + """ + # 1. Find the files + files = glob.glob(os.path.join(data_folder, '*.nc')) + files.sort() + self.observation_dates = {} + 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.observation_dates[this_date] = fich + + # 2. Store the emulator(s) + self.emulators = emulators + + def _read_backscatter(self, fname, polarisation): + """ + 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 + + """ + fpath = 'NETCDF:"{:s}":sigma0_{:s}'.format(fname, polarisation) + g = gdal.Open(fpath) + backscatter = g.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) + self.unc = 1. + + # need to find a good way to calculate ENL + # self.ENL = (backscatter.mean()**2) / (backscatter.std()**2) + + return self.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.where(backscatter == WRONG_VALUE) + 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.observation_dates[timestep] + observations = self._read_backscatter(this_file, polarisation) + uncertainty = self._calculate_uncertainty(observations) + mask = self._get_mask(observations) + emulator = self.emulators[polarisation] + sardata = SARdata(observations, uncertainty, mask, None, emulator) + return sardata From 613a4873b9b452d5638f69fad96963b1e31eb92c Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 18:32:03 +0000 Subject: [PATCH 07/18] Sentinel1 reader class ready The class seems to work, but untested within KaFKA --- kafka/input_output/Sentinel1_Observations.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/kafka/input_output/Sentinel1_Observations.py b/kafka/input_output/Sentinel1_Observations.py index 9de4a7e..eb2b3d4 100644 --- a/kafka/input_output/Sentinel1_Observations.py +++ b/kafka/input_output/Sentinel1_Observations.py @@ -30,13 +30,13 @@ def __init__(self, data_folder, # 1. Find the files files = glob.glob(os.path.join(data_folder, '*.nc')) files.sort() - self.observation_dates = {} + self.dates = {} 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.observation_dates[this_date] = fich + self.dates[this_date] = fich # 2. Store the emulator(s) self.emulators = emulators @@ -129,10 +129,15 @@ def get_band_data(self, timestep, band): polarisation = 'VV' elif band == 1: polarisation = 'VH' - this_file = self.observation_dates[timestep] + this_file = self.dates[timestep] observations = self._read_backscatter(this_file, polarisation) uncertainty = self._calculate_uncertainty(observations) mask = self._get_mask(observations) emulator = self.emulators[polarisation] sardata = SARdata(observations, uncertainty, mask, None, emulator) return sardata + + +if __name__ == "__main__": + data_folder = "/media/ucfajlg/WERKLY/jose/new/" + sentinel1 = S1Observations(data_folder) From 28e1e29bd3286cbfc926e9a5b97f3b90235e7434 Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Tue, 16 Jan 2018 18:32:49 +0000 Subject: [PATCH 08/18] Minor code cleanups --- .../sar_forward_model.py | 22 +++++++++---------- 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index 47b6346..93110e1 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -2,11 +2,13 @@ """ """ +import logging import numpy as np import scipy.sparse as sp +LOG = logging.getLogger(__name__) def sar_observation_operator(x, polarisation): @@ -49,7 +51,7 @@ def sar_observation_operator(x, polarisation): # conversion of incidence angle to radiant # the incidence angle itself should probably implemented in x) theta = np.deg2rad(23.) - theta = np.deg2rad(np.arange(1.,80.)) + theta = np.deg2rad(np.arange(1., 80.)) # Simpler definition of cosine of theta mu = np.cos(theta) @@ -78,12 +80,12 @@ def sar_observation_operator(x, polarisation): for i in range(n_elems): tau_value = np.exp(-2 * B / mu * x[0, i]) grad[0, i] = A * E * mu * (x[0, i] ** (E - 1)) * (1 - tau_value) + \ - 2 * A * B * (x[0, i] ** E) * tau_value - (\ - (2 ** (1/10. * (C + D * x[1, i]) + 1)) * \ - (5 ** (1/10. * (C + D * x[1, i])) * B * tau_value) \ + 2 * A * B * (x[0, i] ** E) * tau_value - ( + (2 ** (1/10. * (C + D * x[1, i]) + 1)) * + (5 ** (1/10. * (C + D * x[1, i])) * B * tau_value) ) / mu - grad[1, i] = D * np.log(10) * tau_value * 10 ** (1/10. * (C + D * x[1, i]) - 1) - + grad[1, i] = D * np.log(10) * tau_value * 10 ** ( + 1/10. * (C + D * x[1, i]) - 1) # returned values are linear scaled not dB!!! # return sigma_0, grad, sigma_veg, sigma_surf, tau @@ -91,7 +93,7 @@ def sar_observation_operator(x, polarisation): def create_sar_observation_operator(n_params, forward_model, metadata, - mask, state_mask, x_forecast, band): + mask, state_mask, x_forecast, band): """Creates the SAR observation operator using the Water Cloud SAR forward model (defined above). @@ -119,12 +121,11 @@ def create_sar_observation_operator(n_params, forward_model, metadata, """ LOG.info("Creating the ObsOp for band %d" % band) n_times = x_forecast.shape[0] / n_params - good_obs = mask.sum() + # good_obs = mask.sum() H_matrix = sp.lil_matrix((n_times, n_params * n_times), dtype=np.float32) H0 = np.zeros(n_times, dtype=np.float32) - # So the model has spectral components. if band == 0: # VV @@ -132,8 +133,6 @@ def create_sar_observation_operator(n_params, forward_model, metadata, elif band == 1: # VH polarisation = "VH" - - # This loop can be JIT'ed x0 = np.zeros((n_times, n_params)) for i, m in enumerate(mask[state_mask].flatten()): @@ -154,7 +153,6 @@ def create_sar_observation_operator(n_params, forward_model, metadata, H0[i] = H0_[n] n += 1 LOG.info("\tDone!") - return (H0, H_matrix.tocsr()) # # Calculate Gradient without conversion of sigma_soil from dB to linear From c476bf631dbc7a4ae27019388cb528b7e6cb8986 Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 12:06:23 +0000 Subject: [PATCH 09/18] Testing stuff --- kafka/input_output/observations.py | 41 +++++++++++++++++++++++++----- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/kafka/input_output/observations.py b/kafka/input_output/observations.py index 82c1b2a..3b8a4fd 100644 --- a/kafka/input_output/observations.py +++ b/kafka/input_output/observations.py @@ -299,20 +299,47 @@ class BHRObservationsTest(object): """A class to test BHR data "one pixel at a time". In essence, one only needs to define a self.dates dictionary (keys are datetime objects), and a 2 element list or array with the VIS/NIR albedo. then we need the get_band_method...""" - def __init__(self, dates, vis_albedo, nir_albedo): + def __init__(self, emulator, dates, vis_albedo, nir_albedo): assert (len(dates) == len(vis_albedo)) assert (len(dates) == len(nir_albedo)) - self.dates = {} + self.emulator = self._get_emulator(emulator) + self.dates = dates + self.data = {} for ii, the_date in enumerate(dates): - self.dates[the_date] = [vis_albedo[ii], nir_albedo[ii]] - + self.data[the_date] = [vis_albedo[ii], nir_albedo[ii]] + + + def _get_emulator(self, emulator): + if not os.path.exists(emulator): + raise IOError("The emulator {} doesn't exist!".format(emulator)) + # Assuming emulator is in an pickle file... + return cPickle.load(open(emulator, 'rb')) + + def get_band_data(self, the_date, band_no): - bhr = self.dates[the_date][band_no] - mask = np.array(1, dtype=np.bool) + bhr = self.data[the_date][band_no] + mask = np.ones((1,1), dtype=np.bool) R_mat = 1./(np.maximum(2.5e-3, bhr * 0.05))**2 + R_mat_sp = sp.lil_matrix((1,1)) + R_mat_sp.setdiag(R_mat) + R_mat_sp = R_mat_sp.tocsr() + + + bhr_data = BHR_data(bhr, mask, R_mat_sp, None, self.emulator) + return bhr_data - +class KafkaOutputBHRTest(object): + """A very simple class to output the state.""" + def __init__(self): + self.output = {} + def dump_data(self, timestep, x_analysis, P_analysis, P_analysis_inv, + state_mask): + solution = {} + for ii, param in enumerate(["w_vis", "x_vis", "a_vis", + "w_nir", "x_nir", "a_nir", "TeLAI"]): + solution[param] = x_analysis[ii::7] + self.output[timestep] = solution class KafkaOutput(object): """A very simple class to output the state.""" From 200e204e2f042f459cc2b47cf7c12489beef9852 Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 12:10:05 +0000 Subject: [PATCH 10/18] Make SAR codes discoverable to the library --- kafka/input_output/__init__.py | 3 ++- kafka/observation_operators/__init__.py | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) create mode 100644 kafka/observation_operators/__init__.py diff --git a/kafka/input_output/__init__.py b/kafka/input_output/__init__.py index 8c36463..ea76916 100644 --- a/kafka/input_output/__init__.py +++ b/kafka/input_output/__init__.py @@ -1,3 +1,4 @@ -__all__ = ["observations"] +__all__ = ["observations", "Sentinel1_Observations"] from .observations import * +from .Sentinel1_Observations import * diff --git a/kafka/observation_operators/__init__.py b/kafka/observation_operators/__init__.py new file mode 100644 index 0000000..9063770 --- /dev/null +++ b/kafka/observation_operators/__init__.py @@ -0,0 +1,3 @@ +__all__ = ["sar_forward_model"] + +from .sar_forward_model import * From b9b721c1ee7f2d3619c0ecd2a30d37583403b32f Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 13:00:16 +0000 Subject: [PATCH 11/18] Some changes to match the observations to the state mask --- kafka/input_output/Sentinel1_Observations.py | 43 ++++++++++++++++--- .../sar_forward_model.py | 1 + 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/kafka/input_output/Sentinel1_Observations.py b/kafka/input_output/Sentinel1_Observations.py index eb2b3d4..5ce0087 100644 --- a/kafka/input_output/Sentinel1_Observations.py +++ b/kafka/input_output/Sentinel1_Observations.py @@ -11,17 +11,45 @@ import numpy as np +import osr + 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, + def __init__(self, data_folder, state_mask, emulators={'VV': 'SOmething', 'VH': 'Other'}): """ @@ -30,6 +58,7 @@ def __init__(self, data_folder, # 1. Find the files files = glob.glob(os.path.join(data_folder, '*.nc')) files.sort() + self.state_mask = state_mask self.dates = {} for fich in files: fname = os.path.basename(fich) @@ -41,7 +70,7 @@ def __init__(self, data_folder, # 2. Store the emulator(s) self.emulators = emulators - def _read_backscatter(self, fname, polarisation): + def _read_backscatter(self, obs_ptr): """ Read backscatter from NetCDF4 File Should return a 2D array that should overlap with the "state grid" @@ -56,9 +85,7 @@ def _read_backscatter(self, fname, polarisation): backscatter values stored within NetCDF4 file for given polarisation """ - fpath = 'NETCDF:"{:s}":sigma0_{:s}'.format(fname, polarisation) - g = gdal.Open(fpath) - backscatter = g.ReadAsArray() + backscatter = obs_ptr.ReadAsArray() return backscatter def _calculate_uncertainty(self, backscatter): @@ -130,10 +157,14 @@ def get_band_data(self, timestep, band): elif band == 1: polarisation = 'VH' this_file = self.dates[timestep] - observations = self._read_backscatter(this_file, polarisation) + 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) emulator = self.emulators[polarisation] + # TODO read in angle of incidence from netcdf file + # metadata['incidence_angle_deg'] = sardata = SARdata(observations, uncertainty, mask, None, emulator) return sardata diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index 93110e1..a40dc1c 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -50,6 +50,7 @@ def sar_observation_operator(x, polarisation): # conversion of incidence angle to radiant # the incidence angle itself should probably implemented in x) + # TODO needs to come from the data theta = np.deg2rad(23.) theta = np.deg2rad(np.arange(1., 80.)) From dcec80bd8ec38a94f30cfd6faea4aeaaa3da0ec3 Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 13:04:26 +0000 Subject: [PATCH 12/18] Added travis file. Tests 3.6 and 2.7. Hopefully --- .travis.yml | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 .travis.yml diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..77125d2 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,25 @@ +language: python +python: + - "2.7" + - "3.6" +install: + - sudo apt-get update + - echo $TRAVIS_PYTHON_VERSION + - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then + wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh -O miniconda.sh; + else + wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; + fi + - bash miniconda.sh -b -p $HOME/miniconda + - export PATH="$HOME/miniconda/bin:$PATH" + - hash -r + - conda config --set always_yes yes --set changeps1 no + - conda update -q conda + - conda info -a + - conda create -q -n test python=$TRAVIS_PYTHON_VERSION numpy numba scipy pytest pytest-cov gdal coverage sphinx nose + - source activate test + - pip install -U codecov + - python setup.py install + +script: + - pytest From 60bacd45648bd6b8067184df08e8d5c7e439c52e Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 13:29:16 +0000 Subject: [PATCH 13/18] Added setup.py --- setup.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 setup.py diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..b95177b --- /dev/null +++ b/setup.py @@ -0,0 +1,19 @@ +#!/usr/bin/env python + +from setuptools import setup + +requirements = [ + 'pytest', + 'numpy', + 'scipy', + 'gdal', + 'BRDF_descriptors', + 'matplotlib' +] + +setup(name='KaFKA', + description='MULTIPLY KaFKA inference engine', + author='MULTIPLY Team', + packages=['kafka'], + install_requires=requirements +) From 50ad7656c0728e03b45b82538b8cd8550119118e Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 14:34:07 +0000 Subject: [PATCH 14/18] Now extracts incidence angle from files --- kafka/input_output/Sentinel1_Observations.py | 5 ++++- kafka/observation_operators/sar_forward_model.py | 8 ++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/kafka/input_output/Sentinel1_Observations.py b/kafka/input_output/Sentinel1_Observations.py index 5ce0087..2374871 100644 --- a/kafka/input_output/Sentinel1_Observations.py +++ b/kafka/input_output/Sentinel1_Observations.py @@ -165,7 +165,10 @@ def get_band_data(self, timestep, band): emulator = self.emulators[polarisation] # TODO read in angle of incidence from netcdf file # metadata['incidence_angle_deg'] = - sardata = SARdata(observations, uncertainty, mask, None, emulator) + fname = 'NETCDF:"{:s}":{:s}'.format(this_file, "theta") + obs_ptr = reproject_image(fname, self.state_mask) + metadata = {'incidence_angle': obs_ptr.ReadAsArray()} + sardata = SARdata(observations, uncertainty, mask, metadata, emulator) return sardata diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index a40dc1c..e2692f0 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -10,7 +10,7 @@ LOG = logging.getLogger(__name__) -def sar_observation_operator(x, polarisation): +def sar_observation_operator(x, theta, polarisation): """ For the sar_observation_operator a simple Water Cloud Model (WCM) is used @@ -51,8 +51,7 @@ def sar_observation_operator(x, polarisation): # conversion of incidence angle to radiant # the incidence angle itself should probably implemented in x) # TODO needs to come from the data - theta = np.deg2rad(23.) - theta = np.deg2rad(np.arange(1., 80.)) + theta = np.deg2rad(theta) # Simpler definition of cosine of theta mu = np.cos(theta) @@ -140,10 +139,11 @@ def create_sar_observation_operator(n_params, forward_model, metadata, if m: x0[i, :] = x_forecast[(n_params * i): (n_params*(i+1))] LOG.info("Running SAR forward model") + theta = metadata['incidence_angle'][mask[state_mask]] # Calls the run_emulator method that only does different vectors # It might be here that we do some sort of clustering - H0_, dH = forward_model(x0[mask[state_mask]], polarisation) + H0_, dH = forward_model(x0[mask[state_mask]], theta, polarisation) LOG.info("Storing emulators in H matrix") # This loop can be JIT'ed too From 092fb3c66861cb3d719a765c4fbf02f37dc5fb0b Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Wed, 17 Jan 2018 18:12:15 +0000 Subject: [PATCH 15/18] *Nearly* there with SAR Current version has had SAR coupled. One thing that I didn't envisage was that we need the Hessian of the RT model to fully calculate the non-linear terms in the uncertainty, and in this case, we do not have it available, so currently, the system breaks when it tries to calculate the Hessian correction using the GP's Hessian method... Other than that... 1. Observations are read and warped to match the state grid 2. The observation operator results in a linearised problem 3. The linear problem is solved 4. The linear part of the uncertainty is solved 5. It then crashes, but we know why! @McWhity This is where we left it --- kafka/__init__.py | 4 ++- kafka/inference/kf_tools.py | 5 +++- kafka/input_output/Sentinel1_Observations.py | 29 ++++++++++++++----- .../sar_forward_model.py | 25 ++++++++-------- 4 files changed, 41 insertions(+), 22 deletions(-) diff --git a/kafka/__init__.py b/kafka/__init__.py index c2aa732..687e357 100644 --- a/kafka/__init__.py +++ b/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 diff --git a/kafka/inference/kf_tools.py b/kafka/inference/kf_tools.py index cc7960a..8d6c091 100644 --- a/kafka/inference/kf_tools.py +++ b/kafka/inference/kf_tools.py @@ -20,7 +20,10 @@ def band_selecta(band): def hessian_correction_pixel(gp, x0, C_obs_inv, innovation, band, nparams): selecta = band_selecta(band) - ddH = gp.hessian(np.atleast_2d(x0[selecta])) + try: + ddH = gp.hessian(np.atleast_2d(x0[selecta])) + except AttributeError: + return 0. # In case the hessian is not implemented big_ddH = np.zeros((nparams, nparams)) for i, ii in enumerate(selecta): for j, jj in enumerate(selecta): diff --git a/kafka/input_output/Sentinel1_Observations.py b/kafka/input_output/Sentinel1_Observations.py index 2374871..c43bf4c 100644 --- a/kafka/input_output/Sentinel1_Observations.py +++ b/kafka/input_output/Sentinel1_Observations.py @@ -13,6 +13,8 @@ import osr +import scipy.sparse as sp + WRONG_VALUE = -999.0 # TODO tentative missing value SARdata = namedtuple('SARdata', @@ -59,14 +61,15 @@ def __init__(self, data_folder, state_mask, files = glob.glob(os.path.join(data_folder, '*.nc')) files.sort() self.state_mask = state_mask - self.dates = {} + 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[this_date] = fich - + self.dates.append(this_date) + self.date_data[this_date] = fich # 2. Store the emulator(s) self.emulators = emulators @@ -108,12 +111,13 @@ def _calculate_uncertainty(self, backscatter): """ # first approximation of uncertainty (1 dB) - self.unc = 1. + unc = backscatter*0.05 + # need to find a good way to calculate ENL # self.ENL = (backscatter.mean()**2) / (backscatter.std()**2) - return self.unc + return unc def _get_mask(self, backscatter): """ @@ -131,7 +135,8 @@ def _get_mask(self, backscatter): """ - mask = np.where(backscatter == WRONG_VALUE) + mask = np.ones_like(backscatter, dtype=np.bool) + mask[backscatter == WRONG_VALUE] = False return mask def get_band_data(self, timestep, band): @@ -156,19 +161,27 @@ def get_band_data(self, timestep, band): polarisation = 'VV' elif band == 1: polarisation = 'VH' - this_file = self.dates[timestep] + 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) metadata = {'incidence_angle': obs_ptr.ReadAsArray()} - sardata = SARdata(observations, uncertainty, mask, metadata, emulator) + sardata = SARdata(observations, R_mat_sp, mask, metadata, emulator) return sardata diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index e2692f0..c6780d8 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -67,9 +67,9 @@ def sar_observation_operator(x, theta, polarisation): raise ValueError('Only VV and VH polarisations available!') # Calculate Model - tau = np.exp(-2 * B / mu * x[0, :]) - sigma_veg = A * (x[0, :] ** E) * mu * (1 - tau) - sigma_surf = 10 ** ((C + D * x[1, :]) / 10.) + tau = np.exp(-2 * B / mu * x[:, 0]) + sigma_veg = A * (x[:, 0] ** E) * mu * (1 - tau) + sigma_surf = 10 ** ((C + D * x[:, 1]) / 10.) sigma_0 = sigma_veg + tau * sigma_surf #pdb.set_trace() @@ -78,14 +78,14 @@ def sar_observation_operator(x, theta, polarisation): grad = x*0 n_elems = x.shape[1] for i in range(n_elems): - tau_value = np.exp(-2 * B / mu * x[0, i]) - grad[0, i] = A * E * mu * (x[0, i] ** (E - 1)) * (1 - tau_value) + \ - 2 * A * B * (x[0, i] ** E) * tau_value - ( - (2 ** (1/10. * (C + D * x[1, i]) + 1)) * - (5 ** (1/10. * (C + D * x[1, i])) * B * tau_value) - ) / mu - grad[1, i] = D * np.log(10) * tau_value * 10 ** ( - 1/10. * (C + D * x[1, i]) - 1) + tau_value = np.exp(-2 * B / mu[i] * x[i, 0]) + grad[i, 0] = A * E * mu[i] * (x[i, 0] ** (E - 1)) * (1 - tau_value) + \ + 2 * A * B * (x[i, 0] ** E) * tau_value - ( + (2 ** (1/10. * (C + D * x[i, 1]) + 1)) * + (5 ** (1/10. * (C + D * x[i, 1])) * B * tau_value) + ) / mu[i] + grad[i, 1] = D * np.log(10) * tau_value * 10 ** ( + 1/10. * (C + D * x[i, 1]) - 1) # returned values are linear scaled not dB!!! # return sigma_0, grad, sigma_veg, sigma_surf, tau @@ -135,11 +135,12 @@ def create_sar_observation_operator(n_params, forward_model, metadata, polarisation = "VH" # This loop can be JIT'ed x0 = np.zeros((n_times, n_params)) + theta = np.zeros((n_times)) for i, m in enumerate(mask[state_mask].flatten()): if m: x0[i, :] = x_forecast[(n_params * i): (n_params*(i+1))] + theta[i] = 23.#metadata['incidence_angle'] LOG.info("Running SAR forward model") - theta = metadata['incidence_angle'][mask[state_mask]] # Calls the run_emulator method that only does different vectors # It might be here that we do some sort of clustering From 4847d69f4653a5960581c16c27e12661b4956066 Mon Sep 17 00:00:00 2001 From: Jose Gomez-Dans Date: Thu, 18 Jan 2018 13:50:13 +0000 Subject: [PATCH 16/18] Hessian fix updated for master --- kafka/inference/kf_tools.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/kafka/inference/kf_tools.py b/kafka/inference/kf_tools.py index 8d6c091..753c468 100644 --- a/kafka/inference/kf_tools.py +++ b/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: @@ -20,10 +25,7 @@ def band_selecta(band): def hessian_correction_pixel(gp, x0, C_obs_inv, innovation, band, nparams): selecta = band_selecta(band) - try: - ddH = gp.hessian(np.atleast_2d(x0[selecta])) - except AttributeError: - return 0. # In case the hessian is not implemented + ddH = gp.hessian(np.atleast_2d(x0[selecta])) big_ddH = np.zeros((nparams, nparams)) for i, ii in enumerate(selecta): for j, jj in enumerate(selecta): @@ -36,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 = [] @@ -62,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 From 4f6ca96f5f86fef11eacfbc12e45011e08aa9882 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20G=C3=B3mez-Dans?= Date: Thu, 18 Jan 2018 18:48:53 +0000 Subject: [PATCH 17/18] Fixed issue with gradient calculation Fixed issue with gradient calculation. Inversion of SAR data for one single time step has now been tested. The results appear plausible (e.g. the observations are closely fitted and the LAI and SM values are within the expected ranges). I haven't yet tested whether the uncertainty is sensible, but want to merge this back into master. --- .../sar_forward_model.py | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/kafka/observation_operators/sar_forward_model.py b/kafka/observation_operators/sar_forward_model.py index c6780d8..2cebea1 100644 --- a/kafka/observation_operators/sar_forward_model.py +++ b/kafka/observation_operators/sar_forward_model.py @@ -67,25 +67,25 @@ def sar_observation_operator(x, theta, polarisation): raise ValueError('Only VV and VH polarisations available!') # Calculate Model - tau = np.exp(-2 * B / mu * x[:, 0]) - sigma_veg = A * (x[:, 0] ** E) * mu * (1 - tau) - sigma_surf = 10 ** ((C + D * x[:, 1]) / 10.) + tau = np.exp(-2. * B / mu * x[:, 0]) + sigma_veg = A * (x[:, 0] ** E) * mu * (1. - tau) + sigma_surf = 10. ** ((C + D * x[:, 1]) / 10.) sigma_0 = sigma_veg + tau * sigma_surf #pdb.set_trace() # Calculate Gradient (grad has same dimension as x) grad = x*0 - n_elems = x.shape[1] + n_elems = x.shape[0] for i in range(n_elems): - tau_value = np.exp(-2 * B / mu[i] * x[i, 0]) - grad[i, 0] = A * E * mu[i] * (x[i, 0] ** (E - 1)) * (1 - tau_value) + \ - 2 * A * B * (x[i, 0] ** E) * tau_value - ( - (2 ** (1/10. * (C + D * x[i, 1]) + 1)) * - (5 ** (1/10. * (C + D * x[i, 1])) * B * tau_value) + tau_value = np.exp(-2. * B / mu[i] * x[i, 0]) + grad[i, 0] = A * E * mu[i] * (x[i, 0] ** (E - 1.)) * (1. - tau_value) + \ + 2. * A * B * (x[i, 0] ** E) * tau_value - ( + (2. ** (1/10. * (C + D * x[i, 1]) + 1.)) * + (5. ** (1/10. * (C + D * x[i, 1])) * B * tau_value) ) / mu[i] - grad[i, 1] = D * np.log(10) * tau_value * 10 ** ( - 1/10. * (C + D * x[i, 1]) - 1) + grad[i, 1] = D * np.log(10.) * tau_value * 10. ** ( + 1./10. * (C + D * x[i, 1]) - 1.) # returned values are linear scaled not dB!!! # return sigma_0, grad, sigma_veg, sigma_surf, tau @@ -149,6 +149,7 @@ def create_sar_observation_operator(n_params, forward_model, metadata, LOG.info("Storing emulators in H matrix") # This loop can be JIT'ed too n = 0 + for i, m in enumerate(mask[state_mask].flatten()): if m: H_matrix[i, (n_params * i): (n_params*(i+1))] = dH[n] From e684adb87ccea02790ed458b0ecf696f50a23e32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thomas=20Wei=C3=9F?= Date: Sun, 4 Feb 2018 23:47:40 +0100 Subject: [PATCH 18/18] uncomment parts that are not functioning and not needed for SAR --- kafka/input_output/Sentinel1_Observations.py | 7 +- kafka/input_output/observations.py | 180 +++++++++---------- 2 files changed, 95 insertions(+), 92 deletions(-) diff --git a/kafka/input_output/Sentinel1_Observations.py b/kafka/input_output/Sentinel1_Observations.py index c43bf4c..162e2bb 100644 --- a/kafka/input_output/Sentinel1_Observations.py +++ b/kafka/input_output/Sentinel1_Observations.py @@ -10,7 +10,7 @@ import gdal import numpy as np - +import pdb import osr import scipy.sparse as sp @@ -180,7 +180,10 @@ def get_band_data(self, timestep, band): # metadata['incidence_angle_deg'] = fname = 'NETCDF:"{:s}":{:s}'.format(this_file, "theta") obs_ptr = reproject_image(fname, self.state_mask) - metadata = {'incidence_angle': obs_ptr.ReadAsArray()} + 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 diff --git a/kafka/input_output/observations.py b/kafka/input_output/observations.py index 3b8a4fd..34adc1b 100644 --- a/kafka/input_output/observations.py +++ b/kafka/input_output/observations.py @@ -49,7 +49,7 @@ import numpy as np -from BRDF_descriptors import RetrieveBRDFDescriptors +# from BRDF_descriptors import RetrieveBRDFDescriptors #from kernels import Kernels @@ -209,92 +209,92 @@ def get_band_data(self, the_date, band_no): BHR = np.sum(BHR * to_NIR, axis=0) + a_to_NIR -class BHRObservations(RetrieveBRDFDescriptors): - def __init__(self, emulator, tile, mcd43a1_dir, - start_time, ulx=0, uly=0, dx=2400, dy=2400, end_time=None, - mcd43a2_dir=None): - """The class needs to locate the data granules. We assume that - these are available somewhere in the filesystem and that we can - index them by location (MODIS tile name e.g. "h19v10") and - time. The user can give a folder for the MCD43A1 and A2 granules, - and if the second is ignored, it will be assumed that they are - in the same folder. We also need a starting date (either a - datetime object, or a string in "%Y-%m-%d" or "%Y%j" format. If - the end time is not specified, it will be set to the date of the - latest granule found.""" - - # Call the constructor first - # Python2 - RetrieveBRDFDescriptors.__init__(self, tile, - mcd43a1_dir, start_time, end_time, - mcd43a2_dir) - # Python3 - # super().__init__(tile, mcd43a1_dir, start_time, end_time, - # mcd43a2_dir) - self._get_emulator(emulator) - self.dates = sorted(self.a1_granules.keys()) - self.dates = self.dates[::8] - a1_temp = {} - a2_temp = {} - for k in self.dates: - a1_temp[k] = self.a1_granules[k] - a2_temp[k] = self.a2_granules[k] - self.a1_granules = a1_temp - self.a2_granules = a2_temp - self.band_transfer = {0: "vis", - 1: "nir"} - self.ulx = ulx - self.uly = uly - self.dx = dx - self.dy = dy - - def define_output(self): - reference_fname = self.a1_granules[self.dates[0]] - g = gdal.Open('HDF4_EOS:EOS_GRID:' + - '"%s":MOD_Grid_BRDF:BRDF_Albedo_Parameters_vis' % - reference_fname) - proj = g.GetProjection() - geoT = np.array(g.GetGeoTransform()) - new_geoT = geoT*1. - new_geoT[0] = new_geoT[0] + self.ulx*new_geoT[1] - new_geoT[3] = new_geoT[3] + self.uly*new_geoT[5] - return proj, new_geoT.tolist() +# class BHRObservations(RetrieveBRDFDescriptors): +# def __init__(self, emulator, tile, mcd43a1_dir, +# start_time, ulx=0, uly=0, dx=2400, dy=2400, end_time=None, +# mcd43a2_dir=None): +# """The class needs to locate the data granules. We assume that +# these are available somewhere in the filesystem and that we can +# index them by location (MODIS tile name e.g. "h19v10") and +# time. The user can give a folder for the MCD43A1 and A2 granules, +# and if the second is ignored, it will be assumed that they are +# in the same folder. We also need a starting date (either a +# datetime object, or a string in "%Y-%m-%d" or "%Y%j" format. If +# the end time is not specified, it will be set to the date of the +# latest granule found.""" + +# # Call the constructor first +# # Python2 +# RetrieveBRDFDescriptors.__init__(self, tile, +# mcd43a1_dir, start_time, end_time, +# mcd43a2_dir) +# # Python3 +# # super().__init__(tile, mcd43a1_dir, start_time, end_time, +# # mcd43a2_dir) +# self._get_emulator(emulator) +# self.dates = sorted(self.a1_granules.keys()) +# self.dates = self.dates[::8] +# a1_temp = {} +# a2_temp = {} +# for k in self.dates: +# a1_temp[k] = self.a1_granules[k] +# a2_temp[k] = self.a2_granules[k] +# self.a1_granules = a1_temp +# self.a2_granules = a2_temp +# self.band_transfer = {0: "vis", +# 1: "nir"} +# self.ulx = ulx +# self.uly = uly +# self.dx = dx +# self.dy = dy + +# def define_output(self): +# reference_fname = self.a1_granules[self.dates[0]] +# g = gdal.Open('HDF4_EOS:EOS_GRID:' + +# '"%s":MOD_Grid_BRDF:BRDF_Albedo_Parameters_vis' % +# reference_fname) +# proj = g.GetProjection() +# geoT = np.array(g.GetGeoTransform()) +# new_geoT = geoT*1. +# new_geoT[0] = new_geoT[0] + self.ulx*new_geoT[1] +# new_geoT[3] = new_geoT[3] + self.uly*new_geoT[5] +# return proj, new_geoT.tolist() + +# def _get_emulator(self, emulator): +# if not os.path.exists(emulator): +# raise IOError("The emulator {} doesn't exist!".format(emulator)) +# # Assuming emulator is in an pickle file... +# self.emulator = cPickle.load(open(emulator, 'rb')) + +# def get_band_data(self, the_date, band_no): + +# to_BHR = np.array([1.0, 0.189184, -1.377622]) +# retval = self.get_brdf_descriptors(band_no, the_date) +# if retval is None: # No data on this date +# return None +# kernels, mask, qa_level = retval +# mask = mask[self.uly:(self.uly + self.dy), +# self.ulx:(self.ulx + self.dx)] +# qa_level = qa_level[self.uly:(self.uly + self.dy), +# self.ulx:(self.ulx + self.dx)] +# kernels = kernels[:, self.uly:(self.uly + self.dy), +# self.ulx:(self.ulx + self.dx)] +# bhr = np.where(mask, +# kernels * to_BHR[:, None, None], np.nan).sum(axis=0) +# R_mat = np.zeros_like(bhr) + +# R_mat[qa_level == 0] = np.maximum(2.5e-3, bhr[qa_level == 0] * 0.05) +# R_mat[qa_level == 1] = np.maximum(2.5e-3, bhr[qa_level == 1] * 0.07) +# 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() + +# bhr_data = BHR_data(bhr, mask, R_mat_sp, None, self.emulator) +# return bhr_data - def _get_emulator(self, emulator): - if not os.path.exists(emulator): - raise IOError("The emulator {} doesn't exist!".format(emulator)) - # Assuming emulator is in an pickle file... - self.emulator = cPickle.load(open(emulator, 'rb')) - def get_band_data(self, the_date, band_no): - - to_BHR = np.array([1.0, 0.189184, -1.377622]) - retval = self.get_brdf_descriptors(band_no, the_date) - if retval is None: # No data on this date - return None - kernels, mask, qa_level = retval - mask = mask[self.uly:(self.uly + self.dy), - self.ulx:(self.ulx + self.dx)] - qa_level = qa_level[self.uly:(self.uly + self.dy), - self.ulx:(self.ulx + self.dx)] - kernels = kernels[:, self.uly:(self.uly + self.dy), - self.ulx:(self.ulx + self.dx)] - bhr = np.where(mask, - kernels * to_BHR[:, None, None], np.nan).sum(axis=0) - R_mat = np.zeros_like(bhr) - - R_mat[qa_level == 0] = np.maximum(2.5e-3, bhr[qa_level == 0] * 0.05) - R_mat[qa_level == 1] = np.maximum(2.5e-3, bhr[qa_level == 1] * 0.07) - 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() - - bhr_data = BHR_data(bhr, mask, R_mat_sp, None, self.emulator) - return bhr_data - - class BHRObservationsTest(object): """A class to test BHR data "one pixel at a time". In essence, one only needs to define a self.dates dictionary (keys are datetime objects), and a 2 element @@ -307,8 +307,8 @@ def __init__(self, emulator, dates, vis_albedo, nir_albedo): self.data = {} for ii, the_date in enumerate(dates): self.data[the_date] = [vis_albedo[ii], nir_albedo[ii]] - - + + def _get_emulator(self, emulator): if not os.path.exists(emulator): raise IOError("The emulator {} doesn't exist!".format(emulator)) @@ -327,8 +327,8 @@ def get_band_data(self, the_date, band_no): bhr_data = BHR_data(bhr, mask, R_mat_sp, None, self.emulator) return bhr_data - - + + class KafkaOutputBHRTest(object): """A very simple class to output the state.""" def __init__(self): @@ -339,7 +339,7 @@ def dump_data(self, timestep, x_analysis, P_analysis, P_analysis_inv, for ii, param in enumerate(["w_vis", "x_vis", "a_vis", "w_nir", "x_nir", "a_nir", "TeLAI"]): solution[param] = x_analysis[ii::7] - self.output[timestep] = solution + self.output[timestep] = solution class KafkaOutput(object): """A very simple class to output the state."""