CRL.view_CoRegPoints get a SRS issue

[hulius7]

• arosics version: 1.9.2
• Python version: 3.9.16
• Operating System: Windows 10 64bit

Description

I managed to successfully register a 5m panchromatic image on a 10m sentinel-2 band, using global and local registration functions. Everything worked perfectly until I tried to plot the CoRegPoints.

I got this error: RuntimeError: Translating source or target SRS failed: EPSG:32632

This is the code:

CRL.view_CoRegPoints(hide_filtered=False, shapes2plot='vectors', figsize=(20, 20))

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[60], line 1
----> 1 CRL.view_CoRegPoints(hide_filtered=False, shapes2plot='vectors', figsize=(20, 20))

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\arosics\CoReg_local.py:566, in COREG_LOCAL.view_CoRegPoints(self, shapes2plot, attribute2plot, cmap, exclude_fillVals, backgroundIm, hide_filtered, figsize, figsize_multiplier, title, vector_scale, savefigPath, savefigDPI, showFig, vmin, vmax, return_map)
    563     figsize = (figsize[0] * figsize_multiplier, figsize[1] * figsize_multiplier)
    565 # get a map showing the background image
--> 566 fig, ax = backgroundIm.show_map(figsize=figsize,
    567                                 nodataVal=self.nodata[1],
    568                                 return_map=True,
    569                                 band=self.COREG_obj.shift.band4match)
    571 # set figure title
    572 dict_attr_title = dict(
    573     X_WIN_SIZE='size of the matching window in x-direction [pixels]',
    574     Y_WIN_SIZE='size of the matching window in y-direction [pixels]',
   (...)
    584     RELIABILITY='reliability of the computed shift vector'
    585 )

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\geoarray\baseclasses.py:1455, in GeoArray.show_map(self, xlim, ylim, band, boundsMap, boundsMapPrj, out_epsg, figsize, interpolation, vmin, vmax, pmin, pmax, cmap, nodataVal, res_factor, return_map)
   1452 # get image to plot
   1453 # (reproject to LonLat as workaround in case self.epsg is None because cartopy relies on an existing EPSG code)
   1454 nodataVal = nodataVal if nodataVal is not None else self.nodata
-> 1455 gA2plot = GeoArray(*self._get_plottable_image(xlim, ylim, band,
   1456                                               boundsMap=boundsMap,
   1457                                               boundsMapPrj=boundsMapPrj,
   1458                                               res_factor=res_factor,
   1459                                               nodataVal=nodataVal,
   1460                                               # FIXME EPSG:4326 fails for extraterrestrial data
   1461                                               out_prj=self.epsg or 4326
   1462                                               ),
   1463                    nodata=nodataVal)
   1464 image2plot = gA2plot[:]
   1466 # create map

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\geoarray\baseclasses.py:1235, in GeoArray._get_plottable_image(self, xlim, ylim, band, boundsMap, boundsMapPrj, res_factor, nodataVal, out_prj, ignore_rotation)
   1233 if xdim or ydim or out_prj:
   1234     from py_tools_ds.geo.raster.reproject import warp_ndarray
-> 1235     image2plot, gt, prj = warp_ndarray(image2plot, self.geotransform, self.projection,
   1236                                        out_XYdims=(xdim, ydim),
   1237                                        in_nodata=in_nodata,
   1238                                        out_nodata=out_nodata,
   1239                                        transformerOptions=transOpt,
   1240                                        out_prj=out_prj,
   1241                                        q=True)
   1242     if transOpt and 'NO_GEOTRANSFORM' in ','.join(transOpt):
   1243         image2plot = np.flipud(image2plot)

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\py_tools_ds\geo\raster\reproject.py:456, in warp_ndarray(ndarray, in_gt, in_prj, out_prj, out_dtype, out_gsd, out_bounds, out_bounds_prj, out_XYdims, rspAlg, in_nodata, out_nodata, in_alpha, out_alpha, targetAlignedPixels, gcpList, polynomialOrder, options, transformerOptions, warpOptions, CPUs, warpMemoryLimit, progress, q)
    447     in_ds = gdal.Translate(
    448         '', in_ds, format='MEM',
    449         outputSRS=get_SRS(out_prj),
    450         GCPs=gcpList,
    451         callback=ProgressBar(prefix='Translating progress', timeout=None) if progress and not q else None
    452     )
    453     # NOTE: options = ['SPARSE_OK=YES'] ## => what is that for?
    454 
    455 # GDAL Warp
