Skip to content

cerbere
cerbere
Commit 2204fc4f authored by Gilles Guitton's avatar Gilles Guitton
Browse files
Options

Merge branch 'feature/sentinel2' into develop

parents bc66c927 90965fbc
Hide whitespace changes
Inline Side-by-side
Showing with 128 additions and 13 deletions
cerbere/mapper/safemsil1cgranulefile.py
View file @ 2204fc4f
......@@ -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())
......
cerbere/mapper/safemsil1cstitchedfile.py
View file @ 2204fc4f
......@@ -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())
......
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