deep-crop/convert-svg-to-geotiff.py

import os
import argparse
import lxml.etree as ET
import subprocess
import flirimageextractor
import cv2
import numpy as np
from pathlib import Path
from selenium import webdriver
import rasterio
# import shapefile

arg_parser = argparse.ArgumentParser(description='Export SVG composition of FLIR images as TIFF with thermo layer')

arg_parser.add_argument('Input',
                       metavar='input_svg',
                       type=str,
                       help='Path to the input SVG file cotaining xlinks to FLIR images')

arg_parser.add_argument('Output',
                       metavar='output_tiff',
                       type=str,
                       help='Output filename')


args = arg_parser.parse_args()
dirname = os.path.dirname(__file__)
INPUT_PATH = os.path.join(dirname, args.Input)
INPUT_DIR = os.path.split(INPUT_PATH)[0]

TEMP_MAP_THERMALPNG_SVG_PATH = os.path.join(INPUT_DIR, 'map_thermalpng.svg')
TEMP_MAP_THERMALPNG_PATH = os.path.join(INPUT_DIR, 'map_thermalpng.png')
TEMP_MAP_PREVIEW_PATH = os.path.join(INPUT_DIR, 'map_preview.png')
TEMP_MAP_UNGROUPED_SVG_PATH = os.path.join(INPUT_DIR, 'map_ungrouped.svg')

TEMP_MAP_ALIGNMENTPROOF_PATH = os.path.join(INPUT_DIR, 'map_thermalpng_proof.png')


THERMALPNG_DIR = 'thermalpngs'

# VECTOR_SHAPEFILE_PATH = os.path.join(INPUT_DIR, 'shapefile')


OUTPUT_PATH = os.path.join(dirname, args.Output)



def make_thermalpng_tiles():
  """
  Extract thermal infomration as greyscale PNG-16 
  (temp * 1000 to retain some decimals)
  and save the png tiles.
  Assuming, that the data will always have a positive value.
  """
  print('Building PNG tiles representing the thermal data ...')
  Path(os.path.join(INPUT_DIR, THERMALPNG_DIR)).mkdir(parents=True, exist_ok=True)
  png_output_dir = os.path.join(INPUT_DIR, THERMALPNG_DIR)

  done_files_count = 0
  for root_path, directories, file in os.walk(os.path.join(dirname, INPUT_DIR)):
    for file in file:
      if(file.endswith(".jpg")):
        output_file_path = os.path.join(png_output_dir, file + '.thermal.png')
        if os.path.isfile(output_file_path):
          done_files_count = done_files_count+1
        else:
          print('  Processing ' + file)
          full_filepath = os.path.join(root_path, file)
          flir = flirimageextractor.FlirImageExtractor()
          flir.process_image(full_filepath)
          thermal_img_np = flir.thermal_image_np
          multiplied_image = cv2.multiply(thermal_img_np, 1000)
          cv2.imwrite(output_file_path, multiplied_image.astype(np.uint16))
  if done_files_count != 0:
    print('Using {} pre-built tiles.'.format(str(done_files_count)))




def make_thermalpng_svg():
  """
  replaces the image paths with the thermal pngs
  and creates new SVG file
  """
  print('Replacing the images inside the SVG with the PNG tiles ...')
  # print("svg_file")
  # print(dir(svg_file))
  tree = ET.parse(INPUT_PATH)
  root = tree.getroot()
  # print(ET.tostring(root))
  # tile_rows = root.xpath('//image', namespaces={'n': "http://www.w3.org/2000/svg"})
  # print(dir(root))
  tile_elements = root.xpath('//*[@class="thermal_image"]')
  linkattrib ='{http://www.w3.org/1999/xlink}href'
  for tile in tile_elements:
    new_file_path = os.path.join(THERMALPNG_DIR, tile.attrib[linkattrib] + '.thermal.png')
    tile.attrib[linkattrib] = new_file_path
    # Post Production
    # tile.attrib["mask"] = 'url(#tilefademask)'
    # tile.attrib["style"] = 'opacity:.7'


  # newxml = ET.tostring(tree, encoding="unicode")
  # print(newxml)
  # return newxml




  with open(TEMP_MAP_THERMALPNG_SVG_PATH, 'wb') as f:
    tree.write(f, encoding='utf-8')

  return tree



