search

GFZ / arosics

AROSICS - Automated and Robust Open-Source Image Co-Registration Software
https://git.gfz-potsdam.de/danschef/arosics
Apache License 2.0
137 stars 26 forks source link

One of the input images does not have sufficient gray value information (non-no-data values) for placing a matching window at the position (795768.75, 3793443.75). Matching failed. #11

Closed shannonmgbrown closed 1 year ago

shannonmgbrown commented 2 years ago

Description

Hello, I am trying to coregister L5,L7,L8, and S2 images downloaded from google earth engine. I am using Killian Vos' coastsat tool (https://github.com/kvos/CoastSat) to track shoreline change. The site that I am having trouble with is in Wrightsville Beach and in this location, some of the images have UTM zone 17N and some have UTM zone 18M. I also had to make a mask for the pixels with the value -Inf on the border.

The error I am getting is: One of the input images does not have sufficient gray value information (non-no-data values) for placing a matching window at the position (795768.75, 3793443.75). Matching failed.

I pulled the images into arc, and they overlap for most of the images so I am a little lost on this error. I attached them as well as the meta files. Below is my code. Thanks for your help! images.zip

What I Did

site_dir = 'C:\\Users\\RDCHLSMB\\Documents\\CoastSat-master\\data\\Shortterm\\WrightsvilleBeach\\'

im_reference = site_dir+'mask_ref\\im\\2017-11-26-15-47-48_L8_WrightsvilleBeach_pan.tif'

sat = ['L5\\30m\\','L7\\pan\\','L8\\pan\\',] #'\\L5\\30m\\','\\L7\\pan\\','\\L8\\pan\\',
for mission in sat:
    ims = os.listdir(site_dir+mission)
    for images in ims:
        im_target = os.path.join(site_dir,mission,images)
        np.seterr(divide='ignore', invalid='ignore')
        im1 = gdal.Open(im_target)
        driver = gdal.GetDriverByName("GTiff")
        outdata = driver.Create(site_dir+'arosics\\'+mission+images,int(im1.RasterXSize),int(im1.RasterYSize),int(im1.RasterCount),gdal.GDT_Float32)
        outdata.SetGeoTransform(im1.GetGeoTransform())
        outdata.SetProjection(im1.GetProjection())
        for i in range(im1.RasterCount):
            myarray = np.array(im1.GetRasterBand(i+1).ReadAsArray())
            myarray[myarray == -math.inf] = 600
            outdata.GetRasterBand(i+1).WriteArray(myarray)
            outdata.GetRasterBand(i+1).SetNoDataValue(600)
            myarray = []
        outdata = None
        im1 = None
        nodata = 600
        geoArr = gdal.Open(site_dir+'arosics\\'+mission+images)
        tgt_ndarray = geoArr.ReadAsArray()
        tgt_gt = geoArr.GetGeoTransform()
        tgt_prj = geoArr.GetProjection()
        dimen = tgt_ndarray.shape
        if ref_prj == tgt_prj:
            geoArr_reference = GeoArray(ref_ndarray,ref_gt,ref_prj)
            geoArr_target = GeoArray(tgt_ndarray,tgt_gt,tgt_prj)
            kwargs = {'grid_res' : 15, 'window_size' : dimen, 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir': site_dir+'arosics\\'+mission,'q':False,'path_out' : site_dir+mission+images}
            CRL = COREG_LOCAL(GeoArray(ref_ndarray,ref_gt,ref_prj),GeoArray(tgt_ndarray,tgt_gt,tgt_prj), nodata = [600,600],**kwargs)
            CRL.correct_shifts()
        else:
            ref = GeoArray(im_reference)
            tgt = GeoArray(im_target)
            ref.reproject_to_new_grid(prototype=tgt)
            ref_gt = ref.geotransform
            ref_prj = ref.projection
            kwargs = {'grid_res' : 30, 'window_size' : (256,256), 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir' : site_dir+site+'arosics\\'+mission,'q' : False,'path_out' : site_dir+site+mission+images}
            CRL = COREG_LOCAL(GeoArray(ref_ndarray,ref_gt,ref_prj),GeoArray(tgt_ndarray,tgt_gt,tgt_prj), nodata = [600,600],**kwargs)
            CRL.correct_shifts()

This is what arosics printed before the crash:

site_dir = 'C:\\Users\\RDCHLSMB\\Documents\\CoastSat-master\\data\\Shortterm\\WrightsvilleBeach\\'

im_reference = site_dir+'mask_ref\\im\\2017-11-26-15-47-48_L8_WrightsvilleBeach_pan.tif'

sat = ['L5\\30m\\','L7\\pan\\','L8\\pan\\',] #'\\L5\\30m\\','\\L7\\pan\\','\\L8\\pan\\',
for mission in sat:
    ims = os.listdir(site_dir+mission)
    for images in ims:
        im_target = os.path.join(site_dir,mission,images)
        np.seterr(divide='ignore', invalid='ignore')
        im1 = gdal.Open(im_target)
        driver = gdal.GetDriverByName("GTiff")
        outdata = driver.Create(site_dir+'arosics\\'+mission+images,int(im1.RasterXSize),int(im1.RasterYSize),int(im1.RasterCount),gdal.GDT_Float32)
        outdata.SetGeoTransform(im1.GetGeoTransform())
        outdata.SetProjection(im1.GetProjection())
        for i in range(im1.RasterCount):
            myarray = np.array(im1.GetRasterBand(i+1).ReadAsArray())
            myarray[myarray == -math.inf] = 600
            outdata.GetRasterBand(i+1).WriteArray(myarray)
            outdata.GetRasterBand(i+1).SetNoDataValue(600)
            myarray = []
        #outdata.FlushCache()
        outdata = None
        im1 = None
        nodata = 600
        #im_target = 'C:\\Users\\RDCHLSMB\\Documents\\Brown\\CoastSat\\WrightsvilleBeach\\L8\\L8\\pan\\2014-08-14-15-47-42_L8_WrightsvilleBeach_pan.tif'
        geoArr = gdal.Open(site_dir+'arosics\\'+mission+images)
        tgt_ndarray = geoArr.ReadAsArray()
        tgt_gt = geoArr.GetGeoTransform()
        tgt_prj = geoArr.GetProjection()
        dimen = tgt_ndarray.shape
        if ref_prj == tgt_prj:
            geoArr_reference = GeoArray(ref_ndarray,ref_gt,ref_prj)
            geoArr_target = GeoArray(tgt_ndarray,tgt_gt,tgt_prj)
            #CR = COREG(geoArr_reference,geoArr_target,ws = (219, 201))
            #kwargs = {'grid_res' : 15, 'window_size' : (277,411), 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir': os.path.join('C:\\Users\\RDCHLSMB\\Documents\\Brown\\CoastSat\\Arosics',images),'q':False,'path_out' : 'C:\\Users\\RDCHLSMB\\Documents\\Brown\\CoastSat\\Arosics'}
            kwargs = {'grid_res' : 15, 'window_size' : dimen, 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir': site_dir+'arosics\\'+mission,'q':False,'path_out' : site_dir+mission+images}
            CRL = COREG_LOCAL(GeoArray(ref_ndarray,ref_gt,ref_prj),GeoArray(tgt_ndarray,tgt_gt,tgt_prj), nodata = [600,600],**kwargs)
            CRL.correct_shifts()
        else:
            ref = GeoArray(im_reference)
            tgt = GeoArray(im_target)
            ref.reproject_to_new_grid(prototype=tgt)
            ref_gt = ref.geotransform
            ref_prj = ref.projection
            #CR = COREG(geoArr_reference,geoArr_target,ws = (219, 201))
            kwargs = {'grid_res' : 30, 'window_size' : (256,256), 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir' : site_dir+site+'arosics\\'+mission,'q' : False,'path_out' : site_dir+site+mission+images}
            CRL = COREG_LOCAL(GeoArray(ref_ndarray,ref_gt,ref_prj),GeoArray(tgt_ndarray,tgt_gt,tgt_prj), nodata = [600,600],**kwargs)
            CRL.correct_shifts()
site_dir = 'C:\\Users\\RDCHLSMB\\Documents\\CoastSat-master\\data\\Shortterm\\WrightsvilleBeach\\'
​
im_reference = site_dir+'mask_ref\\im\\2017-11-26-15-47-48_L8_WrightsvilleBeach_pan.tif'
​
sat = ['L5\\30m\\','L7\\pan\\','L8\\pan\\',] #'\\L5\\30m\\','\\L7\\pan\\','\\L8\\pan\\',
for mission in sat:
    ims = os.listdir(site_dir+mission)
    for images in ims:
        im_target = os.path.join(site_dir,mission,images)
        np.seterr(divide='ignore', invalid='ignore')
        im1 = gdal.Open(im_target)
        driver = gdal.GetDriverByName("GTiff")
        outdata = driver.Create(site_dir+'arosics\\'+mission+images,int(im1.RasterXSize),int(im1.RasterYSize),int(im1.RasterCount),gdal.GDT_Float32)
        outdata.SetGeoTransform(im1.GetGeoTransform())
        outdata.SetProjection(im1.GetProjection())
        for i in range(im1.RasterCount):
            myarray = np.array(im1.GetRasterBand(i+1).ReadAsArray())
            myarray[myarray == -math.inf] = 600
            outdata.GetRasterBand(i+1).WriteArray(myarray)
            outdata.GetRasterBand(i+1).SetNoDataValue(600)
            myarray = []
        #outdata.FlushCache()
        outdata = None
        im1 = None
        nodata = 600
        #im_target = 'C:\\Users\\RDCHLSMB\\Documents\\Brown\\CoastSat\\WrightsvilleBeach\\L8\\L8\\pan\\2014-08-14-15-47-42_L8_WrightsvilleBeach_pan.tif'
        geoArr = gdal.Open(site_dir+'arosics\\'+mission+images)
        tgt_ndarray = geoArr.ReadAsArray()
        tgt_gt = geoArr.GetGeoTransform()
        tgt_prj = geoArr.GetProjection()
        dimen = tgt_ndarray.shape
        if ref_prj == tgt_prj:
            geoArr_reference = GeoArray(ref_ndarray,ref_gt,ref_prj)
            geoArr_target = GeoArray(tgt_ndarray,tgt_gt,tgt_prj)
            #CR = COREG(geoArr_reference,geoArr_target,ws = (219, 201))
            #kwargs = {'grid_res' : 15, 'window_size' : (277,411), 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir': os.path.join('C:\\Users\\RDCHLSMB\\Documents\\Brown\\CoastSat\\Arosics',images),'q':False,'path_out' : 'C:\\Users\\RDCHLSMB\\Documents\\Brown\\CoastSat\\Arosics'}
            kwargs = {'grid_res' : 15, 'window_size' : dimen, 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir': site_dir+'arosics\\'+mission,'q':False,'path_out' : site_dir+mission+images}
            CRL = COREG_LOCAL(GeoArray(ref_ndarray,ref_gt,ref_prj),GeoArray(tgt_ndarray,tgt_gt,tgt_prj), nodata = [600,600],**kwargs)
            CRL.correct_shifts()
        else:
            ref = GeoArray(im_reference)
            tgt = GeoArray(im_target)
            ref.reproject_to_new_grid(prototype=tgt)
            ref_gt = ref.geotransform
            ref_prj = ref.projection
            #CR = COREG(geoArr_reference,geoArr_target,ws = (219, 201))
            kwargs = {'grid_res' : 30, 'window_size' : (256,256), 'path_out' : 'auto', 'fmt_out' : 'GTIFF', 'projectDir' : site_dir+site+'arosics\\'+mission,'q' : False,'path_out' : site_dir+site+mission+images}
            CRL = COREG_LOCAL(GeoArray(ref_ndarray,ref_gt,ref_prj),GeoArray(tgt_ndarray,tgt_gt,tgt_prj), nodata = [600,600],**kwargs)
            CRL.correct_shifts()
Calculating footprint polygon and actual data corner coordinates for reference image...
Bounding box of calculated footprint for reference image:
    (793687.5, 3787357.5, 797842.5, 3793522.5)
Calculating footprint polygon and actual data corner coordinates for image to be shifted...
Bounding box of calculated footprint for image to be shifted:
    (793695.0, 3793365.0, 799905.0, 3793545.0)
Matching window position (X,Y): 795768.75/3793443.75
C:\Users\RDCHLSMB\Anaconda3\envs\arosics\lib\site-packages\IPython\core\interactiveshell.py:3361: UserWarning: 
First attempt to check the functionality of co-registration failed. Check your input data and parameters. The following error occurred:
  exec(code_obj, self.user_global_ns, self.user_ns)
danschef commented 2 years ago

I think you should first break down all this code to the actual input arguments of the COREG_LOCAL class. Just save the input images as you pass them to AROSICS, open them in a GIS and have a look at the matching window position at 795768.75/3793443.75. I guess this should already make clear why you get the above error message.

If not, upload these images somewhere and post the link here.

danschef commented 1 year ago

I close this due to inactivity. Feel free to re-open if needed.