From 648b59c91190fecacf8fd803927b4c10296b0e6d Mon Sep 17 00:00:00 2001 From: "stephen.worsley" Date: Tue, 29 Mar 2022 15:50:49 +0100 Subject: [PATCH] add explanation to the mathematics. --- esmf_regrid/_esmf_sdo.py | 42 ++++++++++++++++++++++++++++------------ 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/esmf_regrid/_esmf_sdo.py b/esmf_regrid/_esmf_sdo.py index 613cbce5..3f6df614 100644 --- a/esmf_regrid/_esmf_sdo.py +++ b/esmf_regrid/_esmf_sdo.py @@ -397,20 +397,19 @@ def _collapse_weights(self, tgt): # The column indices represent each of the cells in the refined grid. column_indices = np.arange(self._extended_size) - # The row indices represent the cells of the represented grid. These are broadcast - # so that each index coincides with all column indices of the refined cells which - # make the represented cell is split into. + # The row indices represent the cells of the unrefined grid. These are broadcast + # so that each row index coincides with all column indices of the refined cells + # which the unrefined cell is split into. if self.lat_expansion > 1: - row_indices = np.empty( - [ - self.n_lons_orig, - self.n_lats_orig * self.lat_expansion, - ] - ) - row_indices[:] = np.arange(self.n_lons_orig * self.n_lats_orig)[ - :, np.newaxis - ] + # The latitudes are expanded only in the case where there is one latitude + # bound from -90 to 90. In this case, there is no longitude expansion. + row_indices = np.empty([self.n_lons_orig, self.lat_expansion]) + row_indices[:] = np.arange(self.n_lons_orig)[:, np.newaxis] else: + # The row indices are broadcast across a dimension representing the expansion + # of the longitude. Each row index is broadcast and flattened so that all the + # row indices representing the unrefined cell match up with the column indices + # representing the refined cells it is split into. row_indices = np.empty( [self.n_lons_orig, self.lon_expansion, self.n_lats_orig] ) @@ -427,9 +426,28 @@ def _collapse_weights(self, tgt): shape=matrix_shape, ) if tgt: + # When the RefinedGridInfo is the target of the regridder, we want to take + # the average of the weights of each refined target cell. This is because + # these weights represent the proportion of area of the target cells which + # is covered by a given source cell. Since the refined cells are divided in + # such a way that they have equal area, the weights for the unrefined cells + # can be reconstructed by taking an average. This is done via matrix + # multiplication, with the returned matrix pre-multiplying the weight matrix + # so that it operates on the rows of the weight matrix (representing the + # target cells). At this point the returned matrix consists of ones, so we + # divided by the number of refined cells per unrefined cell. refinement_weights = refinement_weights / ( self.lon_expansion * self.lat_expansion ) else: + # When the RefinedGridInfo is the source of the regridder, we want to take + # the sum of the weights of each refined target cell. This is because those + # weights represent the proportion of the area of a given target cell which + # is covered by each refined source cell. The total proportion covered by + # each unrefined source cell is then the sum of the weights from each of its + # refined cells. This sum is done by matrix multiplication, the returned + # matrix post-multiplying the weight matrix so that it operates on the columns + # of the weight matrix (representing the source cells). In order for the + # post-multiplication to work, the returned matrix must be transposed. refinement_weights = refinement_weights.T return refinement_weights