def make_thermalpng():
  """
  exports the SVG canvas as Gray_16 PNG
  """
  print('Creating a big 16-bit grayscale PNG image representing the thermal data out of the SVG file ...')
  command = [
              'inkscape',
              '--pipe',
              '--export-type=png',
              '--export-png-color-mode=Gray_16'
            ],
  input_file = open(TEMP_MAP_THERMALPNG_SVG_PATH, "rb")
  output_file = open(TEMP_MAP_THERMALPNG_PATH, "wb")
  completed = subprocess.run(
        *command,
        cwd=INPUT_DIR, # needed for reative image links
        stdin=input_file,
        stdout=output_file
      )
  return completed

# def make_thermalpreview():
#   """
#   exports the preview image
#   """
#   command = [
#               '/snap/bin/inkscape',
#               '--pipe',
#               '--export-type=png',
#               '--export-png-color-mode=Gray_8'
#             ],
#   input_file = open(TEMP_MAP_THERMALPNG_SVG_PATH, "rb")
#   output_file = open(TEMP_MAP_PREVIEW_PATH, "wb")
#   completed = subprocess.run(
#         *command,
#         cwd=INPUT_DIR, # needed for reative image links
#         stdin=input_file,
#         stdout=output_file
#       )
#   return completed



def get_thermal_numpy_array():
  print('Converting the PNG into NumPy Array and normalize temperature values ...')
  image = cv2.imread(TEMP_MAP_THERMALPNG_PATH, cv2.IMREAD_ANYDEPTH)
  image_float = image.astype(np.float32)
  image_float_normalized = cv2.divide(image_float, 1000)
  # print(image_float_normalized[1000][905]) # looking what's the value of some pixel
  return image_float_normalized

def get_used_tiles_relpaths():
  """
  outputs an array of all used tile filenames in the input SVG 
  (relative filepaths like they appear in the svg.)
  """
  images = []
  tree = ET.parse(INPUT_PATH)
  root = tree.getroot()
  tile_elements = root.xpath('//*[@class="thermal_image"]')
  linkattrib ='{http://www.w3.org/1999/xlink}href'
  for tile in tile_elements:
    images.append(tile.attrib[linkattrib])
  return images



# def deg_coordinates_to_decimal(coordStr):
#   coordArr = coordStr.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]
#   decimal = (degrees + (minutes * 1/60) + (seconds * 1/60 * 1/60))
#   # print(decimal)
#   return decimal


# def read_coordinates_from_tile(filename):
#   full_filepath = os.path.join(INPUT_DIR, filename)
#   with Image(filename=full_filepath) as image:
#     for key, value in image.metadata.items():
#       if key == 'exif:GPSLatitude':
#         # print('latstr', value)
#         lat = deg_coordinates_to_decimal(value) # lat -> Y vertical
#       if key == 'exif:GPSLongitude':
#         # print('lonstr', value)
#         lon = deg_coordinates_to_decimal(value) # lon -> X horizontal

#   return [lat, lon]


# def get_coordinate_boundaries():
#   image_names = get_used_tiles_relpaths()
#   coordinates = {
#     'lat': [],
#     'lon': []
#   }
#   for filename in image_names:
#     tile_coordinates = read_coordinates_from_tile(filename)
#     coordinates['lat'].append(tile_coordinates[0])
#     coordinates['lon'].append(tile_coordinates[1])

#   boundaries = {
#     'xmin': min(coordinates['lon']),
#     'xmax': max(coordinates['lon']),
#     'ymin': min(coordinates['lat']),
#     'ymax': max(coordinates['lat']),
#   }
#   return boundaries 


# def create_ungrouped_svg():
#   """
#   exports the SVG without any grouped elements 
#   (the quick and dirty way)
#   """
#   print('Create an SVG without groups ...')
#   inkscape_actions = "select-all:groups; SelectionUnGroup; select-all:groups; SelectionUnGroup; select-all:groups; SelectionUnGroup; select-all:groups; SelectionUnGroup; select-all:groups; SelectionUnGroup; select-all:groups; SelectionUnGroup; select-all:groups; SelectionUnGroup; export-filename: {}; export-plain-svg; export-do;".format(TEMP_MAP_UNGROUPED_SVG_PATH)
#   command = [
#               '/snap/bin/inkscape',
#               '--pipe',
#               '--actions={}'.format(inkscape_actions),
#             ],
#   input_file = open(INPUT_PATH, "rb")
#   completed = subprocess.run(
#         *command,
#         cwd=INPUT_DIR, # needed for reative image links
#         stdin=input_file,
#       )
#   print('completed', completed)
#   return completed


