Commit 90965fbc authored by Gilles Guitton's avatar Gilles Guitton

- mapper SAFEMSIL1CGranuleFile : improve viewing angles interpolation by...

- mapper SAFEMSIL1CGranuleFile : improve viewing angles interpolation by extrapolating original low res grid.
parent 12a02e25
......@@ -14,6 +14,7 @@ from osgeo import gdal, ogr
from datetime import datetime
import tempfile
import shutil
from scipy.ndimage.morphology import binary_dilation
from cerbere.mapper.abstractmapper import AbstractMapper
from cerbere import READ_ONLY
......@@ -182,7 +183,121 @@ class SAFEMSIL1CGranuleFile(AbstractMapper):
ydim = grid['row_step']
if np.sign(ydim) != np.sign(self.read_global_attribute('ydim')):
ydim = -ydim
values = grid['values']
values = grid['values'].copy()
isfin = np.isfinite(values)
if 'viewing' in angle_name and isfin.any():
#if self._extrapolate == True and 'viewing' in angle_name and isfin.any():
# Extrapolate
max_deg = 3
deg = [min([np.sum(isfin.sum(axis=1) > 0) - 1, max_deg]),
min([np.sum(isfin.sum(axis=0) > 0) - 1, max_deg])]
indfin = np.where(isfin)
vdr = np.polynomial.polynomial.polyvander2d(indfin[0],
indfin[1],
deg)
coeff, rs, rk, s = np.linalg.lstsq(vdr, values[indfin])
coeff = coeff.reshape((deg[0] + 1, deg[1] + 1))
# print vdr.shape, rs, rk, s
# print coeff
struct = np.ones((3, 3), dtype='bool')
dilisfin = binary_dilation(isfin, structure=struct, border_value=0)
ind2fill = np.where(~isfin & dilisfin)
values[ind2fill] = np.polynomial.polynomial.polyval2d(
ind2fill[0], ind2fill[1], coeff
)
isfin = np.isfinite(values)
## Extrapolate test
# print angle_name
# max_deg = 5
# coeff = np.zeros((max_deg + 1, max_deg + 1))
# counti = isfin.sum(axis=1).reshape((-1, 1))
# countj = isfin.sum(axis=0).reshape((1, -1))
# for i in range(0, min([max_deg + 1, np.sum(counti > 0)])):
# for j in range(0, min([max_deg + 1, np.sum(countj > 0)])):
# indcount = np.where((counti > j) & (countj > i))
# if np.unique(indcount[0]).size > i and \
# np.unique(indcount[1]).size > j:
# coeff[i, j] = 1
# indcoeff = np.where(coeff > 0)
# indfin = np.where(isfin)
# g = np.zeros((indfin[0].size, indcoeff[0].size))
# for k, (i, j) in enumerate(zip(indcoeff[0], indcoeff[1])):
# g[:, k] = indfin[0] ** i * indfin[1] ** j
# c, rs, rk, s = np.linalg.lstsq(g, values[indfin])
# print g.shape, rs, rk, s
# coeff[indcoeff] = c
# print coeff
# coeff = coeff[0:indcoeff[0].max() + 1, 0:indcoeff[1].max() + 1]
# struct = np.ones((3, 3), dtype='bool')
# dilisfin = binary_dilation(isfin, structure=struct, border_value=0)
# ind2fill = np.where(~isfin & dilisfin)
# values[ind2fill] = np.polynomial.polynomial.polyval2d(
# ind2fill[0], ind2fill[1], coeff
# )
# isfin = np.isfinite(values)
## Extrapolate test
# isfin = np.isfinite(values)
# cnty = isfin.sum(axis=0)
# for x in np.where((cnty >= 2) & (cnty < isfin.shape[0]))[0]:
# indy = np.where(isfin[:, x])[0]
# ny = indy.size
# y0, y1 = indy[0], indy[-1]
# if ny != (y1 - y0 + 1):
# raise Exception('Holes in angle grid ?')
# elif ny == 2:
# deg = 1
# elif ny == 3:
# deg = 2
# else:
# deg = 3
# coef = np.polyfit(indy, values[indy, x], deg)
# if y0 != 0:
# values[y0 - 1, x] = np.polyval(coef, y0 - 1)
# if y1 != values.shape[0] - 1:
# values[y1 + 1, x] = np.polyval(coef, y1 + 1)
# # if y0 != 0:
# # values[y0 - 1, x] = 2 * values[y0, x] - values[y0 + 1, x]
# # if y1 != values.shape[0] - 1:
# # values[y1 + 1, x] = 2 * values[y1, x] - values[y1 - 1, x]
# isfin = np.isfinite(values)
# cntx = isfin.sum(axis=1)
# for y in np.where((cntx >= 2) & (cntx < isfin.shape[1]))[0]:
# indx = np.where(isfin[y, :])[0]
# nx = indx.size
# x0, x1 = indx[0], indx[-1]
# if nx != (x1 - x0 + 1):
# raise Exception('Holes in angle grid ?')
# elif nx == 2:
# deg = 1
# elif nx == 3:
# deg = 2
# else:
# deg = 3
# coef = np.polyfit(indx, values[y, indx], deg)
# if x0 != 0:
# values[y, x0 - 1] = np.polyval(coef, x0 - 1)
# if x1 != values.shape[1] - 1:
# values[y, x1 + 1] = np.polyval(coef, x1 + 1)
# # if x0 != 0:
# # values[y, x0 - 1] = 2 * values[y, x0] - values[y, x0 + 1]
# # if x1 != values.shape[1] - 1:
# # values[y, x1 + 1] = 2 * values[y, x1] - values[y, x1 - 1]
## Check extrapolation
# import matplotlib.pyplot as plt
# vmin, vmax = np.nanmin(values), np.nanmax(values)
# plt.figure() ; plt.imshow(grid['values'], interpolation='nearest', vmin=vmin, vmax=vmax) ; plt.colorbar()
# plt.figure() ; plt.imshow(values, interpolation='nearest', vmin=vmin, vmax=vmax) ; plt.colorbar()
# plt.show()
# plt.figure()
# for i in range(values.shape[0]):
# plt.plot(values[i, :], '+-b')
# plt.plot(grid['values'][i, :], '+-r')
# plt.figure()
# for i in range(values.shape[1]):
# plt.plot(values[:, i], '+-b')
# plt.plot(grid['values'][:, i], '+-r')
# plt.show()
#import pdb ; pdb.set_trace()
drv = gdal.GetDriverByName('MEM')
dset = drv.Create(angle_name, values.shape[1], values.shape[0], 1,
gdal.GDT_Float32)
......@@ -193,7 +308,7 @@ class SAFEMSIL1CGranuleFile(AbstractMapper):
dset.SetGeoTransform(geotf)
if nodatavalue is not None:
dset.GetRasterBand(1).SetNoDataValue(nodatavalue)
values[np.where(~np.isfinite(values))] = nodatavalue
values[~isfin] = nodatavalue
dset.GetRasterBand(1).WriteArray(values)
self._angle_handlers[angle_name] = dset
return self._angle_handlers[angle_name]
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment