<jupyter_start><jupyter_text>Stitching thermal drone imagery<jupyter_code># Variable definition
imgDir = 'data/lsm/Flight 1 - Library/thermal-data/'
pathFile = 'data/lsm/Flight 1.csv'
utcDiff = -1
pathColumns = ['timestamp','latitude','longitude','altitude(m)','pitch(deg)','roll(deg)','yaw(deg)']
minAlt, maxAlt, altCol = 119.7, 120.3, 'altitude(m)'
imgTimes = [1619811011923,1619811099149] # Specify times if wanted
utmZone, hemisphere = '31n', 'north'
pxSize = 0.27 # pixel size (m) - for now
distFilt = False
resolutions = [2, 6, 60]
# Import packages
#%pip install opencv-python
import glob
import rasterio as ro
import re
from datetime import datetime
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pyproj import Proj
from tqdm import tqdm
from scipy.interpolate import griddata
from PIL import Image
from skimage.measure import block_reduce
#import cv2 as cv
## Retrieve images and path information
imgs = [_ for _ in glob.glob(imgDir+'*.*') if _.endswith(('.npy'))]
imgs.sort()
pathDf = pd.read_csv(pathFile)
## Get corresponding GPS info matched with images
# Extract date and time from filenames
imgdates = [re.search('/20(.+?)_', path).group(1) for path in imgs] # Extract date from filename
imgtimes = [re.search('_(.+?)_', path).group(1) for path in imgs] # Extract time from filename
# Convert to unix datetime 
imgdatetimes = np.array([(datetime.timestamp(datetime(int('20'+imgdates[i][:2]),int(imgdates[i][2:4]),int(imgdates[i][4:6]),int(imgtimes[i][:2])+utcDiff,int(imgtimes[i][2:4]),int(imgtimes[i][4:6])))) for i in range(len(imgs))])*1000
# Get nearest GPS timestamp
gpstimes = [min(pathDf['timestamp'], key=(lambda list_value : abs(list_value - i))) for i in imgdatetimes]
# Create image dataframe
imgDf = pd.DataFrame(data=np.array([imgs,gpstimes]).transpose(),columns=['imgPath','timestamp'])
imgDf['timestamp'] = imgDf['timestamp'].astype(float)

# Merge with path dataframe
merged = imgDf.merge(pathDf[pathColumns], on='timestamp', how='left')
# Filter for acceptable measurement altitudes
merged = merged.loc[(merged[altCol] > minAlt) & (merged[altCol] < maxAlt)]
if imgTimes: merged = merged.loc[(merged['timestamp'] > imgTimes[0]) & (merged['timestamp'] < imgTimes[1])]

