deep-crop/gis-svg-stitcher.py

import os
import argparse
from wand.image import Image
import lxml.etree as ET
import copy
import math


from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.database import query_utm_crs_info
from pyproj import Transformer



arg_parser = argparse.ArgumentParser(description='Place drone FLIR-tiles into a SVG in order to edit them in Inkscape')


arg_parser.add_argument('Input',
                       metavar='input_directory',
                       type=str,
                       help='Path where the FLIR tiles are')


arg_parser.add_argument(
    '--base_rotation', action='store', default=115,
    help="Base orientation of drone in degrees (0-360) Defaults to 115",
    type=int, dest='base_rotation'
)

arg_parser.add_argument(
    '--scale', action='store', default=15,
    help="Scaling (higher number leads to bigger canvas and less dense tiles) (defaults to 15)",
    type=int, dest='scale'
)




args = arg_parser.parse_args()

dirname = os.path.dirname(__file__)

working_dir = args.Input
OUTPUT_PATH = os.path.join(working_dir,'map.svg')


filename = os.path.join(dirname, 'canvas.svg')
tree = ET.parse(filename)

root = tree.getroot()
d = root.nsmap

main_layer = root.xpath('//*[@id="tiles"]', namespaces={'n': "http://www.w3.org/2000/svg"})[0]

tile_rows = root.xpath('//*[@id="tile_rows"]', namespaces={'n': "http://www.w3.org/2000/svg"})[0]


def deg_coordinates_to_decimal(coordStr):
  coordArr = value.split(', ')
  calculatedCoordArray = []
  for calculation in coordArr:
    calculationArr = calculation.split('/')
    calculatedCoordArray.append(int(calculationArr[0]) / int(calculationArr[1]))
  degrees = calculatedCoordArray[0]
  minutes = calculatedCoordArray[1]
  seconds = calculatedCoordArray[2]
  return (degrees + (minutes * 1/60) + (seconds * 1/60 * 1/60))

# finding the boundaries of the whole canvas
latsArr = []
lonsArr = []

for root_path, directories, file in os.walk(os.path.join(dirname, working_dir)):
  for file in file:
    if(file.endswith(".jpg")):
      print(os.path.join(root_path, file))
      full_filepath = os.path.join(root_path, file)
      with Image(filename=full_filepath) as image:
        print(image.width)
        print(image.height)
        for key, value in image.metadata.items():
          if key == 'exif:GPSLatitude':
            lat = deg_coordinates_to_decimal(value) # lat -> Y vertical
            latsArr.append(lat)
            print("{}: {}".format(key, value))
            print('lat '+ str(lat))
          if key == 'exif:GPSLongitude':
            lon = deg_coordinates_to_decimal(value) # lon -> X horizontal
            lonsArr.append(lon)
            print("{}: {}".format(key, value))
            print('lon '+ str(lon))


minLat = min(latsArr)
minLon = min(lonsArr)
maxLat = max(latsArr)
maxLon = max(lonsArr)
midLon = (minLon + maxLon) /2
midLat = (minLat + maxLat) /2



# find CRS system
utm_crs_list = query_utm_crs_info(
    datum_name="WGS 84",
    area_of_interest=AreaOfInterest(
        west_lon_degree=minLon,
        south_lat_degree=minLat,
        east_lon_degree=maxLon,
        north_lat_degree=maxLat,
    ),
)
utm_crs = CRS.from_epsg(utm_crs_list[0].code)
transformer = Transformer.from_crs("EPSG:4326", utm_crs, always_xy=True)

min_transformed_lon, min_transformed_lat = transformer.transform(minLon, minLat)
max_transformed_lon, max_transformed_lat = transformer.transform(maxLon, maxLat)


width  = max_transformed_lon - min_transformed_lon
height = max_transformed_lat - min_transformed_lat



# def latlngToGlobalXY(lat, lng):
#     earth_radius = 6371
#     # Calculates x based on cos of average of the latitudes
#     x = earth_radius * lng * math.cos((minLat + maxLat)/2)
#     # Calculates y based on latitude
#     y = earth_radius * lat
#     return {x: x, y: y}

# def latlngToScreenXY(lat, lng):
#     topLeft_corner = latlngToGlobalXY(minLat, minLon)
#     bottomRight_corner = latlngToGlobalXY(maxLat, maxLon)
#     # Calculate global X and Y for projection point
#     pos = latlngToGlobalXY(lat, lng)
#     # Calculate the percentage of Global X position in relation to total global width
#     pos.perX = ((pos.x - topLeft_corner.x) / (bottomRight_corner.x - topLeft_corner.x))
#     # Calculate the percentage of Global Y position in relation to total global height
#     pos.perY = ((pos.y - topLeft_corner.y) / (bottomRight_corner.y - topLeft_corner.y))

#     # Returns the screen position based on reference points
#     return {
#         x: p0.scrX + (p1.scrX - p0.scrX)*pos.perX,
#         y: p0.scrY + (p1.scrY - p0.scrY)*pos.perY
#     }






# placing the images into the svg






# image_rotation_up = rotation #32
# image_rotation_down = rotation + 180 #192

