deep-crop/export-gis-jpg.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 wand.image import Image
from osgeo import gdal
from osgeo import osr

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')
THERMALPNG_DIR = 'thermalpngs'

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
  """
  Path(os.path.join(INPUT_DIR, THERMALPNG_DIR)).mkdir(parents=True, exist_ok=True)
  png_output_dir = os.path.join(INPUT_DIR, THERMALPNG_DIR)

  for root_path, directories, file in os.walk(os.path.join(dirname, INPUT_DIR)):
    for file in file:
      if(file.endswith(".jpg")):
        print('Extracting thermal info from ' + 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)
        output_file_path = os.path.join(png_output_dir, file + '.thermal.png')
        print(output_file_path)
        cv2.imwrite(output_file_path, multiplied_image.astype(np.uint16))




def make_thermalpng_svg():
  """
  replaces the image paths with the thermal pngs
  and creates new SVG file
  """
  # 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:
    tile.attrib[linkattrib] = os.path.join(THERMALPNG_DIR, tile.attrib[linkattrib] + '.thermal.png')
  # 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
  """
  command = [
              '/snap/bin/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 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():
  # input_file = open(TEMP_MAP_THERMALPNG_PATH, "rb")
  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])
  # cv2.imshow("OpenCV Image Reading", image)
  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
      if key == 'exif:GPSImgDirection':
        direction = value.split('/')
        print(int(direction[0])/int(direction[1])/2, '    ', (value))
        
  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 make_geotiff_image():
  thermal_numpy_array = get_thermal_numpy_array()

  # 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])

  # set geotransform
  nx = image_size[0]
  ny = image_size[1]
  xres = (geo_bound['xmax'] - geo_bound['xmin']) / float(nx)
  yres = (geo_bound['ymax'] - geo_bound['ymin']) / float(ny)
  geotransform = (geo_bound['xmin'], xres, 0, geo_bound['ymax'], 0, -yres)

  # create the 3-band raster file
  dst_ds = gdal.GetDriverByName('GTiff').Create(OUTPUT_PATH, ny, nx, 1, gdal.GDT_Float32)
  dst_ds.SetGeoTransform(geotransform)    # specify coords
  srs = osr.SpatialReference()            # establish encoding
  res = srs.SetWellKnownGeogCS( "WGS84" ) # WGS84 lat/long
  dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file
  dst_ds.GetRasterBand(1).WriteArray(thermal_numpy_array)   # write thermal-band to the raster
  dst_ds.FlushCache()                     # write to disk




# make_thermalpng_tiles()
# make_thermalpng_svg()
# make_thermalpreview()
# make_thermalpng()
make_geotiff_image()


# dataset = gdal.Open("working_result_example.tif", gdal.GA_ReadOnly)
# print(dir(dataset))
# print(dataset.GetMetadata_List())