Merge branch 'enhancement/add_anyproj_compatibility' into 'master'

Provide full support for projections other than UTM or geographic coordinates.

Closes #37

See merge request !17

arosics/CoReg.py

@@ -537,13 +537,10 @@ class COREG(object):
         if not prj_equal(self.ref.prj, self.shift.prj):
             from pyproj import CRS
 
-            def get_prjdesc(proj):
-                crs = CRS.from_user_input(proj)
-                return "%s (EPSG: %d)" % (crs.name, crs.to_epsg())
-
             raise RuntimeError(
                 'Input projections are not equal. Different projections are currently not supported. '
-                'Got %s / %s.' % (get_prjdesc(self.ref.prj), get_prjdesc(self.shift.prj)))
+                'Got %s vs. %s.' % (CRS.from_user_input(self.ref.prj).name,
+                                    CRS.from_user_input(self.shift.prj).name))
 
     def _get_overlap_properties(self) -> None:
         overlap_tmp = get_overlap_polygon(self.ref.poly, self.shift.poly, self.v)
@@ -594,10 +591,10 @@ class COREG(object):
         import folium
         import geojson
 
-        refPoly = reproject_shapelyGeometry(self.ref.poly, self.ref.epsg, 4326)
-        shiftPoly = reproject_shapelyGeometry(self.shift.poly, self.shift.epsg, 4326)
-        overlapPoly = reproject_shapelyGeometry(self.overlap_poly, self.shift.epsg, 4326)
-        matchBoxPoly = reproject_shapelyGeometry(self.matchBox.mapPoly, self.shift.epsg, 4326)
+        refPoly = reproject_shapelyGeometry(self.ref.poly, self.ref.prj, 4326)
+        shiftPoly = reproject_shapelyGeometry(self.shift.poly, self.shift.prj, 4326)
+        overlapPoly = reproject_shapelyGeometry(self.overlap_poly, self.shift.prj, 4326)
+        matchBoxPoly = reproject_shapelyGeometry(self.matchBox.mapPoly, self.shift.prj, 4326)
 
         m = folium.Map(location=tuple(np.array(overlapPoly.centroid.coords.xy).flatten())[::-1])
         for poly in [refPoly, shiftPoly, overlapPoly, matchBoxPoly]:

arosics/CoReg_local.py

@@ -408,7 +408,7 @@ class COREG_LOCAL(object):
     def CoRegPoints_table(self) -> GeoDataFrame:
         """Return a GeoDataFrame containing all the results from coregistration for all points in the tie point grid.
 
-        Columns of the GeoDataFrame: 'geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM','X_WIN_SIZE', 'Y_WIN_SIZE',
+        Columns of the GeoDataFrame: 'geometry','POINT_ID','X_IM','Y_IM','X_MAP','Y_MAP','X_WIN_SIZE', 'Y_WIN_SIZE',
                                      'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M', 'Y_SHIFT_M', 'ABS_SHIFT' and 'ANGLE'
         """
         return self.tiepoint_grid.CoRegPoints_table
@@ -692,7 +692,7 @@ class COREG_LOCAL(object):
         folium.GeoJson(points_values).add_to(map_osm)
 
         # add overlap polygon
-        overlapPoly = reproject_shapelyGeometry(self.COREG_obj.overlap_poly, self.im2shift.epsg, 4326)
+        overlapPoly = reproject_shapelyGeometry(self.COREG_obj.overlap_poly, self.im2shift.prj, 4326)
         gjs = geojson.Feature(geometry=overlapPoly, properties={})
         folium.GeoJson(gjs).add_to(map_osm)
 

arosics/Tie_Point_Grid.py

@@ -37,7 +37,7 @@ from shapely.geometry import Point
 
 # internal modules
 from .CoReg import COREG
-from py_tools_ds.geo.projection import isProjectedOrGeographic, isLocal, get_UTMzone
+from py_tools_ds.geo.projection import isLocal
 from py_tools_ds.io.pathgen import get_generic_outpath
 from py_tools_ds.processing.progress_mon import ProgressBar
 from py_tools_ds.geo.vector.conversion import points_to_raster
