New exif epsg

main_dg.py

@@ -1,25 +1,33 @@
 import os
-import numpy as np
 import time
-from module.ExifData import *
-from module.EoData import *
+import numpy as np
+#from module.ExifData import *
+from module.ExifData import get_exif, restoreOrientation
+from module.EoData import rpy_to_opk, Rot3D
 from module.Boundary import boundary
 from module.BackprojectionResample import rectify_plane_parallel, createGeoTiff
-from rich.console import Console
+from pyproj import Transformer
 from rich.table import Table
+from rich.console import Console
+import cv2
 
 console = Console()
 
-input_folder = 'Data'
-ground_height = 0   # unit: m
-sensor_width = 6.3  # unit: mm, Mavic
-# sensor_width = 13.2  # unit: mm, P4RTK
-# sensor_width = 17.3  # unit: mm, Inspire
-epsg = 5186     # editable
-gsd = 0.1   # unit: m, set 0 to compute automatically
+INPUT_FOLDER = 'Data'
+GROUND_HEIGHT = 0   # unit: m
+gsd = 0.1   # unit: m
+AUTO_GSD = True  # True to compute automatically gsd per image
+
+
+# Modify next line with local epsg of your data
+# Find your local EPSG (https://epsg.io/)
+EPSG_OUT = 5186  # Korea 2000 / Central Belt 2010
+
+# Most of the time data location is GPS based with WGS-84
+ESPSG_IN = 4326  # WGS-84 [dd.mm.sssd]
 
 if __name__ == '__main__':
-    for root, dirs, files in os.walk(input_folder):
+    for root, dirs, files in os.walk(INPUT_FOLDER):
         files.sort()
         for file in files:
             image_start_time = time.time()
@@ -29,79 +37,95 @@
             file_path = root + '/' + file
             dst = './' + filename
 
-            if extension == '.JPG' or extension == '.jpg':
+            if extension in ['.JPG', '.jpg']  :
                 print('Georeferencing - ' + file)
                 start_time = time.time()
                 image = cv2.imread(file_path, -1)
 
                 # 1. Extract metadata from a image
-                focal_length, orientation, eo, maker = get_metadata(file_path)  # unit: m, _, ndarray
+                Camera = get_exif(file_path)
+                print(f"[Camera['Roll'] {Camera['Roll']} Camera['Pitch'] {Camera['Pitch']} Camera['Yaw'] {Camera['Yaw']}")
 
                 # 2. Restore the image based on orientation information
-                restored_image = restoreOrientation(image, orientation)
-
-                image_rows = restored_image.shape[0]
-                image_cols = restored_image.shape[1]
+                restored_image = restoreOrientation(image, Camera['Orientation'])
+                image_rows, image_cols, _ = restored_image.shape
 
-                pixel_size = sensor_width / image_cols  # unit: mm/px
+                pixel_size = Camera['Size'][0] / image_cols  # unit: mm/px
                 pixel_size = pixel_size / 1000  # unit: m/px
 
-                eo = geographic2plane(eo, epsg)
-                opk = rpy_to_opk(eo[3:], maker)
-                eo[3:] = opk * np.pi / 180   # degree to radian
+                # Get Local coordinates
+                coordonates_transformer = Transformer.from_crs(f'epsg:{ESPSG_IN}', f'epsg:{EPSG_OUT}')
+
+                Camera['X'], Camera['Y'] = coordonates_transformer.transform(Camera['Lat'], Camera['Lon'], errcheck=True)
+                # OPK from RPY
+                Camera['Omega'], Camera['Phi'], Camera['Kappa'] = rpy_to_opk(Camera)
+
+                # already converted to math.radians in rpy_to_opk
+                # eo[3:] = opk * np.pi / 180   # degree to radian
+                eo = np.array([Camera['X'], Camera['Y'], Camera['RelAltitude'],
+                              Camera['Omega'], Camera['Phi'], Camera['Kappa']])
+
                 R = Rot3D(eo)
 
                 console.print(
-                    f"EOP: {eo[0]:.2f} | {eo[1]:.2f} | {eo[2]:.2f} | {eo[3]:.2f} | {eo[4]:.2f} | {eo[5]:.2f}\n"
-                    f"Focal Length: {focal_length * 1000:.2f} mm, Maker: {maker}",
+                    f"EOP: Longitude {eo[0]:.2f} | Latitude {eo[1]:.2f} | Altitude {eo[2]:.2f} | Omega {eo[3]:.2f} | Phi {eo[4]:.2f} | Kappa {eo[5]:.2f}\n"
+                    f"Focal Length: {Camera['Focalength'] * 1000:.2f} mm",
                     style="blink bold red underline")
                 georef_time = time.time() - start_time
