Converting JPGs from M2EA to TIFF but Ortho-Mosaic fails

[FilipRothaut]

Hello everyone, I've written the following code to convert the R-JPGs spitted out by the DJI M2EA into TIFFs, read the tmperature and put the output into an mosaic dataset to display an AOI in a small map. 

My problem now is, that the map created in the process seems to turn out 90 degress rotated and I can not find a reason to why it should be displayed wrong. If I insert the original JPGs in Drone2Map it turns out fine, but something is missing in my arcpy process.

Any help would be really appreciated!

import os
import subprocess
import numpy as np
import imageio
import arcpy
import pandas as pd
from tkinter import Tk
from tkinter.filedialog import askdirectory
from tifffile import TiffFile
from PIL import Image
from pathlib import Path
import time

for root, dirs, files in os.walk(os.path.abspath(os.path.join(os.getcwd(), os.pardir))):
    for dir in dirs:
        if "release_x64" in dir:
            sdk_dir = (os.path.join(root, dir))
os.chdir(sdk_dir)

root1 = root
root = Tk()
root.withdraw()
folder_path = askdirectory()
print("askdir: " + folder_path)
root = root1

in_dir = folder_path

out_dir = "ir_calib/"
out_dir = os.path.join(in_dir, out_dir)
os.makedirs(out_dir, exist_ok=True)
for root, dirs, files in os.walk(sdk_dir):
    for file in files:
        if "exiftool" in file:
            exiftool = os.path.join(sdk_dir, file)


def get_exif_data(file_path):
    exif_data = {}
    result = subprocess.run([exiftool, file_path], stdout=subprocess.PIPE)
    output = result.stdout.decode("utf-8")
    # Parse the output into a dictionary
    for line in output.splitlines():
        key = line.strip().split(":")[0]
        value = line.strip().split(":")[1]
        exif_data[key.strip()] = value.strip()

    return exif_data


def get_length(image_path):
    result = subprocess.run([exiftool, "-ImageLength", image_path], stdout=subprocess.PIPE)
    output = result.stdout.decode("utf-8")
    result2 = output.split(":")
    length2 = int(result2[1])
    return length2


def get_width(image_path):
    result = subprocess.run([exiftool, "-ImageWidth", image_path], stdout=subprocess.PIPE)
    output = result.stdout.decode("utf-8")
    result2 = output.split(":")
    width2 = int(result2[1])
    return width2


