pyproj Issue

[disbr007]

• arosics version: 0.9.19
• py_tools_ds version: 0.15.0
• Python version: 3.7.4
• Operating System: Windows

I was able to run a global coregistration and write output that looked great. However, when attempting a local coregistration I am running into a pyproj issue, where during the creation of the CoRegPoints_table geodataframe the crs is not accepted by gpd.GeoDataFrame because of a pyproj issue.

I was able to follow the traceback and hard-code the line defining the crs (line 312 in Tie_Point_Grid) to what I know to be the crs I am using and it worked:

    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'
        # (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']
            crs = "epsg:3413"
        else:
            crs = None
        print('CRS: {}'.format(crs))
        GDF = GeoDataFrame(index=range(len(geomPoints)), crs=crs,
                           columns=['geometry', 'POINT_ID', 'X_IM', 'Y_IM', 'X_UTM', 'Y_UTM'])

Here is the traceback when running it without any modifications:

from pathlib import Path

import arosics

prj_dir = Path(r'V:\pgc\data\scratch\jeff\ms\2020aug22')
dems_dir = prj_dir / 'dems' / 'aoi1'
img_dir = prj_dir / 'img' / 'aoi1'

# Imagery is reference -> move everything to it
im_ref = img_dir / 'WV02_20180403212306_103001007A45AA00_18APR03212306-P1BS-502128082080_01_P008_u16ns3413_aoi1.tif'
# Ortho is target -> move it to imagery
im_tar = dems_dir / 'WV02_20180403_103001007A66E000_103001007A45AA00_2m_lsf_seg1_ortho_aoi1.tif'

global_ortho = dems_dir / 'global' / '{}_global.tif'.format(im_tar.stem)
local_ortho = dems_dir / 'local' / '{}_local.tif'.format(im_tar.stem)

# # Global coreg
# cr = arosics.COREG(str(im_ref), str(im_tar), max_shift=100,
#                    path_out=str(global_ortho))
# ss = cr.calculate_spatial_shifts()
# cr.correct_shifts()

# Try local coreg
kwargs = {
    'grid_res': 25,
    'window_size': (250,250),
    'path_out': str(local_ortho),
    'projectDir': str(prj_dir),
    'q': False
}
crl = arosics.COREG_LOCAL(str(im_ref), str(im_tar), **kwargs)
crl.correct_shifts()

Calculating actual data corner coordinates for reference image...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Bounding box of calculated footprint for reference image:
    (-2013501.3361572332, -437212.6544690558, -2005835.4634647095, -433313.0556767457)
Calculating actual data corner coordinates for image to be shifted...
Polygonize progress     |==================================================| 100.0% Complete  => 0:00:00
Bounding box of calculated footprint for image to be shifted:
    (-2013500.0, -437212.0, -2005836.0, -433312.0)
Matching window position (X,Y): -2009717.9995549403/-435251.72088528174
C:\miniconda\envs\aros\lib\site-packages\py_tools_ds\geo\projection.py:51: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
  return CRS.from_user_input(proj).to_proj4().strip()
Initializing tie points grid...
CRS: {'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'utm', 'zone': 0, 'units': 'm', 'no_defs': True}
Traceback (most recent call last):
  File "C:\miniconda\envs\aros\lib\site-packages\IPython\core\interactiveshell.py", line 3331, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-2-a907f5f744d0>", line 32, in <module>
    crl.correct_shifts()
  File "C:\Users\disbr007\AppData\Roaming\Python\Python37\site-packages\arosics\CoReg_local.py", line 605, in correct_shifts
    coreg_info = self.coreg_info
  File "C:\Users\disbr007\AppData\Roaming\Python\Python37\site-packages\arosics\CoReg_local.py", line 577, in coreg_info
    'GCPList': self.tiepoint_grid.GCPList,
  File "C:\Users\disbr007\AppData\Roaming\Python\Python37\site-packages\arosics\CoReg_local.py", line 333, in tiepoint_grid
    self._tiepoint_grid.get_CoRegPoints_table()
  File "C:\Users\disbr007\AppData\Roaming\Python\Python37\site-packages\arosics\Tie_Point_Grid.py", line 320, in get_CoRegPoints_table
    columns=['geometry', 'POINT_ID', 'X_IM', 'Y_IM', 'X_UTM', 'Y_UTM'])
  File "C:\miniconda\envs\aros\lib\site-packages\geopandas\geodataframe.py", line 65, in __init__
    self.crs = crs
  File "C:\miniconda\envs\aros\lib\site-packages\geopandas\geodataframe.py", line 97, in __setattr__
    super(GeoDataFrame, self).__setattr__(attr, val)
  File "C:\miniconda\envs\aros\lib\site-packages\pandas\core\generic.py", line 5287, in __setattr__
    return object.__setattr__(self, name, value)
  File "C:\miniconda\envs\aros\lib\site-packages\geopandas\base.py", line 153, in crs
    self._crs = None if not value else CRS.from_user_input(value)
  File "C:\miniconda\envs\aros\lib\site-packages\pyproj\crs.py", line 580, in from_user_input
    return cls(value)
  File "C:\miniconda\envs\aros\lib\site-packages\pyproj\crs.py", line 436, in __init__
    super(CRS, self).__init__(projstring)
  File "pyproj/_crs.pyx", line 1738, in pyproj._crs._CRS.__init__
pyproj.exceptions.CRSError: Invalid projection: +proj=utm +ellps=WGS84 +datum=WGS84 +zone=0 +units=m +no_defs +type=crs: (Internal Proj Error: proj_create: Error -35: invalid UTM zone number)

I attempted to do an actual fix it but couldn't follow everything I needed to in py_tools_ds.

[danschef]

Thanks for reporting. AROSICS was initially written with full support for UTM and geographic coordinates (EPSG 4326). I just revised the underlying code in py_tools_ds a bit, also to get more compatibility with other projections. However, I guess there will also be some other places in the code that need to be revised to support any kind of projection.

I will look into that as soon as I have some spare time. In the meantime, I recommend to use UTM for your input data.

[danschef]

I created an issue in the GitLab issue tracker regarding the support of other projections: https://gitext.gfz-potsdam.de/danschef/arosics/-/issues/37

[danschef]

AROSICS 1.4.0 adds full support for any kind of projection. See https://git.gfz-potsdam.de/danschef/arosics/-/merge_requests/17.