Commit 3e56cb8f authored by Jeff Piollé's avatar Jeff Piollé

Merge branch 'develop' of https://git.cersat.fr/cerbere/cerbere into develop

parents e9311903 bfeb12e9
......@@ -58,8 +58,8 @@ class AbstractFeature(object):
# check type of fields
if fields is not None:
if isinstance(fields, dict):
logging.warning("Usage of dict for 'fields' is deprecated")
if not hasattr(fields,'_OrderedDict__root'):# _OrderedDict__marker
logging.warning("Usage of dict for 'fields' is deprecated (%s)",fields)
self._fields = OrderedDict(fields)
elif not isinstance(fields, OrderedDict):
raise TypeError(
......
......@@ -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]
......@@ -483,7 +598,7 @@ class SAFEMSIL1CGranuleFile(AbstractMapper):
yindmin = int(np.maximum(0, np.floor(ymaxind)))
yindmax = int(np.minimum(self._attrs['ncols'], np.ceil(yminind)))
# Create tight slices (equivalent to view slices in other mappers)
slices = [slice(yindmin, yindmax, 1), slice(xindmin, xindmax, 1)]
slices = (slice(yindmin, yindmax, 1), slice(xindmin, xindmax, 1))
dimsizes = [self._attrs['ncols'], self._attrs['nrows']]
self._tight_slices = Slices(slices, dimsizes)
# Update attributes
......@@ -629,9 +744,9 @@ class SAFEMSIL1CGranuleFile(AbstractMapper):
epsg = self.read_global_attribute('horizontal_cs_code')
proj = pyproj.Proj(init=epsg)
shp = slices.shape()
y = np.tile(self.read_values('y', slices=[slices[0]])[:, np.newaxis],
y = np.tile(self.read_values('y', slices=(slices[0],))[:, np.newaxis],
(1, shp[1]))
x = np.tile(self.read_values('x', slices=[slices[1]])[np.newaxis, :],
x = np.tile(self.read_values('x', slices=(slices[1],))[np.newaxis, :],
(shp[0], 1))
lon, lat = proj(x, y, inverse=True)
if fieldname == 'lon':
......@@ -712,8 +827,8 @@ class SAFEMSIL1CGranuleFile(AbstractMapper):
nodatavalue=float(fillvalue))
indmin = [i.min() for i in ind]
indmax = [i.max() + 1 for i in ind]
rel_slices = Slices([slice(indmin[0], indmax[0]),
slice(indmin[1], indmax[1])],
rel_slices = Slices((slice(indmin[0], indmax[0]),
slice(indmin[1], indmax[1])),
det_ind.shape)
abs_slices = rel_slices.absolute_slices(slices)
target_ds = self._get_gdal_dataset_from_slices(abs_slices,
......@@ -954,8 +1069,8 @@ class SAFEMSIL1CGranuleFile(AbstractMapper):
tuple: bbox expressed as (lonmin, latmin, lonmax, latmax)
"""
dimsizes = self.get_full_dimensions('lon').values()
slices = [slice(None, None, dimsizes[0] - 1),
slice(None, None, dimsizes[1] - 1)]
slices = (slice(None, None, dimsizes[0] - 1),
slice(None, None, dimsizes[1] - 1))
lon = self.read_values('lon', slices=slices)
lat = self.read_values('lat', slices=slices)
return (lon.min(), lat.min(), lon.max(), lat.max())
......
......@@ -283,9 +283,9 @@ class SAFEMSIL1CStitchedFile(AbstractMapper):
epsg = self.read_global_attribute('horizontal_cs_code')
proj = pyproj.Proj(init=epsg)
shp = slices.shape()
y = np.tile(self.read_values('y', slices=[slices[0]])[:, np.newaxis],
y = np.tile(self.read_values('y', slices=(slices[0],))[:, np.newaxis],
(1, shp[1]))
x = np.tile(self.read_values('x', slices=[slices[1]])[np.newaxis, :],
x = np.tile(self.read_values('x', slices=(slices[1],))[np.newaxis, :],
(shp[0], 1))
lon, lat = proj(x, y, inverse=True)
if fieldname == 'lon':
......@@ -500,8 +500,8 @@ class SAFEMSIL1CStitchedFile(AbstractMapper):
tuple: bbox expressed as (lonmin, latmin, lonmax, latmax)
"""
dimsizes = self.get_full_dimensions('lon').values()
slices = [slice(None, None, dimsizes[0] - 1),
slice(None, None, dimsizes[1] - 1)]
slices = (slice(None, None, dimsizes[0] - 1),
slice(None, None, dimsizes[1] - 1))
lon = self.read_values('lon', slices=slices)
lat = self.read_values('lat', slices=slices)
return (lon.min(), lat.min(), lon.max(), lat.max())
......
......@@ -257,7 +257,7 @@ class Slices(tuple):
raise Exception('Unexpected type for dimsizes.')
# Check slices input type
if slices is None:
_slices = [slice(None)] * len(dimsizes)
_slices = [slice(None)] * len(_dimsizes)
elif isinstance(slices, (slice, int)):
_slices = [slices]
elif isinstance(slices, tuple):
......@@ -602,17 +602,21 @@ class Slices(tuple):
return None, None
outmin = int(np.maximum(0, outmin))
outmax = int(np.minimum(abssize - 1, outmax))
# TMP check
if outmin < 0 or outmin >= abssize or outmax < 0 or outmax >= abssize:
raise Exception()
# Check outmin / outmax
assert outmin >= 0 and outmin < abssize and \
outmax >= 0 and outmax < abssize
# Find insub min and max
inmin = int(absmin + outmin * abs(absstep) - submin)
inmax = int(absmin + outmax * abs(absstep) - submin)
# TMP check
if inmin < 0 or inmin >= subsize or inmax < 0 or inmax >= subsize:
raise Exception()
# Check inmin / inmax
assert inmin >= 0 and inmin < subsize and \
inmax >= 0 and inmax < subsize
# Make slices
if absstep > 0:
if self._is_index[i]:
assert outmin == outmax and inmin == inmax
outsub_slices.append(outmin)
insub_slices.append(inmin)
elif absstep > 0:
outsub_slices.append(slice(outmin, outmax+1, 1))
insub_slices.append(slice(inmin, inmax+1, absstep))
else:
......@@ -623,10 +627,9 @@ class Slices(tuple):
else:
instop = inmin - 1
insub_slices.append(slice(inmax, instop, absstep))
res = (Slices(outsub_slices, shp), Slices(insub_slices, sub_size))
# TMP check
for size1, size2 in zip(res[0].shape(), res[1].shape()):
if size1 != size2:
raise Exception()
res = (Slices(tuple(outsub_slices), shp),
Slices(tuple(insub_slices), sub_size))
# Check shapes
assert res[0].shape() == res[1].shape()
return res
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