for root_path, directories, file in os.walk(os.path.join(dirname, working_dir)):
  for file in file:
    if(file.endswith(".jpg")):
      print(os.path.join(root_path, file))
      full_filepath = os.path.join(root_path, file)
      with Image(filename=full_filepath) as image:
        # print(image.width)
        # print(image.height)
        for key, value in image.metadata.items():
          # print("{}: {}".format(key, value))
          if key == 'exif:GPSLatitude':
            lat = deg_coordinates_to_decimal(value)
            lat_offset = lat - minLat
          if key == 'exif:GPSLongitude':
            lon = deg_coordinates_to_decimal(value)
            lon_offset = lon - minLon
          if key == 'exif:GPSImgDirection':
            direction = value.split('/')
            rotation = ( int(direction[0]) / int(direction[1]) ) / 2  + args.base_rotation
            print('rotation',rotation)
          transformed_lon, transformed_lat = transformer.transform(lon, lat)
          lon_offset = transformed_lon - min_transformed_lon
          lat_offset = transformed_lat - min_transformed_lat
          # print(transformed_lon, min_transformed_lon, transformed_lat, min_transformed_lat)
          # print('lon_offset, lat_offset', lon_offset, lat_offset)
        g_pos_el_attributes = {
          'transform': "translate({}, {})".format(format(lon_offset*args.scale, '.20f'), format(lat_offset*args.scale*-1, '.20f')),
          'data-lat-offset': format(lat_offset, '.20f'),
          'data-lon-offset': format(lon_offset, '.20f'),
          'class': 'tile',
          'id': 'tile_{}'.format(file.split('.')[0]),
        }
        g_pos_el = ET.SubElement(main_layer, 'g', attrib=g_pos_el_attributes)

        g_offset_corr_el_attributes = {
          'transform': "translate({}, {})".format(-image.width/2, -image.height/2),
          'class': 'tile-offset-corr',
        }
        g_offset_corr_el = ET.SubElement(g_pos_el, 'g', attrib=g_offset_corr_el_attributes)

        g_rot_el_attributes = {
          'class': 'tile-rotate',
          'data-image-rotation': str(rotation),
          'data-image-dimensions': str(image.width) + ' ' + str(image.height),
          'transform': 'rotate({} {} {})'.format(str(rotation), str(image.width/2), str(image.height/2))
          # 'transform': 'rotate({} {} {})'.format(str(rotation), 0,0)
        }
        g_rot_el = ET.SubElement(g_offset_corr_el, 'g', attrib=g_rot_el_attributes)

        xlinkns ="http://www.w3.org/1999/xlink"
        image_el = ET.SubElement(g_rot_el, 'image', {
          "class": 'thermal_image',
          "{%s}href" % xlinkns: file,
          "width": str(image.width),
          "height": str(image.height),
          'data-lat': format(lat, '.20f'),
          'data-lon': format(lon, '.20f'),
        })


# sort elements
def getkey(elem):
    # Used for sorting elements by @LIN.
    # returns a tuple of ints from the exploded @LIN value
    # '1.0' -> (1,0)
    # '1.0.1' -> (1,0,1)
    return float(elem.get('id').split('_')[2])

main_layer[:] = sorted(main_layer, key=getkey)


# find rows
# up/down is actually left/right or right/left
last_state = 'down'
for index, el in enumerate(main_layer):
  if(el.getprevious() is not None):
    if (el.getprevious().attrib['data-lon-offset'] > el.attrib['data-lon-offset'] or (el.getprevious().attrib['data-lon-offset'] == el.attrib['data-lon-offset'] and last_state == 'up')):
      print('up')
      rot_el = el[0][0]
      el.attrib['data-direction'] = 'up'
    else:
      rot_el = el[0][0]
      el.attrib['data-direction'] = 'down'
      print('down')

    # NOT NEEDED SINCE THERE IS A ROTATION INFORMATION
    # merge tiles into groups
    # print(index)
    # print("el.attrib['data-direction'] " + el.attrib['data-direction'])
    # print("last_state " + last_state)
    if index is 1 or last_state is not el.attrib['data-direction']:
      current_row = ET.SubElement(tile_rows, 'g', attrib={ 'class': 'tile-row' })
    copyElem = copy.deepcopy(el)
    current_row.insert(0, copyElem)
    last_state = el.attrib['data-direction']

# remove temporary group
root.remove(main_layer)

# resize canvas to tiles and add some padding


print(width, height, args.scale)
scaled_width = width * args.scale
scaled_height = height * args.scale

padding = 500

canvas_width  = str(scaled_width  + padding*2)
canvas_height = str(scaled_height + padding*2)

viewbox_x = str(padding * -1)
viewbox_y = str((scaled_height + padding) * -1)
viewbox_width = canvas_width
viewbox_height = canvas_height

root.attrib['width'] = canvas_width
root.attrib['height'] = canvas_height
root.attrib['viewBox'] = "{} {} {} {}".format(viewbox_x, viewbox_y, viewbox_width, viewbox_height)



# Finally save the svg
with open(OUTPUT_PATH, 'wb') as f:
    tree.write(f, encoding='utf-8')


print('Done!')