def copy_exif_data(src_file, dest_file):
    exif_data = tag_list
    args = [exiftool, "-TagsFromFile", src_file]
    for data in exif_data:
        if data == "FocalLength":
            args.extend(["-{}={}".format(data, "9")])
        elif data == "FocalLengthIn35mmFormat":
            args.extend(["-{}={}".format(data, "38")])
        #elif data == "Orientation":
            #args.extend(["-{}={}".format(data, "Rotate 180")])
        else:
            args.extend(["-{}<={}".format(data, data)])
    args.append(dest_file)

    result = subprocess.run(args, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    if result.returncode == 0:
        print("EXIF data copied successfully.")
    else:
        print("Error:", result.stderr.decode().strip())


in_files = []
for root, dirs, files in os.walk(in_dir):
    for file in files:
        if "_T" in file:
            in_files.append(file)

path = str(folder_path)


emissivity = 0.957  # Default: 1.0 range: 0.1-1.0 https://royalsocietypublishing.org/doi/pdf/10.1098/rsos.181281 (e.g., average 0.957 for vegetation)
humidity = 38  # Default: 70 %  range: 20-100 %
distance = 25  # Default: 5 m   range: 1-25 m    # Altitude - object height (set to 25 if camera-target-distance > 25 m)

for i in range(len(in_files)):
    # Calibration to celsius
    with open(os.path.join(in_dir, in_files[i]), 'rb') as f:
        in_name = in_files[i]
    out_name = os.path.splitext(in_files[i])[0] + ".raw"

    in_path = path + "\\" + in_name
    os.system(
        f'dji_irp.exe -s {in_path} -a measure -o {out_name} --measurefmt float32 --emissivity {emissivity} --humidity {humidity} --distance {distance}')

    exif_data = get_exif_data(in_path)
    length = int(str(exif_data.get("Image Height")))
    width = int(str(exif_data.get("Image Width")))

    with open(out_name, "rb") as f:
        raw_data = f.read()  # image length und image width müssen zurecht geschnitten werden
    image_matrix = np.frombuffer(raw_data, dtype='float32').reshape(length, width)
    out_name_tif = out_name[:-4] + ".tif"
    out_path = out_dir + "\\" + out_name_tif
    imageio.imwrite(out_path, image_matrix)
    img = imageio.imread(out_path)
    img = img.astype("float16")
    imageio.imwrite(out_path, img)


    with open(out_name, "rb") as f:
        tag_list = ["Model", "Make", "ExifOffset", "FocalLength", "FocalLengthIn35mmFormat", "DigitalZoomRatio",
                    "GPSAltitude", "RtkStdLon", "RtkStdLat", "RtkStdHgt", "FlightRollDegree", "FlightYawDegree", "FlightRollDegree",
                                                             "GimbalDegree","FlightDegree","GPSLongitude", "GPS Longitude Ref", "GPSLatitude",
                    "GPS Latitude Ref", "Orientation", "ApertureValue", "GimbalPitchDegree", "GimbalRollDegree"]

    copy_exif_data(in_path, out_path)

    f.close()

for root, dirs, files in os.walk(out_dir):
    for file in files:
        if "_original" in file:
            os.remove(os.path.join(out_dir, file))

for root, dirs, files in os.walk(sdk_dir):
    for file in files:
        if file.endswith(".raw"):
            os.remove(os.path.join(sdk_dir, file))


folder_path = out_dir
i = 0
tif_list = []
u_half = []
l_half = []
max_list = []
q1 = []
q2 = []
q3 = []
q4 = []
max2 = 0


for filename in os.listdir(folder_path):
    if filename.endswith(".tif"):
        with TiffFile(os.path.join(folder_path, filename)) as tiff:
            page = tiff.pages[0]
            image_data = page.asarray()
            i = i + 1
            pixels = []
            for row in range(image_data.shape[0]):
                for col in range(image_data.shape[1]):
                    # Get the temperature value at the current pixel
                    temperature = image_data[row, col]
                    pixels.append((row, col, temperature))
            dfp = pd.DataFrame(pixels, columns=["row", "col", "temperature"])
            max1 = max(dfp["temperature"])
            max_list.append(max1)
max1 = max(max_list)
i = 0

print("")
print("Maximum is: " + str(max1))
print("")

for filename in os.listdir(folder_path):
    # Only process TIFF files
    if filename.endswith(".tif"):
        # Open the TIFF file
        with TiffFile(os.path.join(folder_path, filename)) as tiff:
            # Get the first page of the TIFF file
            page = tiff.pages[0]

            i = i + 1
            # Get the image data from the page
            image_data = page.asarray()
            # print(tiff)
            # Initialize an empty list to store the pixels for the current TIFF file
            pixels = []

            # Loop over each row and column of the image data
            for row in range(image_data.shape[0]):
                for col in range(image_data.shape[1]):
                    # Get the temperature value at the current pixel
                    temperature = image_data[row, col]
                    pixels.append((row, col, temperature))

            df = pd.DataFrame(pixels, columns=["row", "column", "temperature"])
            median = df['temperature'].median()
            u_half = df[df['temperature'] >= median].values.tolist()
            dfu = pd.DataFrame(u_half, columns=["row", "column", "temperature"])
            medi2 = dfu['temperature'].median()
            q1 = dfu[dfu['temperature'] >= medi2].values.tolist()
            dfq = pd.DataFrame(q1, columns=["row", "column", "temperature"])
            mean = dfq['temperature'].mean()
            std = dfq['temperature'].std()
            threshold = mean + 3 * std
            extremes = dfq[(dfq['temperature'] > threshold) | (dfq['temperature'] < mean - 3 * std)]
            dfe = pd.DataFrame(extremes, columns=["row", "column", "temperature"])
            medi4 = dfe['temperature'].median()
            interest = dfe[dfe['temperature'] >= medi4].values.tolist()
            dfi = pd.DataFrame(interest, columns=["row", "column", "temperature"])
            top = max(dfi["temperature"])
            if top > max1 - 2:
                tif_list.append(str(tiff))


tif2 = []
for elements in tif_list:
    start = elements.rfind("D")
    end = elements.rfind("f") + 1
    filename = elements[start:end]
    tif2.append(filename)
tif2.sort()
print("tif2:")
print(tif2)

interpol = []
interpol = tif2
print("interpol:")
print(interpol)

path2 = os.path.abspath(os.path.join(folder_path, os.pardir))
path3 = os.path.abspath(os.path.join(path2, os.pardir))
os.chdir(os.path.join(path3, "Template_Project"))
count = 0

for i in interpol:
    filename = interpol[count]
    count = count + 1
    jpg_files = sorted([f for f in os.listdir(in_dir) if f.endswith("_T.JPG")])
    tif_files = sorted([f for f in os.listdir(out_dir) if f.endswith("_T.tif")])


    def get_neighbors(filename):
        index = tif_files.index(filename)
        if index == 0:
            neighbors = tif_files[1:3]
            neighbors.append(tif_files[index])
        elif index == len(jpg_files) - 1:
            neighbors = tif_files[-3:-1]
            neighbors.append(tif_files[index])
        else:
            neighbors = tif_files[index - 1:index + 2]
        return sorted(neighbors)


    neighbors = get_neighbors(filename)
    print("nachbarn")
    print(neighbors)

    with arcpy.EnvManager(scratchWorkspace=os.path.join(os.getcwd(), "Template.gdb"),
                          workspace=os.path.join(os.getcwd(), "Template.gdb")):
        # To allow overwriting outputs change overwriteOutput option to True.
        arcpy.env.overwriteOutput = True

        Template_gdb = os.path.join(os.getcwd(), "Template.gdb")
        raster1 = arcpy.Raster(str(os.path.join(out_dir, neighbors[0])))
        raster2 = arcpy.Raster(str(os.path.join(out_dir, neighbors[1])))
        raster3 = arcpy.Raster(str(os.path.join(out_dir, neighbors[2])))
        Area_of_Interest = "in_memory\\feature_set1"
        M2EA_camfile_csv = os.path.join(os.getcwd(), "M2EA_camfile.csv")

        # Process: Create Mosaic Dataset (Create Mosaic Dataset) (management)
        mosaik = arcpy.management.CreateMosaicDataset(in_workspace=Template_gdb, in_mosaicdataset_name="mosaik",
                                                      coordinate_system="PROJCS[\"ETRS_1989_UTM_Zone_32N\",GEOGCS[\"GCS_ETRS_1989\",DATUM[\"D_ETRS_1989\",SPHEROID[\"GRS_1980\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",500000.0],PARAMETER[\"False_Northing\",0.0],PARAMETER[\"Central_Meridian\",9.0],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",0.0],UNIT[\"Meter\",1.0]]",
                                                      num_bands=None, pixel_type="", product_definition="NONE",
                                                      product_band_definitions=[])[0]

        # Process: Add Rasters To Mosaic Dataset (Add Rasters To Mosaic Dataset) (management)
        mosaik_2_ = arcpy.management.AddRastersToMosaicDataset(in_mosaic_dataset=mosaik, raster_type="UAV/UAS",
                                                               input_path=[raster1, raster2, raster3],
                                                               update_cellsize_ranges="UPDATE_CELL_SIZES",
                                                               update_boundary="UPDATE_BOUNDARY",
                                                               update_overviews="NO_OVERVIEWS",
                                                               maximum_pyramid_levels=None, maximum_cell_size=0,
                                                               minimum_dimension=1500, spatial_reference="",
                                                               filter="*.*", sub_folder="SUBFOLDERS",
                                                               duplicate_items_action="ALLOW_DUPLICATES",
                                                               build_pyramids="NO_PYRAMIDS",
                                                               calculate_statistics="NO_STATISTICS",
                                                               build_thumbnails="NO_THUMBNAILS",
                                                               operation_description="",
                                                               force_spatial_reference="NO_FORCE_SPATIAL_REFERENCE",
                                                               estimate_statistics="ESTIMATE_STATISTICS",
                                                               aux_inputs=[["CameraFile", M2EA_camfile_csv],
                                                                           ["IsAltitudeFlightHeight", ""]],
                                                               enable_pixel_cache="NO_PIXEL_CACHE", cache_location="")[
            0]

        # Process: Compute Mosaic Candidates (Compute Mosaic Candidates) (management)
        Derived_Mosaic_Dataset = \
            arcpy.management.ComputeMosaicCandidates(in_mosaic_dataset=mosaik_2_, maximum_overlap=0.6,
                                                     maximum_area_loss=0.05)[0]

        # Process: Build Seamlines (Build Seamlines) (management)
        Updated_Mosaic_Dataset_2_ = \
            arcpy.management.BuildSeamlines(in_mosaic_dataset=Derived_Mosaic_Dataset, cell_size=[],
                                            sort_method="NORTH_WEST", sort_order="ASCENDING", order_by_attribute="",
                                            order_by_base_value=None, view_point="", computation_method="RADIOMETRY",
                                            blend_width=None, blend_type="BOTH", request_size=1000,
                                            request_size_type="PIXELS", blend_width_units="PIXELS",
                                            area_of_interest=Area_of_Interest, where_clause="",
                                            update_existing="IGNORE_EXISTING", min_region_size=100,
                                            min_thinness_ratio=0.05, max_sliver_size=20)[0]

        # Process: Copy Raster (Copy Raster) (management)
        name = "Ausschnitt" + str(count)
        Output_Raster_Dataset = str(os.path.join(Template_gdb, name))
        arcpy.management.CopyRaster(in_raster=Updated_Mosaic_Dataset_2_, out_rasterdataset=Output_Raster_Dataset,
                                    config_keyword="", background_value=None, nodata_value="",
                                    onebit_to_eightbit="NONE", colormap_to_RGB="NONE", pixel_type="",
                                    scale_pixel_value="NONE", RGB_to_Colormap="NONE", format="", transform="NONE",
                                    process_as_multidimensional="CURRENT_SLICE",
                                    build_multidimensional_transpose="NO_TRANSPOSE")
        aprx = arcpy.mp.ArcGISProject(os.path.join(os.getcwd(), "Template.aprx"))
        aprx.saveACopy(os.path.join(os.getcwd(), "Project.aprx"))
        aprx = arcpy.mp.ArcGISProject(os.path.join(os.getcwd(), "Project.aprx"))
        mxd = aprx.listMaps("Map")[0]
        lyt = aprx.listLayouts("Layout")[0]
        Lyr1 = "Lyr1"
        Output_Layer = arcpy.management.MakeRasterLayer(Output_Raster_Dataset, "AreaofInterest")
        lyrx1 = arcpy.management.SaveToLayerFile(Output_Layer, "Lyr1.lyrx", "ABSOLUTE")
        lyrx1 = arcpy.mp.LayerFile(lyrx1)
        mxd.addLayer(lyrx1, 'TOP')
        lyr = mxd.listLayers("AreaofInterest")[0]

        sym = lyr.symbology
        if hasattr(sym, 'colorizer'):
            if sym.colorizer.type == 'RasterStretchColorizer':
                sym.colorizer.colorRamp = aprx.listColorRamps('Inferno')[0]
                lyr.symbology = sym

        elem = lyt.listElements("MAPFRAME_ELEMENT", "*")[0]
        elem.camera.setExtent(elem.getLayerExtent(lyr, False))
        aprx.saveACopy(os.path.join(os.getcwd(), "Project.aprx"))

        output_pdf = out_dir + "karte" + str(count) + ".pdf"
        lyt.exportToPDF(output_pdf)
        for root, dirs, files in os.walk(os.getcwd()):
            for file in files:
                if file.endswith(".lyrx"):
                    os.remove(os.path.join(os.getcwd(), file))