# Convert coordinates to UTM
myProj = Proj('+proj=utm +zone='+utmZone+', +'+hemisphere+' +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
merged['y'], merged['x'] = myProj(merged['longitude'].values, merged['latitude'].values)
merged<jupyter_output><empty_output><jupyter_text>Velocity Correction Attempt<jupyter_code># Coordinates from initial and speed
init = merged.iloc[0]
pi1 = pathDf.loc[pathDf['timestamp']==init['timestamp']].index[0]
meas = merged.iloc[1]
pi2 = pathDf.loc[pathDf['timestamp']==meas['timestamp']].index[0]
vels = pathDf.loc[pi1:pi2][['velocityX(mps)','velocityY(mps)']].mean()
vels
pathDf.loc[pathDf['timestamp']==init['timestamp']].index[0]
init<jupyter_output><empty_output><jupyter_text>Orientation Correction<jupyter_code># Orientation correction
dist = merged['altitude(m)']*np.tan(np.deg2rad(merged['pitch(deg)']))
merged['xc'], merged['yc'] = merged.x+(dist*np.sin(np.deg2rad(merged['yaw(deg)']))), merged.y+(dist*np.cos(np.deg2rad(merged['yaw(deg)'])))
if distFilt: merged = merged[abs(dist) < distFilt]
# Single image
val = 0.27
downsample = int(np.floor(val/pxSize))
imgExArr = np.load(merged.iloc[0]['imgPath'])
ds_array = imgExArr/255
rEx = block_reduce(ds_array, (downsample, downsample), np.mean)

# Get x,y coordinates for each pixel
yCoords, xCoords = [(np.array(range(rEx.shape[i]))-(rEx.shape[i]/2))*val+val/2 for i in [0,1]]
xq, yq = int(np.floor(len(xCoords)/4)), int(np.floor(len(yCoords)/4))
xmCoords, ymCoords = np.meshgrid(xCoords[xq:3*xq],yCoords[yq:3*yq])
xmcr, ymcr = xmCoords.reshape(-1,1), ymCoords.reshape(-1,1)

row = merged.iloc[5]
imgst = np.load(row.imgPath)[yq:3*yq,xq:3*xq]/255
total = np.concatenate((xmcr+row.xc, ymcr+row.yc,imgst.reshape(-1,1)),axis=-1)
extent = (np.ceil(min(total[:,0])), np.floor(max(total[:,0])),np.ceil(min(total[:,1])),np.floor(max(total[:,1])))
rs = np.dstack((total[:,2].reshape(2*yq,2*xq)))
plt.imshow(imgst,extent=extent)
row = merged.iloc[10]
imgst = np.load(row.imgPath)[yq:3*yq,xq:3*xq]/255
total2 = np.concatenate((xmcr+row.xc, ymcr+row.yc,imgst.reshape(-1,1)),axis=-1)
extent = (np.ceil(min(total2[:,0])), np.floor(max(total[:,0])),np.ceil(min(total[:,1])),np.floor(max(total[:,1])))
rs = np.dstack((total[:,2].reshape(2*yq,2*xq)))
plt.imshow(imgst,extent=extent)
alltot = np.vstack((total, total2))
extot = (np.ceil(min(alltot[:,0])), np.floor(max(alltot[:,0])),np.ceil(min(alltot[:,1])),np.floor(max(alltot[:,1])))# minx, maxx, miny, maxy = min(total[:,0]), max(total[:,0]), min(total[:,1]), max(total[:,1])

xsGrid = np.arange(extot[0], extot[1], val)
ysGrid = np.arange(extot[2], extot[3], val)
xsGridm, ysGridm = np.meshgrid(xsGrid, ysGrid)

method = 'nearest'
r,g,b = [griddata(alltot[1:,:2], alltot[1:,i], (xsGridm, ysGridm), method=method) for i in [2,3,4]]

# Plot 
fig, ax = plt.subplots(1,1,figsize=[10,15])
#extent = (np.ceil(min(total[:,0])), np.floor(max(total[:,0])),np.ceil(min(total[:,1])),np.floor(max(total[:,1])))# minx, maxx, miny, maxy = min(total[:,0]), max(total[:,0]), min(total[:,1]), max(total[:,1])

plt.imshow(np.dstack((r, g, b)),extent=extot)
# Loop over resolution to create combined image for each
# Initial setup and downsampling if desired
val = 0.27
downsample = int(np.floor(val/pxSize))
imgExArr = np.load(merged.iloc[0]['imgPath'])
ds_array = imgExArr/255
rEx = block_reduce(ds_array, (downsample, downsample), np.mean)

# Get x,y coordinates for each pixel
yCoords, xCoords = [(np.array(range(rEx.shape[i]))-(rEx.shape[i]/2))*val+val/2 for i in [0,1]]
xq, yq = int(np.floor(len(xCoords)/4)), int(np.floor(len(yCoords)/4))
xmCoords, ymCoords = np.meshgrid(xCoords[xq:3*xq],yCoords[yq:3*yq])
xmcr, ymcr = xmCoords.reshape(-1,1), ymCoords.reshape(-1,1)

# Loop over images
skip=1
alltot = [0,0,0]
for index, row in tqdm(merged[::skip].iterrows()):
    imgst = np.load(row.imgPath)[yq:3*yq,xq:3*xq]
    total = np.concatenate((xmcr+row.xc, ymcr+row.yc,imgst.reshape(-1,1)),axis=-1)
    alltot = np.vstack((alltot,total))

alltot = alltot[1:,:]
extot = (np.ceil(min(alltot[:,0])), np.floor(max(alltot[:,0])),np.ceil(min(alltot[:,1])),np.floor(max(alltot[:,1])))
xsGrid = np.arange(extot[0], extot[1], val)
ysGrid = np.arange(extot[2], extot[3], val)
xsGridm, ysGridm = np.meshgrid(xsGrid, ysGrid)

method = 'nearest'
gridded = griddata(alltot[1:,:2], alltot[1:,2], (xsGridm, ysGridm), method=method)

# Plot 
fig, ax = plt.subplots(1,1,figsize=[10,10])

a = ax.imshow(gridded,extent=extot, cmap='RdYlGn')
fig.colorbar(a, ax=ax)
ax.set_title('Interpolation to grid using every '+str(skip)+' images')
gridded
# val = 6
cmap='hot'
# method = 'nearest'
# gridded = griddata(alltot[1:,:2], alltot[1:,2], (xsGridm, ysGridm), method=method)
# downsample = int(np.floor(val/pxSize))
# ds = block_reduce(gridded, (downsample, downsample), np.mean)

# Plot 
fig, axs = plt.subplots(1,2,figsize=[15,5])

a = axs[0].imshow(gridded,extent=extot, cmap=cmap)
b = axs[1].imshow(ds,extent=extot, cmap=cmap, vmin=gridded.min(),vmax=gridded.max())
plt.tight_layout(h_pad=1)
fig.subplots_adjust(right=0.9)
cbar_ax = fig.add_axes([0.92, 0.1, 0.03, 0.85])
c = fig.colorbar(a, cax=cbar_ax)
c.set_label('Temperature ($^{\circ}$C)')

axs[0].set_title('Full resolution'), axs[1].set_title(str(val)+' m resolution')
fig.savefig('plots/thermal_rollcorr.png')
# Plot example downsampled image
for i, val in enumerate(resolutions):
    downsample = int(np.floor(val/pxSize))
    # first, change to 0-1
    ds_array = imgExArr/255
    r = block_reduce(ds_array[:, :, 0], (downsample, downsample), np.mean)
    g = block_reduce(ds_array[:, :, 1], (downsample, downsample), np.mean)
    b = block_reduce(ds_array[:, :, 2], (downsample, downsample), np.mean)
    ds_array = np.stack((r, g, b), axis=-1)
    if i<2:
        c = axs[1,i].imshow(ds_array)
        axs[1,i].set_title('Image {} downsampled to {} m resolution'.format(str(imageNum),val))
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2057953144266338, 'fraction_numerical': 0.03193588162762022, 'mean_word_length': 6.905458089668616, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '636430', 'license': 'MIT', 'n_tokens_mistral': 3464, 'n_tokens_neox': 3361, 'n_words': 585, 'path': 'numpy-thermal-addition.ipynb', 'repo_name': 'UranKhatolaaa/s4-thermal', 'size': '8110'}