--> 456 res_ds = gdal.Warp(
    457     '', in_ds, format='MEM',
    458     dstSRS=get_SRS(out_prj),
    459     outputType=get_GDT(out_dtype) if out_dtype else get_GDT(in_dtype_np),
    460     xRes=out_gsd[0],
    461     yRes=out_gsd[1],
    462     outputBounds=out_bounds,
    463     outputBoundsSRS=get_SRS(out_bounds_prj),
    464     width=out_XYdims[0],
    465     height=out_XYdims[1],
    466     resampleAlg=rspAlg,
    467     srcNodata=in_nodata,
    468     dstNodata=out_nodata,
    469     srcAlpha=in_alpha,
    470     dstAlpha=out_alpha,
    471     options=options if options else [],
    472     warpOptions=warpOptions or [],
    473     transformerOptions=transformerOptions or [],
    474     targetAlignedPixels=targetAlignedPixels,
    475     tps=True if gcpList else False,
    476     polynomialOrder=polynomialOrder,
    477     warpMemoryLimit=warpMemoryLimit,
    478     callback=ProgressBar(prefix='Warping progress    ', timeout=None) if progress and not q else None,
    479     callback_data=[0],
    480     errorThreshold=0.125,  # this is needed to get exactly the same output like the console version of GDAL warp
    481 )
    483 gdal.SetConfigOption('GDAL_NUM_THREADS', None)
    485 if res_ds is None:

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\osgeo\gdal.py:1032, in Warp(destNameOrDestDS, srcDSOrSrcDSTab, **kwargs)
   1029 _WarnIfUserHasNotSpecifiedIfUsingExceptions()
   1031 if 'options' not in kwargs or isinstance(kwargs['options'], (list, str)):
-> 1032     (opts, callback, callback_data) = WarpOptions(**kwargs)
   1033 else:
   1034     (opts, callback, callback_data) = kwargs['options']

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\osgeo\gdal.py:1012, in WarpOptions(options, format, srcBands, dstBands, outputBounds, outputBoundsSRS, xRes, yRes, targetAlignedPixels, width, height, srcSRS, dstSRS, coordinateOperation, srcAlpha, dstAlpha, warpOptions, errorThreshold, warpMemoryLimit, creationOptions, outputType, workingType, resampleAlg, srcNodata, dstNodata, multithread, tps, rpc, geoloc, polynomialOrder, transformerOptions, cutlineDSName, cutlineLayer, cutlineWhere, cutlineSQL, cutlineBlend, cropToCutline, copyMetadata, metadataConflictValue, setColorInterpretation, overviewLevel, callback, callback_data)
   1009 if return_option_list:
   1010     return new_options
-> 1012 return (GDALWarpAppOptions(new_options), callback, callback_data)

File ~\AppData\Local\miniconda3\envs\arosics\Lib\site-packages\osgeo\gdal.py:5406, in GDALWarpAppOptions.__init__(self, *args)
   5404 def __init__(self, *args):
   5405     r"""__init__(GDALWarpAppOptions self, char ** options) -> GDALWarpAppOptions"""
-> 5406     _gdal.GDALWarpAppOptions_swiginit(self, _gdal.new_GDALWarpAppOptions(*args))

RuntimeError: Translating source or target SRS failed:
EPSG:32632

[danschef]

Sorry for the late reply, now I have some time left to look into this. However, I guess I need some sample data to reproduce it. Could you upload it somewhere?

[hulius7]

Hi, thank you for the reply! I made some changes by switching to point mode and specifying the background as "tgt". However, apparently, everything is working just fine now, I just rerun the old code. It was most likely my mistake. Thank you very much for this amazing tool"

[danschef]

Ok, sounds like the issue is solved for now, so I close this here. Feel free to reopen if it comes up again.