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.

Answer 1

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()