-                console.print(f"Georeferencing time: {georef_time:.2f} sec", style="blink bold red underline")
+                console.print(
+                    f"Georeferencing time: {georef_time:.2f} sec", style="blink bold red underline")
 
                 print('DEM & GSD')
                 start_time = time.time()
                 # 3. Extract a projected boundary of the image
-                bbox = boundary(restored_image, eo, R, ground_height, pixel_size, focal_length)
+                GROUND_HEIGHT = Camera['AbsAltitude']
+                bbox = boundary(restored_image, eo, R, GROUND_HEIGHT,
+                                pixel_size, Camera['Focalength'])
 
                 # 4. Compute GSD & Boundary size
                 # GSD
-                if gsd == 0:
-                    gsd = (pixel_size * (eo[2] - ground_height)) / focal_length  # unit: m/px
+                if AUTO_GSD:
+                    gsd = (pixel_size * (Camera['RelAltitude'] - Camera['AbsAltitude'] )) / Camera['Focalength']  # unit: m/px
+
                 # Boundary size
                 boundary_cols = int((bbox[1, 0] - bbox[0, 0]) / gsd)
                 boundary_rows = int((bbox[3, 0] - bbox[2, 0]) / gsd)
 
-                console.print(f"GSD: {gsd * 100:.2f} cm/px", style="blink bold red underline")
+                console.print(f"GSD: {gsd * 100:.2f} cm/px",style="blink bold red underline")
                 dem_time = time.time() - start_time
                 console.print(f"DEM time: {dem_time:.2f} sec", style="blink bold red underline")
 
                 print('Rectify & Resampling')
                 start_time = time.time()
-                b, g, r, a = rectify_plane_parallel(bbox, boundary_rows, boundary_cols, gsd, eo, ground_height,
-                                                    R, focal_length, pixel_size, image)
+                b, g, r, a = rectify_plane_parallel(bbox, boundary_rows, boundary_cols, gsd, eo, GROUND_HEIGHT,
+                                                    R, Camera['Focalength'], pixel_size, image)
                 rectify_time = time.time() - start_time
                 console.print(f"Rectify time: {rectify_time:.2f} sec", style="blink bold red underline")
 
                 # 8. Create GeoTiff
                 print('Save the image in GeoTiff')
                 start_time = time.time()
-                createGeoTiff(b, g, r, a, bbox, gsd, epsg, boundary_rows, boundary_cols, dst)
+                createGeoTiff(b, g, r, a, bbox, gsd, EPSG_OUT,
+                              boundary_rows, boundary_cols, dst)
                 # create_pnga_optical(b, g, r, a, bbox, gsd, epsg, dst)   # for test
                 write_time = time.time() - start_time
-                console.print(f"Write time: {write_time:.2f} sec", style="blink bold red underline")
+                console.print(
+                    f"Write time: {write_time:.2f} sec", style="blink bold red underline")
 
                 processing_time = time.time() - image_start_time
-                console.print(f"Process time: {processing_time:.2f} sec", style="blink bold red underline")
+                console.print(
+                    f"Process time: {processing_time:.2f} sec", style="blink bold red underline")
 
                 table = Table(show_header=True, header_style="bold magenta")
-                table.add_column("Image", style="dim", width=12)
+                table.add_column("Image", style="dim")
                 table.add_column("Georeferencing", justify="right")
                 table.add_column("DEM", justify="right")
                 table.add_column("Rectify", justify="right")
                 table.add_column("Write", justify="right")
                 table.add_column("Processing", justify="right")
                 table.add_row(
-                    filename, str(round(georef_time, 5)), str(round(dem_time, 5)), str(round(rectify_time, 5)),
+                    filename, str(round(georef_time, 5)), str(
+                        round(dem_time, 5)), str(round(rectify_time, 5)),
                     str(round(write_time, 5)), str(round(processing_time, 5))
                 )
                 console.print(table)

module/EoData.py

@@ -1,7 +1,6 @@
 import numpy as np
 import math
 from osgeo.osr import SpatialReference, CoordinateTransformation
