From d6d0054fd168db4c0b2827b203ce76b7c69b1515 Mon Sep 17 00:00:00 2001 From: Daniel Scheffler Date: Tue, 14 Sep 2021 13:32:29 +0200 Subject: [PATCH 1/2] Implemented automatic reprojection in case of unequal projections (not yet fully working). Signed-off-by: Daniel Scheffler --- arosics/CoReg.py | 68 +++++++++++++++++++++++++----------------------- 1 file changed, 36 insertions(+), 32 deletions(-) diff --git a/arosics/CoReg.py b/arosics/CoReg.py index d27d301..d0600ef 100755 --- a/arosics/CoReg.py +++ b/arosics/CoReg.py @@ -405,12 +405,25 @@ class COREG(object): self.success = None # default self.deshift_results = None # set by self.correct_shifts() - # try: gdal.AllRegister() + + # TODO: wrap this into an image pre-processor class + self.ref = GeoArray(im_ref, nodata=nodata[0], progress=progress) + self.shift = GeoArray(im_tgt, nodata=nodata[0], progress=progress) + self.ref.band4match = r_b4match - 1 + self.shift.band4match = s_b4match - 1 + self.grid2use = 'ref' if self.shift.xgsd <= self.ref.xgsd else 'shift' + + self._equalize_projections() self._check_and_handle_metaRotation() - self._get_image_params() + + self.params['im_ref'] = self.ref + self.params['im_tgt'] = self.shift + self.ref = GeoArray_CoReg(self.params, 'ref') + self.shift = GeoArray_CoReg(self.params, 'shift') + self._set_outpathes(im_ref, im_tgt) - self.grid2use = 'ref' if self.shift.xgsd <= self.ref.xgsd else 'shift' + if self.v: print('resolutions: ', self.ref.xgsd, self.shift.xgsd) @@ -421,9 +434,6 @@ class COREG(object): write_shp(os.path.join(self.path_verbose_out, 'poly_im2shift.shp'), self.shift.poly, self.shift.prj) write_shp(os.path.join(self.path_verbose_out, 'overlap_poly.shp'), self.overlap_poly, self.ref.prj) - # FIXME: transform_mapPt1_to_mapPt2(im2shift_center_map, ds_imref.GetProjection(), ds_im2shift.GetProjection()) - # FIXME später basteln für den fall, dass projektionen nicht gleich sind - # get_clip_window_properties self._get_opt_winpos_winsize() if not self.q: @@ -436,8 +446,6 @@ class COREG(object): self.success = False if self.success is False or not self.matchBox.boxMapYX else None - # except BaseException: - self._coreg_info = None # private attribute to be filled by self.coreg_info property def _handle_error(self, error, warn=False, warnMsg=None): @@ -541,23 +549,6 @@ class COREG(object): if self.path_verbose_out and not os.path.isdir(self.path_verbose_out): os.makedirs(self.path_verbose_out) - def _get_image_params(self) -> None: - self.ref = GeoArray_CoReg(self.params, 'ref') - self.shift = GeoArray_CoReg(self.params, 'shift') - - if not prj_equal(self.ref.prj, self.shift.prj): - from pyproj import CRS - - crs_ref = CRS.from_user_input(self.ref.prj) - crs_shift = CRS.from_user_input(self.shift.prj) - - name_ref, name_shift = \ - (crs_ref.name, crs_shift.name) if not crs_ref.name == crs_shift.name else (crs_ref.srs, crs_shift.srs) - - raise RuntimeError( - 'Input projections are not equal. Different projections are currently not supported. ' - 'Got %s vs. %s.' % (name_ref, name_shift)) - def _get_overlap_properties(self) -> None: with warnings.catch_warnings(): # already warned in GeoArray_CoReg.__init__() @@ -577,25 +568,38 @@ class COREG(object): assert px_covered > 16 * 16, \ 'Overlap area covers only %s pixels. At least 16*16 pixels are needed.' % px_covered + def _equalize_projections(self): + if not prj_equal(self.ref.prj, self.shift.prj): + from pyproj import CRS + + crs_ref = CRS.from_user_input(self.ref.prj) + crs_shift = CRS.from_user_input(self.shift.prj) + + name_ref, name_shift = \ + (crs_ref.name, crs_shift.name) if not crs_ref.name == crs_shift.name else (crs_ref.srs, crs_shift.srs) + + if not self.q: + print("Projections of reference and target image are unequal (%s vs. %s) " + "and will be equalized at runtime." % (name_ref, name_shift)) + + self.equalize_pixGrids() + def _check_and_handle_metaRotation(self): """Check if the provided input data have a metadata rotation and if yes, correct it. In case there is a rotation, the GDAL GeoTransform is not 0 at positions 2 or 4. So far, AROSICS does not handle such rotations, so the resampling is needed to make things work. """ - gA_ref = GeoArray(self.params['im_ref']) - gA_tgt = GeoArray(self.params['im_tgt']) - msg = 'The %s image needs to be resampled because it has a row/column rotation in ' \ 'its map info which is not handled by AROSICS.' - if GEO.has_metaRotation(gA_ref): + if GEO.has_metaRotation(self.ref): warnings.warn(msg % 'reference') - self.params['im_ref'] = GEO.remove_metaRotation(gA_ref) + self.ref = GEO.remove_metaRotation(self.ref) - if GEO.has_metaRotation(gA_tgt): + if GEO.has_metaRotation(self.shift): warnings.warn(msg % 'target') - self.params['im_tgt'] = GEO.remove_metaRotation(gA_tgt) + self.shift = GEO.remove_metaRotation(self.shift) @property def are_pixGrids_equal(self): -- GitLab From de858d505c31e78a94329031d37c3e29686fce81 Mon Sep 17 00:00:00 2001 From: Daniel Scheffler Date: Tue, 21 Sep 2021 21:51:41 +0200 Subject: [PATCH 2/2] Fixed message of assertion. Signed-off-by: Daniel Scheffler --- arosics/CoReg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/arosics/CoReg.py b/arosics/CoReg.py index d0600ef..3df1d33 100755 --- a/arosics/CoReg.py +++ b/arosics/CoReg.py @@ -94,7 +94,7 @@ class GeoArray_CoReg(GeoArray): assert self.bands >= self.band4match + 1 >= 1, \ "The %s has %s %s. So its band number to match must be %s%s. Got %s." \ % (self.imName, self.bands, 'bands' if self.bands > 1 else - 'band', 'between 1 and ' if self.bands > 1 else '', self.bands, self.band4match) + 'band', 'between 1 and ' if self.bands > 1 else '', self.bands, self.band4match + 1) # set footprint_poly given_footprint_poly = CoReg_params['footprint_poly_%s' % ('ref' if imID == 'ref' else 'tgt')] -- GitLab