@@ -191,7 +191,7 @@ class Tie_Point_Grid(object):
     def CoRegPoints_table(self):
         """Return a GeoDataFrame containing all the results from coregistration for all points in the tie point grid.
 
-        Columns of the GeoDataFrame: 'geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM','X_WIN_SIZE', 'Y_WIN_SIZE',
+        Columns of the GeoDataFrame: 'geometry','POINT_ID','X_IM','Y_IM','X_MAP','Y_MAP','X_WIN_SIZE', 'Y_WIN_SIZE',
                                      'X_SHIFT_PX','Y_SHIFT_PX', 'X_SHIFT_M', 'Y_SHIFT_M', 'ABS_SHIFT' and 'ANGLE'
         """
         if self._CoRegPoints_table is not None:
@@ -249,7 +249,7 @@ class Tie_Point_Grid(object):
     def _exclude_bad_XYpos(self, GDF):
         """Exclude all points outside of the image overlap area and where the bad data mask is True (if given).
 
-        :param GDF:     <geopandas.GeoDataFrame> must include the columns 'X_UTM' and 'Y_UTM'
+        :param GDF:     <geopandas.GeoDataFrame> must include the columns 'X_MAP' and 'Y_MAP'
         :return:
         """
         from skimage.measure import points_in_poly  # import here to avoid static TLS ImportError
@@ -264,7 +264,7 @@ class Tie_Point_Grid(object):
 
         # exclude all points where bad data mask is True (e.g. points on clouds etc.)
         orig_len_GDF = len(GDF)  # length of GDF after dropping all points outside the overlap polygon
-        mapXY = np.array(GDF.loc[:, ['X_UTM', 'Y_UTM']])
+        mapXY = np.array(GDF.loc[:, ['X_MAP', 'Y_MAP']])
         GDF['REF_BADDATA'] = self.COREG_obj.ref.mask_baddata.read_pointData(mapXY) \
             if self.COREG_obj.ref.mask_baddata is not None else False
         GDF['TGT_BADDATA'] = self.COREG_obj.shift.mask_baddata.read_pointData(mapXY) \
@@ -330,30 +330,20 @@ class Tie_Point_Grid(object):
     def get_CoRegPoints_table(self):
         assert self.XY_points is not None and self.XY_mapPoints is not None
 
-        # create a dataframe containing 'geometry','POINT_ID','X_IM','Y_IM','X_UTM','Y_UTM'
+        # create a dataframe containing 'geometry','POINT_ID','X_IM','Y_IM','X_MAP','Y_MAP'
         # (convert imCoords to mapCoords
         XYarr2PointGeom = np.vectorize(lambda X, Y: Point(X, Y), otypes=[Point])
         geomPoints = np.array(XYarr2PointGeom(self.XY_mapPoints[:, 0], self.XY_mapPoints[:, 1]))
 
-        if isLocal(self.COREG_obj.shift.prj):
-            crs = None
-        elif isProjectedOrGeographic(self.COREG_obj.shift.prj) == 'geographic':
-            crs = dict(ellps='WGS84', datum='WGS84', proj='longlat')
-        elif isProjectedOrGeographic(self.COREG_obj.shift.prj) == 'projected':
-            UTMzone = abs(get_UTMzone(prj=self.COREG_obj.shift.prj))
-            south = get_UTMzone(prj=self.COREG_obj.shift.prj) < 0
-            crs = dict(ellps='WGS84', datum='WGS84', proj='utm', zone=UTMzone, south=south, units='m', no_defs=True)
-            if not south:
-                del crs['south']
-        else:
-            crs = None
+        crs = self.COREG_obj.shift.prj if not isLocal(self.COREG_obj.shift.prj) else None
 