def get_ground_control_points():
  """
  Using selenium with firefox for rendering the SVG and 
  getting the positional matching between GPS 
  and x/y coords of SVG-image.
  image tags need to have the gps data attached as data attributes.
  """
  print('Getting ground control points ...')
  options = webdriver.firefox.options.Options()
  options.headless = True
  driver = webdriver.Firefox(options=options)
  driver.get("file://{}".format(INPUT_PATH))
  images = driver.find_elements_by_class_name("thermal_image")

  gcps = []
  for image in images:
    location = image.location
    size = image.size
    raster_y = float(location['y'] + size['height']/2)
    raster_x = float(location['x'] + size['width']/2)
    reference_lon = float(image.get_attribute('data-lon'))
    reference_lat = float(image.get_attribute('data-lat'))
    imageMapping = rasterio.control.GroundControlPoint(row=raster_y, col=raster_x, x=reference_lon, y=reference_lat)
    gcps.append(imageMapping)

  driver.quit()
  return gcps


# def make_vector_shapefile():
#   w = shapefile.Writer(VECTOR_SHAPEFILE_PATH)
#   w.field('name', 'C')

#   tree = ET.parse(INPUT_PATH)
#   root = tree.getroot()
#   tiles = root.xpath('//*[@class="thermal_image"]')
#   for index, tile in enumerate(tiles):
#     w.point(float(tile.attrib['data-lon']), float(tile.attrib['data-lat']))
#     w.record('point{}'.format(index))
#   w.close()
#   return True

def verify_coordinate_matching():
  """
  During development, i want to proof 
  that the point/coordinate matching is right.
  Producing a png file with red dots for visual proof.
  """
  img = cv2.imread(TEMP_MAP_THERMALPNG_PATH, cv2.IMREAD_GRAYSCALE)
  img = cv2.cvtColor(img,cv2.COLOR_GRAY2RGB)

  options = webdriver.firefox.options.Options()
  # options.headless = True
  driver = webdriver.Firefox(options=options)
  driver.get("file://{}".format(INPUT_PATH))
  images = driver.find_elements_by_class_name("thermal_image")

  gcps = []
  for image in images:
    location = image.location
    size = image.size
    raster_y = int(location['y'] + size['height']/2)
    raster_x = int(location['x'] + size['width']/2)
    reference_lon = float(image.get_attribute('data-lon'))
    reference_lat = float(image.get_attribute('data-lat'))
    print(raster_x, raster_y)
    img = cv2.circle(img, (raster_x, raster_y), radius=10, color=(0, 0, 255), thickness=-1)
  driver.quit()
  cv2.imwrite(TEMP_MAP_ALIGNMENTPROOF_PATH, img)




def make_geotiff_image():

  thermal_numpy_array = get_thermal_numpy_array()
  thermal_numpy_array[thermal_numpy_array == 0] = np.NaN # zeros to NaN

  # # coordinates of all tiles
  # geo_bound = get_coordinate_boundaries()
  # print('boundaries', geo_bound)

  np_shape = thermal_numpy_array.shape
  image_size = (np_shape[0], np_shape[1])

  gcps = get_ground_control_points()
  print('Applying affine transform ...')
  gcp_transform = rasterio.transform.from_gcps(gcps)
  print(gcp_transform)

  print('Generating the GeoTiff ...')

  raster_io_dataset = rasterio.open(
      OUTPUT_PATH,
      'w',
      driver='GTiff',
      height=thermal_numpy_array.shape[0],
      width=thermal_numpy_array.shape[1],
      count=1,
      dtype=thermal_numpy_array.dtype,
      transform=gcp_transform,
      crs='+proj=latlong',
  )
  raster_io_dataset.write(thermal_numpy_array, 1)

  # # try to get rid of the black frame 
  # src = rasterio.open(OUTPUT_PATH)
  # src[src == 0] = np.NaN # zeros to NaN
  # raster_io_dataset2 = rasterio.open(
  #     OUTPUT_PATH,
  #     'w',
  #     driver='GTiff',
  #     height=src.shape[0],
  #     width=src.shape[1],
  #     count=1,
  #     dtype=src.dtype,
  #     crs='+proj=latlong',
  # )
  # raster_io_dataset2.write(src, 1)

  print('Saved to ', OUTPUT_PATH)



make_thermalpng_tiles()
make_thermalpng_svg()
make_thermalpng()
make_geotiff_image()

# # Helpers for debugging 
# verify_coordinate_matching()
# make_vector_shapefile()