Affine transformation of with Python

Question

Translating my code from MATLAB I am trying to apply an affine transformation to georeference a TERRA-ASTER satellite image using Python 3,

import numpy as np
import matplotlib.pyplot as plt
from skimage import transform as tf
from skimage import io, exposure

naivasha_rgb = io.imread('naivasha.tif')

using a tranformation matrix from pairs of coordinates (src,dst) of the four scene corners.

# upper left
# lower left
# upper right
# lower right
movingpoints = np.array([
    [36.214332,-0.319922],
    [36.096003,-0.878267],
    [36.770406,-0.400443],
    [36.652213,-0.958743]])
fixedpoints = np.array([
    [0,0],
    [0,4200],
    [4100,0],
    [4100,4200]])

using functions from the Skimage package

tform = tf.estimate_transform('affine',movingpoints,fixedpoints)
newnaivasha_rgb = tf.warp(naivasha_rgb,tform.inverse)

The tform is the same as in MATLAB but transposed but the warping creates a plane green image, not the one expected as in MATLAB (see WeTransfer link below).

newnaivasha_rgb = exposure.rescale_intensity(newnaivasha_rgb,
    out_range='uint8')
    
plt.figure()
io.imshow(newnaivasha_rgb)
plt.axis('off')
plt.show()

Any ideas? The numerical values indicate a low-contrast or no image, depending whether I use tform or tform.inverse with tf.warp, according to solutions I googled.

Download of images from WeTransfer (50 MB) including input image naivasha.tif and output images naivasha_georef_python.png and naivasha_georef_matlab.jpg.

Not that I have much knowledge of the `skimage` package, but I'm wondering if your `moving_points` need to be relative to the `np.array` indices of `naivasha_rgb` rather than lat/lon. — rickhg12hs, Dec 14 '21 at 11:46
Possibly, at least I then get an image – thank you! So the Python functions seem to work differently from those I used with MATLAB. — Martin H. Trauth, Dec 14 '21 at 12:43
I don't know if `skimage` would use them, but that TIFF doesn't seem to have any GeoTIFF metadata. — rickhg12hs, Dec 14 '21 at 12:57
That's right, but the TIFF isn't a GeoTIFF. It has been processed and saved from the original TERRA-ASTER HDF files, which was quite a challenge already using Python, as compared with the simple MATLAB solution. http://mres.uni-potsdam.de/index.php/2017/02/23/importing-and-georeferencing-terra-aster-images-with-matlab/ — Martin H. Trauth, Dec 14 '21 at 13:26
I only mentioned the GeoTIFF because Python has no means of determining lat/lon from `naivasha.tif`. If Matlab is more convenient, why not just use that? — rickhg12hs, Dec 14 '21 at 13:56
Well, good question :o) I am using MATLAB since ~30 years now, wrote three MATLAB-based textbooks on data analysis in geosciences, but young people seem to prefer Python. Therefore I started teaching a course using MATLAB, Python and Julia (probably replacing Python in the next years) and got fascinated by the idea to translate the 5th edition of my most popular book into Python. I will certainly never use Python in my research. http://mres.uni-potsdam.de — Martin H. Trauth, Dec 14 '21 at 14:08
There are perhaps better Python packages for GIS applications (GDAL?). Maybe [GIS Stack Exchange](https://gis.stackexchange.com/) or [Earth Science Stack Exchange](https://earthscience.stackexchange.com/) will have some Python package recommendations. Yeah, Matlab is expensive, and Julia is impressive but has a so very long way to go to catch up to Python's package repertoire. — rickhg12hs, Dec 14 '21 at 14:29
BTW, Python+Numpy+Numba can go at Julia speed with very little effort. Both Numba and Julia use LLVM. Julia's syntax/grammar is IMHO much better than Python's but I miss the breadth and depth of the Python packages. — rickhg12hs, Dec 14 '21 at 14:39
Yes, true, but the documentation of these packages is really bad. They are pretty cryptic and the examples, if any, far too complex to be used for someone with limited time. I am already using GDAL but maybe I will follow your advice and dive a little deeper into it. It is really difficult to find out about the quality of packages. Someone has warned me not to use Pandas. Today I found a promissing package for spectral data but this seems to have too many bugs. — Martin H. Trauth, Dec 14 '21 at 16:39
Yes, the huge number of packages for Python (as compared with Julia) is a strength of Python but also one of the disadvantages over MATLAB, having >90 toolboxes maintained by professionals at a company. — Martin H. Trauth, Dec 14 '21 at 16:40
Yes, quality can be an issue, but can usually be nearly verified in the scope of your particular application by researching previous use or trying your own typical and corner cases with known solutions. Don't use Pandas? Wow, never heard of that before! — rickhg12hs, Dec 14 '21 at 19:05

score 2 · Answer 1 · answered Dec 15 '21 at 16:45

Thanks to the advice of rickhg12hs I was able to solve the problem now. I converted the longitude and latitude values into pixel coordinates with respect to the image and then the affine transformation yields the expected result. Thanks rickhg12hs!

import numpy as np
import matplotlib.pyplot as plt
from skimage import transform as tf
from skimage import io

#%%
naivasha_rgb = io.imread('naivasha.tif')

#%%
plt.figure()
io.imshow(naivasha_rgb)
plt.axis('off')
plt.show()

#%%
UL = [-0.319922,36.214332]
LL = [-0.878267,36.096003]
UR = [-0.400443,36.770406]
LR = [-0.958743,36.652213]

#%%
# 4200 rows correspond to the latitudes  [1]
# 4100 cols correspond to the longitudes [0]
movingpoints = np.array([
    [0.,4100.*(UL[1]-LL[1])/(UR[1]-LL[1])],
    [4200.*(UL[0]-LL[0])/(UL[0]-LR[0]),0.],
    [4200.*(UL[0]-UR[0])/(UL[0]-LR[0]),4100.],
    [4200.,4100.*(LR[1]-LL[1])/(UR[1]-LL[1])]])

#%%
fixedpoints = np.array([
    [0.,0.],
    [4200.,0.],
    [0.,4100.],
    [4200.,4100.]]) 

#%% 
fixedpoints = np.fliplr(fixedpoints)
movingpoints = np.fliplr(movingpoints)
    
#%%
tform = tf.estimate_transform('affine',movingpoints,fixedpoints)

#%%
newnaivasha_rgb = tf.warp(naivasha_rgb,tform)

#%%
plt.figure()
plt.imshow(np.flipud(newnaivasha_rgb),
    extent=[36.096003,36.770406,-0.319922,-0.958743])
ax=plt.gca()
ax.invert_yaxis()
plt.savefig('newnaivasha_rgb.png')
plt.show()

Affine transformation of with Python

1 Answers1