-        GDF = GeoDataFrame(index=range(len(geomPoints)), crs=crs,
-                           columns=['geometry', 'POINT_ID', 'X_IM', 'Y_IM', 'X_UTM', 'Y_UTM'])
+        GDF = GeoDataFrame(index=range(len(geomPoints)),
+                           crs=crs,
+                           columns=['geometry', 'POINT_ID', 'X_IM', 'Y_IM', 'X_MAP', 'Y_MAP'])
         GDF['geometry'] = geomPoints
         GDF['POINT_ID'] = range(len(geomPoints))
         GDF.loc[:, ['X_IM', 'Y_IM']] = self.XY_points
-        GDF.loc[:, ['X_UTM', 'Y_UTM']] = self.XY_mapPoints
+        GDF.loc[:, ['X_MAP', 'Y_MAP']] = self.XY_mapPoints
 
         # exclude offsite points and points on bad data mask
         GDF = self._exclude_bad_XYpos(GDF)
@@ -649,10 +639,10 @@ class Tie_Point_Grid(object):
                               'out of the %s available tie points.' % avail_TP)
 
             # calculate GCPs
-            GDF['X_UTM_new'] = GDF.X_UTM + GDF.X_SHIFT_M
-            GDF['Y_UTM_new'] = GDF.Y_UTM + GDF.Y_SHIFT_M
-            GDF['GCP'] = GDF.apply(lambda GDF_row: gdal.GCP(GDF_row.X_UTM_new,
-                                                            GDF_row.Y_UTM_new,
+            GDF['X_MAP_new'] = GDF.X_MAP + GDF.X_SHIFT_M
+            GDF['Y_MAP_new'] = GDF.Y_MAP + GDF.Y_SHIFT_M
+            GDF['GCP'] = GDF.apply(lambda GDF_row: gdal.GCP(GDF_row.X_MAP_new,
+                                                            GDF_row.Y_MAP_new,
                                                             0,
                                                             GDF_row.X_IM,
                                                             GDF_row.Y_IM),
@@ -827,7 +817,7 @@ class Tie_Point_Grid(object):
         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_UTM'], GDF2pass['Y_UTM'], GDF2pass[attrName]
+        X_coords, Y_coords, ABS_SHIFT = GDF2pass['X_MAP'], GDF2pass['Y_MAP'], GDF2pass[attrName]
 
         xmin, ymin, xmax, ymax = GDF2pass.total_bounds
 
@@ -997,7 +987,7 @@ class Tie_Point_Refiner(object):
         """Detect geometric outliers between point cloud of source and estimated coordinates using RANSAC algorithm."""
         # from skimage.transform import PolynomialTransform  # import here to avoid static TLS ImportError
 
-        src_coords = np.array(inGDF[['X_UTM', 'Y_UTM']])
+        src_coords = np.array(inGDF[['X_MAP', 'Y_MAP']])
         xyShift = np.array(inGDF[['X_SHIFT_M', 'Y_SHIFT_M']])
         est_coords = src_coords + xyShift
 

docs/usage/input_data_requirements.rst

@@ -25,17 +25,21 @@ header file.
 Supported projections
 ~~~~~~~~~~~~~~~~~~~~~
 
-AROSICS provides full support for UTM projections and geographic coordinates. Providing support for other projections
-is currently work in progress (see `here <https://git.gfz-potsdam.de/danschef/arosics/-/issues/37>`__ for the
-status) and may lead to some incompatibities.
+AROSICS was initially written with support for UTM and geographic coordinates only. Full support for any other
+projection was added in version 1.4.0. However, make sure your input images have the same projection. Different
+projections for the reference and target image are currently not supported.
+
+AROSICS can also be applied to images without any projection and geocoding information. In this case, however,
+the input images need to have the same pixel size and must cover more or less the same spatial area
+(with a shift a few pixels at most).
 
 
 Geographic overlap
 ~~~~~~~~~~~~~~~~~~
 
-The input images must have a geographic overlap but clipping them to same geographical extent is NOT neccessary.
-The image overlap area is automatically calculated. Thereby, no-data regions within the images are automatically
-respected.
+The input images must have a geographic overlap but clipping them to same geographical extent is NOT neccessary
+The image overlap area is automatically calculated, given that your input images have valid geocoding and projection
+information). Thereby, no-data regions within the images are automatically respected.
 
 
 Spatial resolution

requirements.txt

 cartopy
 cmocean
 gdal
-geoarray>=0.9.0
+geoarray>=0.11.0
 geojson
 folium>=0.6.0,!=0.12.0
 geopandas

setup.py

@@ -43,7 +43,7 @@ req = [
     'folium>=0.6.0,!=0.12.0',
     'gdal',
     'geojson',
-    'geoarray>=0.9.0',
+    'geoarray>=0.11.0',
     'geopandas',
     'matplotlib',
     'numpy',

tests/CI_docker/context/environment_arosics.yml

@@ -9,7 +9,7 @@ dependencies:
   - cartopy
   - cmocean
   - gdal
-  - geoarray>=0.9.0
+  - geoarray>=0.11.0
   - geopandas
   - holoviews
   - ipython  # needed to test interactive plotting

tests/test_COREG_LOCAL.py

@@ -119,6 +119,61 @@ class CompleteWorkflow_INTER1_S2A_S2A(unittest.TestCase):
         CRL = COREG_LOCAL(ref, tgt, **dict(CPUs=32,
                                            **self.coreg_kwargs))
         CRL.calculate_spatial_shifts()
+        # CRL.view_CoRegPoints()
+
+    def test_calculation_of_tie_point_grid_noepsg(self):
+        """Test local coregistration with a proj. other than LonLat and UTM and a WKT which has no EPSG code (FORCE)."""
+        wkt_noepsg = \
+            """
+            PROJCRS["BU MEaSUREs Lambert Azimuthal Equal Area - SA - V01",
+                BASEGEOGCRS["WGS 84",
+                    DATUM["World Geodetic System 1984",
+                        ELLIPSOID["WGS 84",6378137,298.257223563,
+                            LENGTHUNIT["metre",1]]],
+                    PRIMEM["Greenwich",0,
+                        ANGLEUNIT["degree",0.0174532925199433]],
+                    ID["EPSG",4326]],
+                CONVERSION["unnamed",
+                    METHOD["Lambert Azimuthal Equal Area",
+                        ID["EPSG",9820]],
+                    PARAMETER["Latitude of natural origin",-15,
+                        ANGLEUNIT["degree",0.0174532925199433],
+                        ID["EPSG",8801]],
+                    PARAMETER["Longitude of natural origin",-60,
+                        ANGLEUNIT["degree",0.0174532925199433],
+                        ID["EPSG",8802]],
+                    PARAMETER["False easting",0,
+                        LENGTHUNIT["metre",1],
+                        ID["EPSG",8806]],
+                    PARAMETER["False northing",0,
+                        LENGTHUNIT["metre",1],
+                        ID["EPSG",8807]]],
+                CS[Cartesian,2],
+                    AXIS["easting",east,
+                        ORDER[1],
+                        LENGTHUNIT["metre",1]],
+                    AXIS["northing",north,
+                        ORDER[2],
+                        LENGTHUNIT["metre",1]]]
+            """
+        wkt_noepsg = ' '.join(wkt_noepsg.split())
+
+        # overwrite prj
+        ref = GeoArray(self.ref_path)
+        ref.to_mem()
+        ref.filePath = None
+        tgt = GeoArray(self.tgt_path)
+        tgt.to_mem()
+        tgt.filePath = None
+
+        ref.prj = wkt_noepsg
+        tgt.prj = wkt_noepsg
+
+        # get instance of COREG_LOCAL object
+        CRL = COREG_LOCAL(ref, tgt, **dict(CPUs=32,
+                                           **self.coreg_kwargs))
+        CRL.calculate_spatial_shifts()
+        # CRL.view_CoRegPoints()
 
 
 if __name__ == '__main__':