-import osgeo
 
 def readEO(path):
     eo_line = np.genfromtxt(path, delimiter='\t',
@@ -18,33 +17,6 @@ def readEO(path):
 
     return eo
 
-def geographic2plane(eo, epsg=5186):
-    # Define the Plane Coordinate System (e.g. 5186)
-    plane = SpatialReference()
-    plane.ImportFromEPSG(epsg)
-
-    # Define the wgs84 system (EPSG 4326)
-    geographic = SpatialReference()
-    geographic.ImportFromEPSG(4326)
-
-    coord_transformation = CoordinateTransformation(geographic, plane)
-
-    # Check the transformation for a point close to the centre of the projected grid
-    if int(osgeo.__version__[0]) >= 3:  # version 3.x
-        if str(epsg).startswith("51"):  # for Korean CRS only (temporarily) ... TODO: for whole CRS
-            # Transform(y,x) will return y, x (Northing, Easting)
-            yx = coord_transformation.TransformPoint(float(eo[1]), float(eo[0]))  # The order: Lat, Lon
-            eo[0:2] = yx[0:2][::-1]
-        else:
-            # Transform(y,x) will return x,y (Easting, Northing)
-            xy = coord_transformation.TransformPoint(float(eo[1]), float(eo[0]))  # The order: Lat, Lon
-            eo[0:2] = xy[0:2]
-    else:  # version 2.x
-        # Transform(x,y) will return x,y (Easting, Northing)
-        xy = coord_transformation.TransformPoint(float(eo[0]), float(eo[1]))  # The order: Lon, Lat
-        eo[0:2] = xy[0:2]
-
-    return eo
 
 def tmcentral2latlon(eo):
     # Define the TM central coordinate system (EPSG 5186)
@@ -118,9 +90,18 @@ def rot_2d(theta):
     return np.array([[np.cos(theta), np.sin(theta)],
                      [-np.sin(theta), np.cos(theta)]])
 
-def rpy_to_opk(rpy, maker=""):
-    if maker == "samsung":
-        roll_pitch = np.empty_like(rpy[0:2])
+def rpy_to_opk(Camera):
+    """
+    Convert Roll Pitch Yaw from Camera in degrees
+    to OPK (Omega Phi Kappa) np.array in radians 
+    return the OPK np.array
+    """
+    rpy = [Camera['Roll'], Camera['Pitch'],Camera['Yaw']]
+    maker = Camera['Make']
+
+    roll_pitch = np.empty_like(rpy[0:2])
+
+    if "samsung" in maker:
 
         roll_pitch[0] = -rpy[1]
         roll_pitch[1] = -rpy[0]
@@ -129,7 +110,6 @@ def rpy_to_opk(rpy, maker=""):
         kappa = -rpy[2] - 90
         return np.array([float(omega_phi[0, 0]), float(omega_phi[1, 0]), kappa])
     else:
-        roll_pitch = np.empty_like(rpy[0:2])
         roll_pitch[0] = 90 + rpy[1]
         if 180 - abs(rpy[0]) <= 0.1:
             roll_pitch[1] = 0
@@ -138,4 +118,5 @@ def rpy_to_opk(rpy, maker=""):
 
         omega_phi = np.dot(rot_2d(rpy[2] * np.pi / 180), roll_pitch.reshape(2, 1))
         kappa = -rpy[2]
-        return np.array([float(omega_phi[0, 0]), float(omega_phi[1, 0]), kappa])
+
+    return np.array(np.radians([float(omega_phi[0, 0]), float(omega_phi[1, 0]), kappa]))

module/ExifData.py

@@ -1,6 +1,8 @@
 import cv2
 import numpy as np
 import pyexiv2
+import exiftool
+from os.path import exists
 
 
 def restoreOrientation(image, orientation):
@@ -15,6 +17,7 @@ def restoreOrientation(image, orientation):
 
     return restored_image
 
+
 def rotate(image, angle):
     # https://www.pyimagesearch.com/2017/01/02/rotate-images-correctly-with-opencv-and-python/
 
@@ -41,13 +44,15 @@ def rotate(image, angle):
     rotated_mat = cv2.warpAffine(image, rotation_mat, (bound_w, bound_h))
     return rotated_mat
 
+
 def get_metadata(input_file):
     img = pyexiv2.Image(input_file)
 
     exif = img.read_exif()
     xmp = img.read_xmp()
 
-    focal_length = convert_string_to_float(exif['Exif.Photo.FocalLength']) / 1000    # unit: m
+    focal_length = convert_string_to_float(
+        exif['Exif.Photo.FocalLength']) / 1000    # unit: m
     orientation = int(exif['Exif.Image.Orientation'])
     maker = exif["Exif.Image.Make"]
 
@@ -82,14 +87,100 @@ def get_metadata(input_file):
 
     return focal_length, orientation, eo, maker
 
-def convert_dms_to_deg(dms):
-    dms_split = dms.split(" ")
-    d = convert_string_to_float(dms_split[0])
-    m = convert_string_to_float(dms_split[1]) / 60
-    s = convert_string_to_float(dms_split[2]) / 3600
-    deg = d + m + s
-    return deg
-
-def convert_string_to_float(string):
-    str_split = string.split('/')
-    return int(str_split[0]) / int(str_split[1])    # unit: mm
+
+def get_exif(path_name):
+
+    if (not exists(path_name)):
+        print(f'Error {path_name} is not a path to an existing file')
+        exit(0)
+
+    print(f'Extracting Exif data  from {path_name}')
+    with exiftool.ExifToolHelper() as et:
+        metadata = et.get_metadata(path_name)[0]
+
+    Camera = {}
+    Camera['FileName'] = path_name
+    if 'File:ImageWidth' in metadata:
+        Camera['ImageWidth'] = int(metadata['File:ImageWidth'])
+        Camera['ImageHeight'] = int(metadata['File:ImageHeight'])
+    elif 'EXIF:ImageWidth' in metadata:
+        Camera['ImageWidth'] = int(metadata['EXIF:ImageWidth'])
+        Camera['ImageHeight'] = int(metadata['EXIF:ImageHeight'])
+
+    Camera['Focalength'] = float(metadata['EXIF:FocalLength']) / 1000
+    Camera['Lat'] = float(metadata['EXIF:GPSLatitude'])
+    Camera['Lon'] = float(metadata['EXIF:GPSLongitude'])
+    Camera['LatRef'] = metadata['EXIF:GPSLatitudeRef']
+    Camera['LonRef'] = metadata['EXIF:GPSLongitudeRef']
+    Camera['Make'] = metadata['EXIF:Make']
+
+    # __import__("IPython").embed()
+
+    Camera['Orientation'] = int(metadata['EXIF:Orientation'])
+
+    if 'S' in Camera['LatRef']:
+        Camera['Lat'] = -Camera['Lat']
+
+    if 'W' in Camera['LonRef']:
+        Camera['Lon'] = -Camera['Lon']
+
+    if 'DJI' in Camera['Make']:
+        Camera['AbsAltitude'] = float(metadata['XMP:AbsoluteAltitude'])
+        Camera['RelAltitude'] = float(metadata['XMP:RelativeAltitude'])
+        Camera['Yaw'] = float(metadata['XMP:GimbalYawDegree'])
+        Camera['Pitch'] = float(metadata['XMP:GimbalPitchDegree'])
+        Camera['Roll'] = float(metadata['XMP:GimbalRollDegree'])
+
+        # Sensor Size
+        # Phantom 4 RTK
+        if 'XMP:Model' in metadata:
+            # Phantom 4 RTK
+            if 'FC6310R' in metadata['XMP:Model']:
+                Camera['Focalength'] = 8.8 / 1000
+                # 5472 x 3648 (3:2) mode
+                if Camera['ImageWidth'] == 5472:
+                    Camera['Size'] = (13.2, 8.8)   # in mm
+                # 4864 x 3648 (4:3) mode
+                elif Camera['ImageWidth'] == 4864:
+                    Camera['Size'] = (13.2 * 4864 / 5472,
+                                      8.8 * 4864 / 5472)  # unit: mm
+            # Mavic Pro
+            elif 'FC220' in metadata['XMP:Model']:
+                Camera['Focalength'] = 4.74 / 1000
+                Camera['Size'] = (6.16, 4.55)
+            else :
+                print('Drone not yet in the code.')
+                print('https://www.djzphoto.com/blog/2018/12/5/dji-drone-quick-specs-amp-comparison-page')
+
+    elif 'samsung' in metadata['EXIF:Make']:
+        #Camera['AbsAltitude'] = float(metadata['XMP:AbsoluteAltitude'])
+        Camera['RelAltitude'] = float(metadata['EXIF:GPSAltitude'])
+        Camera['Yaw'] = float(metadata['Xmp:Yaw']) * 180 / np.pi
+        Camera['Pitch'] = float(metadata['Xmp:Pitch']) * 180 / np.pi
+        Camera['Roll'] = float(metadata['Xmp:Roll']) * 180 / np.pi
+
+    else:
+        print('Your UAV/Camera is not supported yet')
+        exit(0)
+#        altitude = 0
+#        roll = 0
+#        pitch = 0
+#        yaw = 0
+    #for tag in metadata:
+    #    print(f'{tag} {metadata[tag]}')
+
+    return Camera
+    #eo = np.array([longitude, latitude, altitude, roll, pitch, yaw])
+
+
+def main():
+    file = './Data/DJI_0386.JPG'
+    file = './tests/20191011_074853.JPG'
+
+    #file = './mayo/100_0138_0001.JPG'
+    print(f'Reading EXIF of {file}')
+    print(get_exif(file))
+
+
+if __name__ == "__main__":
+    main()