import numpy as np
from scipy.ndimage import zoom
array_test = zoom(np.arange(4*4).reshape(4,4),2,order=0)
array_test
array([[ 0, 0, 1, 1, 2, 2, 3, 3],
[ 0, 0, 1, 1, 2, 2, 3, 3],
[ 4, 4, 5, 5, 6, 6, 7, 7],
[ 4, 4, 5, 5, 6, 6, 7, 7],
[ 8, 8, 9, 9, 10, 10, 11, 11],
[ 8, 8, 9, 9, 10, 10, 11, 11],
[12, 12, 13, 13, 14, 14, 15, 15],
[12, 12, 13, 13, 14, 14, 15, 15]])
import xsar
import xsarsea # need pip install git+https://gitlab.ifremer.fr/sarlib/xsarsea
import time
import rasterio
import xarray as xr
array_test_dataarray = xr.DataArray(array_test, dims=("atrack", "xtrack"), coords={"atrack":range(0, 8, 1), "xtrack":range(0,8,1)})
array_test_dataarray = array_test_dataarray.chunk({"xtrack": 1, 'atrack':1})
array_test_dataarray
<xarray.DataArray (atrack: 8, xtrack: 8)> dask.array<xarray-<this-array>, shape=(8, 8), dtype=int64, chunksize=(1, 1), chunktype=numpy.ndarray> Coordinates: * atrack (atrack) int64 0 1 2 3 4 5 6 7 * xtrack (xtrack) int64 0 1 2 3 4 5 6 7
|
array([0, 1, 2, 3, 4, 5, 6, 7])
array([0, 1, 2, 3, 4, 5, 6, 7])
array_test_dataarray.values
array([[ 0, 0, 1, 1, 2, 2, 3, 3],
[ 0, 0, 1, 1, 2, 2, 3, 3],
[ 4, 4, 5, 5, 6, 6, 7, 7],
[ 4, 4, 5, 5, 6, 6, 7, 7],
[ 8, 8, 9, 9, 10, 10, 11, 11],
[ 8, 8, 9, 9, 10, 10, 11, 11],
[12, 12, 13, 13, 14, 14, 15, 15],
[12, 12, 13, 13, 14, 14, 15, 15]])
from scipy import signal
boundary = 'symm'
def resampl_oliv(image,reducx,reducy, mode='constant', constant_values=np.nan, **kwargs):
# reduc= coeff de reduc
nl,nc = image.shape
# shape is ceiled for potential padding
shape = [int(np.ceil(nl/reducy)), int(np.ceil(nc/reducx))]
# how many lines/rows to add
pad_width = [
(0, int(np.ceil(nl/reducy)) * reducy - nl ),
(0, int(np.ceil(nc/reducx)) * reducx - nc) ]
# pad
padded = np.pad(image, pad_width=pad_width, mode=mode, constant_values=constant_values, **kwargs)
# split by box
sh = shape[0],padded.shape[0]//shape[0],shape[1],padded.shape[1]//shape[1]
padded_bybox = padded.reshape(sh)
# return nanmean by boxes
return np.nanmean(np.nanmean(padded_bybox,axis=-1),axis=1)
def R2_oliv(image,reducx,reducy):
# reduc= coeffecifient de reduction
B2 = np.mat('[1,2,1; 2,4,2; 1,2,1]',float) *1/16
B4 = signal.convolve(B2,B2)
_image = signal.convolve2d(image,B4,'same', boundary=boundary)
num = signal.convolve2d(np.ones_like(_image),B4,'same', boundary=boundary)
image = _image/num
image = resampl_oliv(image.copy(),reducx,reducy)
_image = signal.convolve2d(image,B2,'same', boundary=boundary)
num = signal.convolve2d(np.ones_like(_image),B2,'same', boundary=boundary)
image = _image/num
return image
#original
def resampl(image,reducx,reducy):
# reduc= coeff de reduc
nl,nc = image.shape
shape = [nl//reducy, nc//reducx]
sh = shape[0],image.shape[0]//shape[0],shape[1],image.shape[1]//shape[1]
return image.reshape(sh).mean(-1).mean(1) # reduit et lissage par moyenne
def R2(image,reducx,reducy):
# reduc= coeffecifient de reduction
B2 = np.mat('[1,2,1; 2,4,2; 1,2,1]',float) *1/16
B4 = signal.convolve(B2,B2)
_image = signal.convolve2d(image,B4,'same', boundary=boundary)
num = signal.convolve2d(np.ones_like(_image),B4,'same', boundary=boundary)
image = _image/num
image = resampl(image.copy(),reducx,reducy)
_image = signal.convolve2d(image,B2,'same', boundary=boundary)
num = signal.convolve2d(np.ones_like(_image),B2,'same', boundary=boundary)
image = _image/num
return image
resampl(array_test,2,2)
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])
resampl_oliv(array_test,2,2)
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])
array_test_dataarray.coarsen({'atrack':2,'xtrack':2}, side='left' ,boundary='pad').mean().values
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])
R2(array_test,2,2)
array([[ 1.953125, 2.609375, 3.515625, 4.171875],
[ 4.578125, 5.234375, 6.140625, 6.796875],
[ 8.203125, 8.859375, 9.765625, 10.421875],
[10.828125, 11.484375, 12.390625, 13.046875]])
R2_oliv(array_test,2,2)
array([[ 1.953125, 2.609375, 3.515625, 4.171875],
[ 4.578125, 5.234375, 6.140625, 6.796875],
[ 8.203125, 8.859375, 9.765625, 10.421875],
[10.828125, 11.484375, 12.390625, 13.046875]])
xsarsea.R2(array_test_dataarray, {'xtrack':2, 'atrack':2}).values
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]])