diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 5f8ac4933a9f83beae20dc82e254dc55efa438c6..d9538d6763eca0ae8e20ffbd17249d8c3f0ac158 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -17,6 +17,8 @@ test_arosics: - pip install -U py_tools_ds -q - pip install -U geoarray -q + - pip install scipy # TODO remove as soon as CI runner is rebuilt + # run tests - make pytest diff --git a/arosics/Tie_Point_Grid.py b/arosics/Tie_Point_Grid.py index 4d05653621b15fb76154f21bbb8961a88c8322b7..1f675ac336f638b1b94104329bc3ebfc1261a596 100755 --- a/arosics/Tie_Point_Grid.py +++ b/arosics/Tie_Point_Grid.py @@ -36,6 +36,7 @@ import numpy as np from geopandas import GeoDataFrame from pandas import DataFrame, Series from shapely.geometry import Point +from scipy.interpolate import Rbf, RegularGridInterpolator # internal modules from .CoReg import COREG @@ -171,7 +172,7 @@ class Tie_Point_Grid(object): self.XY_points, self.XY_mapPoints = self._get_imXY__mapXY_points(self.grid_res) self._CoRegPoints_table = None # set by self.CoRegPoints_table self._GCPList = None # set by self.to_GCPList() - self.kriged = None # set by Raster_using_Kriging() + self.kriged = None # set by to_Raster_using_Kriging() @property def mean_x_shift_px(self): @@ -954,78 +955,27 @@ class Tie_Point_Grid(object): return out_GA - def to_Raster_using_Kriging(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None, - fName_out=None, tilepos=None, tilesize=500, mp=None): - - mp = False if self.CPUs == 1 else True - self._Kriging_sp(attrName, skip_nodata=skip_nodata, skip_nodata_col=skip_nodata_col, - outGridRes=outGridRes, fName_out=fName_out, tilepos=tilepos) - - # if mp: - # tilepositions = UTL.get_image_tileborders([tilesize,tilesize],self.tgt_shape) - # args_kwargs_dicts=[] - # for tp in tilepositions: - # kwargs_dict = {'skip_nodata':skip_nodata,'skip_nodata_col':skip_nodata_col,'outGridRes':outGridRes, - # 'fName_out':fName_out,'tilepos':tp} - # args_kwargs_dicts.append({'args':[attrName],'kwargs':kwargs_dict}) - # # self.kriged=[] - # # for i in args_kwargs_dicts: - # # res = self.Kriging_mp(i) - # # self.kriged.append(res) - # # print(res) - # - # with multiprocessing.Pool() as pool: - # self.kriged = pool.map(self.Kriging_mp,args_kwargs_dicts) - # pool.close() # needed to make coverage work in multiprocessing - # pool.join() - # else: - # self.Kriging_sp(attrName,skip_nodata=skip_nodata,skip_nodata_col=skip_nodata_col, - # outGridRes=outGridRes,fName_out=fName_out,tilepos=tilepos) - res = self.kriged if mp else None - return res - - def _Kriging_sp(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None, - fName_out=None, tilepos=None): - GDF = self.CoRegPoints_table - GDF2pass = GDF if not skip_nodata else GDF[GDF[skip_nodata_col] != self.outFillVal] - - X_coords, Y_coords, ABS_SHIFT = GDF2pass['X_MAP'], GDF2pass['Y_MAP'], GDF2pass[attrName] - - xmin, ymin, xmax, ymax = GDF2pass.total_bounds - - grid_res = outGridRes if outGridRes else int(min(xmax - xmin, ymax - ymin) / 250) - grid_x, grid_y = np.arange(xmin, xmax + grid_res, grid_res), np.arange(ymax, ymin - grid_res, -grid_res) - - # Reference: P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, - # (Cambridge University Press, 1997) 272 p. - from pykrige.ok import OrdinaryKriging - OK = OrdinaryKriging(X_coords, Y_coords, ABS_SHIFT, variogram_model='spherical', verbose=False) - zvalues, sigmasq = OK.execute('grid', grid_x, grid_y, backend='C', n_closest_points=12) - - if self.CPUs is None or self.CPUs > 1: - fName_out = fName_out if fName_out else \ - "Kriging__%s__grid%s_ws%s_%s.tif" % (attrName, self.grid_res, self.COREG_obj.win_size_XY, tilepos) - else: - fName_out = fName_out if fName_out else \ - "Kriging__%s__grid%s_ws%s.tif" % (attrName, self.grid_res, self.COREG_obj.win_size_XY) - path_out = get_generic_outpath(dir_out=self.dir_out, fName_out=fName_out) - # add a half pixel grid points are centered on the output pixels - xmin = xmin - grid_res / 2 - # ymin = ymin - grid_res / 2 - # xmax = xmax + grid_res / 2 - ymax = ymax + grid_res / 2 + def to_interpolated_raster(self, + metric: str = 'ABS_SHIFT', + method: str = 'Rbf' + ) -> np.ndarray: + TPGI = Tie_Point_Grid_Interpolator(self) - GeoArray(zvalues, - geotransform=(xmin, grid_res, 0, ymax, 0, -grid_res), - projection=self.COREG_obj.shift.prj).save(path_out) + return TPGI.interpolate(metric=metric, method=method) - return zvalues + def to_Raster_using_Kriging(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None, + fName_out=None, tilepos=None, tilesize=500, mp=None): - def _Kriging_mp(self, args_kwargs_dict): - args = args_kwargs_dict.get('args', []) - kwargs = args_kwargs_dict.get('kwargs', []) + TPGI = Tie_Point_Grid_Interpolator(self) + self.kriged = TPGI._interpolate_via_kriging(attrName, + skip_nodata=skip_nodata, + skip_nodata_col=skip_nodata_col, + outGridRes=outGridRes, + fName_out=fName_out, + tilepos=tilepos, + CPUs=self.CPUs) - return self._Kriging_sp(*args, **kwargs) + return self.kriged class Tie_Point_Refiner(object): @@ -1287,3 +1237,157 @@ class Tie_Point_Refiner(object): self.ransac_model_robust = model_robust return outseries + + +class Tie_Point_Grid_Interpolator(object): + def __init__(self, tiepointgrid: Tie_Point_Grid): + self.tpg = tiepointgrid + + def interpolate(self, metric: str, method: str = 'Rbf'): + if method not in ['Rbf', 'Kriging']: + raise ValueError(method) + + if method == 'Rbf': + return self._interpolate_via_rbf(metric=metric, outshape=self.tpg.shift.shape) + + else: + raise NotImplementedError() + + def _interpolate_via_rbf(self, metric: str, outshape: tuple): + tiepoints = self.tpg.CoRegPoints_table[self.tpg.CoRegPoints_table.OUTLIER.__eq__(False)].copy() + + rows = np.array(tiepoints.Y_IM) + cols = np.array(tiepoints.X_IM) + data = np.array(tiepoints[metric]) + + from time import time + t0 = time() + + # https://github.com/agile-geoscience/xlines/blob/master/notebooks/11_Gridding_map_data.ipynb + + # f = Rbf(cols, rows, data, function='linear') + # f = Rbf(cols, rows, data) + # data_full = f(*np.meshgrid(np.arange(outshape[1]), + # np.arange(outshape[0]))) + + # rows_lowres = np.arange(0, outshape[0] + 10, 10) + # cols_lowres = np.arange(0, outshape[1] + 10, 10) + rows_lowres = np.arange(0, outshape[0] + 5, 5) + cols_lowres = np.arange(0, outshape[1] + 5, 5) + f = Rbf(cols, rows, data) + data_interp_lowres = f(*np.meshgrid(cols_lowres, rows_lowres)) + + # https://stackoverflow.com/questions/24978052/interpolation-over-regular-grid-in-python + # from sklearn.gaussian_process import GaussianProcess + # gp = GaussianProcess(theta0=0.1, thetaL=.001, thetaU=1., nugget=0.01) + # gp.fit(X=np.column_stack([rr[vals], cc[vals]]), y=M[vals]) + # rr_cc_as_cols = np.column_stack([rr.flatten(), cc.flatten()]) + # interpolated = gp.predict(rr_cc_as_cols).reshape(M.shape) + + RGI = RegularGridInterpolator(points=[cols_lowres, rows_lowres], + values=data_interp_lowres.T, # must be in shape [x, y] + method='linear', + bounds_error=False) + rows_full = np.arange(outshape[0]) + cols_full = np.arange(outshape[1]) + data_full = RGI(np.dstack(np.meshgrid(cols_full, rows_full))) + + print('interpolation runtime: %.2fs' % (time() - t0)) + + # from matplotlib import pyplot as plt + # plt.figure() + # im = plt.imshow(data_full) + # plt.colorbar(im) + # plt.scatter(cols, rows, c=data, edgecolors='black') + # plt.title(metric) + # plt.show() + + return data_full + + def _interpolate_via_kriging(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None, + fName_out=None, tilepos=None, tilesize=500, CPUs=1): + + CPUs = 1 + if CPUs == 1: + data_full = self._Kriging_sp(attrName, skip_nodata=skip_nodata, skip_nodata_col=skip_nodata_col, + outGridRes=outGridRes, fName_out=fName_out, tilepos=tilepos) + + else: + raise NotImplementedError() + + # tilepositions = UTL.get_image_tileborders([tilesize, tilesize], self.tgt_shape) + # args_kwargs_dicts = [] + # for tp in tilepositions: + # kwargs_dict = {'skip_nodata': skip_nodata, + # 'skip_nodata_col': skip_nodata_col, + # 'outGridRes': outGridRes, + # 'fName_out': fName_out, + # 'tilepos': tp} + # args_kwargs_dicts.append({'args': [attrName], 'kwargs': kwargs_dict}) + # + # # data_full=[] + # # for i in args_kwargs_dicts: + # # res = self.Kriging_mp(i) + # # data_full.append(res) + # # print(res) + # + # with multiprocessing.Pool() as pool: + # data_full = pool.map(self._Kriging_mp, args_kwargs_dicts) + # pool.close() # needed to make coverage work in multiprocessing + # pool.join() + + return data_full + + def _Kriging_sp(self, attrName, skip_nodata=1, skip_nodata_col='ABS_SHIFT', outGridRes=None, + fName_out=None, tilepos=None): + # Reference: P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, + # (Cambridge University Press, 1997) 272 p. + from pykrige.ok import OrdinaryKriging + + GDF = self.tpg.CoRegPoints_table + GDF2pass = GDF if not skip_nodata else GDF[GDF[skip_nodata_col] != self.tpg.outFillVal] + + X_coords, Y_coords, ABS_SHIFT = GDF2pass['X_MAP'], GDF2pass['Y_MAP'], GDF2pass[attrName] + + xmin, ymin, xmax, ymax = GDF2pass.total_bounds + + grid_res = outGridRes or int(min(xmax - xmin, ymax - ymin) / 250) + grid_x = np.arange(xmin, xmax + grid_res, grid_res) + grid_y = np.arange(ymax, ymin - grid_res, -grid_res) + + OK = OrdinaryKriging(X_coords, Y_coords, ABS_SHIFT, + variogram_model='spherical', + verbose=False) + zvalues, sigmasq = OK.execute('grid', grid_x, grid_y, + backend='C', + n_closest_points=12) + + # FIXME: revise this save-routine + if self.tpg.CPUs is None or self.tpg.CPUs > 1: + fName_out = fName_out if fName_out else \ + "Kriging__%s__grid%s_ws%s_%s.tif" \ + % (attrName, self.tpg.grid_res, self.tpg.COREG_obj.win_size_XY, tilepos) + else: + fName_out = fName_out if fName_out else \ + "Kriging__%s__grid%s_ws%s.tif" \ + % (attrName, self.tpg.grid_res, self.tpg.COREG_obj.win_size_XY) + + path_out = get_generic_outpath(dir_out=self.tpg.dir_out, fName_out=fName_out) + + # add half a pixel as grid points are centered on the output pixels + xmin = xmin - grid_res / 2 + # ymin = ymin - grid_res / 2 + # xmax = xmax + grid_res / 2 + ymax = ymax + grid_res / 2 + + GeoArray(zvalues, + geotransform=(xmin, grid_res, 0, ymax, 0, -grid_res), + projection=self.tpg.COREG_obj.shift.prj).save(path_out) + + return zvalues + + def _Kriging_mp(self, args_kwargs_dict): + args = args_kwargs_dict.get('args', []) + kwargs = args_kwargs_dict.get('kwargs', []) + + return self._Kriging_sp(*args, **kwargs) diff --git a/requirements.txt b/requirements.txt index eda4876223969abeb600224b537e39c7887f119a..55c917837f580b71ec7a6ba7e558ced26866ca32 100644 --- a/requirements.txt +++ b/requirements.txt @@ -15,4 +15,5 @@ pykrige pyproj>2.2.0 py_tools_ds>=0.18.0 scikit-image>=0.16.0 +scipy shapely diff --git a/setup.py b/setup.py index 9720a2371c793a971b8ebe2509ffe207c9f3d753..ec7e1688500b307d52fe01b6797a15dcc6c2f318 100644 --- a/setup.py +++ b/setup.py @@ -57,6 +57,7 @@ req = [ 'pyproj>2.2.0', 'py_tools_ds>=0.18.0', 'scikit-image>=0.16.0', + 'scipy', 'shapely', ] diff --git a/tests/CI_docker/context/environment_arosics.yml b/tests/CI_docker/context/environment_arosics.yml index 9ab3e8dd118c024377b1db759fbce43b8a0a066e..1006bf33c04198e1100907205ce6f42406388ce4 100644 --- a/tests/CI_docker/context/environment_arosics.yml +++ b/tests/CI_docker/context/environment_arosics.yml @@ -29,6 +29,7 @@ dependencies: - dill - folium>=0.6.0,!=0.12.0 - geojson + - scipy # doc requirements - flake8 diff --git a/tests/test_tie_point_grid.py b/tests/test_tie_point_grid.py index 3773ee703ede0d70a929a64c5252f89281c909dd..14e08a557d86e24251dad70b4889c1126b571e60 100644 --- a/tests/test_tie_point_grid.py +++ b/tests/test_tie_point_grid.py @@ -34,6 +34,8 @@ import shutil import warnings # custom +import numpy as np + from .cases import test_cases from arosics import COREG_LOCAL, Tie_Point_Grid @@ -125,6 +127,11 @@ class Test_Tie_Point_Grid(unittest.TestCase): self.TPG.to_vectorfield(outpath, fmt='ENVI', mode='uv') self.assertTrue(os.path.isfile(outpath)) + def test_interpolate_to_raster(self): + arr_interp = self.TPG.to_interpolated_raster('ABS_SHIFT', 'Rbf') + + self.assertIsInstance(arr_interp, np.ndarray) + def test_to_Raster_using_Kriging(self): if find_loader('pykrige.ok'): with tempfile.TemporaryDirectory() as tmpdir: