Error when adding the -g option

[whyjz]

Hi,

I have been trying testautoRIFT.py and it runs ok with the most basic command line arguments:

$ python testautoRIFT.py -m ../../Landsat8/LC08_L1TP_062017_20180818_20180829_01_T1_B8_s.TIF -s ../../Landsat8/LC08_L1TP_062017_20180903_20180912_01_T1_B8_s.TIF -fo 1

The command will generate a .mat file but what I want is a geotiff file that I can load on my GIS software. I took a look at the documentation and used testGeogridOptical.py to generate some files, but when I rerun testautoRIFT.py using the following arguments, it throws out an error:

$ python testautoRIFT.py -m ../../Landsat8/LC08_L1TP_062017_20180818_20180829_01_T1_B8_s.TIF -s ../../Landsat8/LC08_L1TP_062017_20180903_20180912_01_T1_B8_s.TIF -fo 1 -g window_location.tif

Traceback (most recent call last): File "testautoRIFT.py", line 378, in srs = ds.GetSpatialRef() File "/data/whyj/Software/anaconda3/envs/autorift/lib/python3.7/site-packages/osgeo/gdal.py", line 2014, in getattr = lambda self, name: _swig_getattr(self, Dataset, name) File "/data/whyj/Software/anaconda3/envs/autorift/lib/python3.7/site-packages/osgeo/gdal.py", line 74, in _swig_getattr return _swig_getattr_nondynamic(self, class_type, name, 0) File "/data/whyj/Software/anaconda3/envs/autorift/lib/python3.7/site-packages/osgeo/gdal.py", line 69, in _swig_getattr_nondynamic return object.getattr(self, name) AttributeError: type object 'object' has no attribute 'getattr'

Any ideas for fixing this error?

[leiyangleon]

Hello, thanks for testing the software and your valuable feedback. There is no problem regarding your command line use. The error basically implies the -g input "window_location.tif" does not have a valid spatial reference system (srs) attribute. So the problem might be in that file. To double check, make sure you

1. open "window_location.tif" in GIS software and click on one or few pixels for the values (x- and y-index in the Landsat images) to check if they are consistent.
2. check the srs of that file by typing "gdalsrsinfo window_location.tif"
3. also "window_location.tif" is from testGeogridOptical.py by using a DEM, which must be in northing/easting (unit of m; not lat/lon). If your DEM is in lat/lon, use gdal to transform it to polar stereographic or UTM

If all the above checks are done and the error persists, please come back and post the results. We will go from there.

[whyjz]

Thanks for the message! I followed your suggestions and here are the test results:

1. yes, there are two bands in window_location.tif and they correspond to the x-index and y-index respectively
2. here's the result: $ gdalsrsinfo window_location.tif

PROJ.4 : +proj=utm +zone=7 +datum=WGS84 +units=m +no_defs

OGC WKT : PROJCS["WGS 84 / UTM zone 7N", GEOGCS["WGS 84", DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree",0.0174532925199433, AUTHORITY["EPSG","9122"]], AUTHORITY["EPSG","4326"]], PROJECTION["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",-141], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1, AUTHORITY["EPSG","9001"]], AXIS["Easting",EAST], AXIS["Northing",NORTH], AUTHORITY["EPSG","32607"]]

The SRS agrees with the projection my input geotiff files, which is EPSG:32607 (UTM 7N)

3. The DEM surely uses the same SRS as well because I created it using $ gdal_calc.py -A ../../Landsat8/LC08_L1TP_062017_20180818_20180829_01_T1_B8_s.TIF --calc=1 --outfile=fakedem.tif which is a fake DEM with all pixels of 1. I am guessing this might be why the script doesn't work, but can't figure out how. I thought it should work if I just want to get the pixel offset in x and y axes in geotiff format.

Thanks!!

[leiyangleon]

Hmm, that is weird. All the inputs and command lines seem okay. For me to re-run autoRIFT on your files, can you put the input files (2 Landsat images + 1 window_location) on Google drive and share the link to me?

[leiyangleon]

I managed to run autoRIFT thru without errors using the same command line as yours. See attached figure for displacement in x direction (blue is -5 and red is +5).

The only problem I saw is your window_locations.tif file is too big that took 2 hrs to finish, which is from the fine DEM grid with Landsat8 spacing (15m). Usually, to do feature tracking, you want a sparse grid with the grid spacing >= pixel_size smallest_window_size. In that way, the local windows for normalized cross-correlation (NCC) will not overlap thus assuring independent displacement results; otherwise, the displacement results as well as the noise are spatially correlated and autoRIFT's NDC filter won't work to remove the noise (see the large blue and red noise in the attached figure). FYI, in your example, pixel_size=15m, smallest_window_size=16 (adjustable in testautoRIFT.py), and grid_spacing is also 15m, which violates the above relation. Try to downsample your DEM at least to a spacing of pixel_size smallest_window_size=15*16=240m or you can adjust the smallest_window_size (see the GitHub manual for variables ChipSizeMinX, ChipSizeMaxX and ChipSize0X, and also line 206-208 in testautoRIFT.py to adjust them) to avoid overlapping of the windows.

As for your error message, I have never seen it before. Since it is related to gdal, my suggestions would be to reinstall the updated version of gdal. You can verify the gdal install by opening a geotiff file using gdal and importing the spatial reference system using the Python commands like ds = gdal.Open($Your GeoTIFF File Name$) srs = ds.GetSpatialRef() If the error persists, you can post that problem in the gdal user forum as autoRIFT assumes and requests a successful install of gdal.

[whyjz]

Thanks for the update. I was not aware that the DEM has something to do with pixel skip, but it makes sense to me.

Regarding the gdal issue, I just sort of solved it myself based on this page https://gis.stackexchange.com/questions/60371/gdal-python-how-do-i-get-coordinate-system-name-from-spatialreference I thought GetSpatialRef() is an ogr function that works for shapefiles instead of geotiffs? I tested srs = ds.GetSpatialRef() with gdal 2.2 + python 3.7 and gdal 2.0 + python 3.5 but have had no luck. But maybe I am wrong -- anyway my solution now is to replace line 378 in testautoRIFT.py with srs = osr.SpatialReference(wkt=proj) and then the script seems to work well. Just wanted to post it here for your future reference.

Thanks again for your help!