r_interp_opc.py 6.41 KB
#!/usr/bin/env python
# -*- coding: utf-8 -*-

import os
import glob
import scipy.io
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
import datetime as dt
import time
from math import radians, cos, sin, asin, sqrt
from mpl_toolkits.basemap import Basemap
import subprocess


def haversine(point1, point2, miles=False):
    """ Calculate the great-circle distance bewteen two points on the Earth surface.

    :input: two 2-tuples, containing the latitude and longitude of each point
    in decimal degrees.

    Example: haversine((45.7597, 4.8422), (48.8567, 2.3508))

    :output: Returns the distance bewteen the two points.
    The default unit is kilometers. Miles can be returned
    if the ``miles`` parameter is set to True.

    """
    # average earth radius
    avg_earth_rad = 6371.2  # km 

    # unpack latitude/longitude
    lat1, lng1 = point1
    lat2, lng2 = point2

    # convert all latitudes/longitudes from decimal degrees to radians
    lat1, lng1, lat2, lng2 = map(radians, (lat1, lng1, lat2, lng2))

    # calculate haversine
    lat = lat2 - lat1
    lng = lng2 - lng1
    d = sin(lat * 0.5) ** 2 + cos(lat1) * cos(lat2) * sin(lng * 0.5) ** 2
    h = 2 * avg_earth_rad * asin(sqrt(d))
    if miles:
        return h * 0.621371  # in miles
    else:
        return h  # in kilometers


def figResolution(resol, lon1, lat1, lon2, lat2):
    distW = haversine((lat1, lon1), (lat1, lon2), miles=False) * 1000
    pixW = distW / resol
    distH = haversine((lat1, lon1), (lat2, lon1), miles=False) * 1000
    pixH = distH / resol
    return pixW, pixH


def writeReadMe():
    """
    write a ReadMe File with georectification informations
    """
    f = open('README.txt','w')
    f.write('README\n')
    f.write('This file contains informations about the georectification performed\n ')
    f.write('with the python script r_interp_opc\n')
    f.write('----------------------\n')
    f.write('Date of processing : ' + dt.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + '\n')
    f.write('Lower Left corner : ' + str(lat1) + '/' + str(lon1) + '\n')
    f.write('Upper Right corner : ' + str(lat2) + ' / ' + str(lon2) + '\n')
    f.write('Approximate resolution in x and y : ' + str(resolution) + ' m/pixel \n')
    f.write('First georectified image : ' + glob.glob('*.jpg')[0] + '\n')
    f.write('Last georectified image : ' + glob.glob('*.jpg')[-1] + '\n')
    f.write('----------------------' + '\n')
    f.write('To produce the geotiff images from the png, use this bash command line : \n')
    f.write('gdal_translate -of Gtiff -a_ullr lon2 lat1 lon2 -a_srs "WGS84" filein fileout \n')
    f.write('*** upper right corner / lower left corner coordinates *** \n')
    f.write('----------------------\n')
    f.write('Please note that the resolution may vary of 10% of the requested value\n')
    f.write('The best way to be sure of the precise resolution is to check\n')
    f.write('on a GIS plateform with the tif images.\n')
    f.close()


def georect_opc(file, m_file):
    """
    Georectification des images du Pic Champlain. 
    choisir le bon fichier de georect matlab (correction de la distorsion)
    :input: file, georectified matlab file
    :output: retruns a png of the georectified area
    """

    filename = os.path.basename(file)
    datef = filename[:13]
    date_img = dt.datetime.strptime(datef, '%Y%m%d_%H%M')

    # read image
    img = mpimg.imread(file)

    # B&W
    if img.ndim == 3:
        img_mean = np.transpose(img.mean(axis=(2)))
    else:
        img_mean = np.transpose(img)

    # read LON & LAT from mfile
    M_data = scipy.io.loadmat(m_file)
    LON = M_data['LON']
    LAT = M_data['LAT']

    # masking upper and lower part of image
    mask = np.ones((LON.shape))
    # mask[:] = 1

    # apply mask over image and LAT/LON
    I = img_mean * mask
    lat = LAT * mask
    lon = LON * mask

    # figure
    fig = plt.figure(frameon=False)
    # set
    a = fig.gca()
    a.set_frame_on(False)

    # define the map settings
    lat_mean = np.mean([lat1, lat2])
    lon_mean = np.mean([lon1, lon2])
    # ~ +proj=lcc +lat_1=60 +lat_2=46 +lat_0=44 +lon_0=-68.5 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs
    map = Basemap(projection='tmerc',
                  lat_0=lat_mean, lon_0=lon_mean,
                  llcrnrlat=lat1, urcrnrlat=lat2,
                  llcrnrlon=lon1, urcrnrlon=lon2,
                  resolution='c', fix_aspect=True)
    # ~ map = Basemap(projection='stere',
    # ~ lon_0=lon_mean, lat_0=lat_mean,
    # ~ llcrnrlat=lat1, urcrnrlat=lat2,
    # ~ llcrnrlon=lon1, urcrnrlon=lon2,
    # ~ rsphere=(avg_earth_rad), resolution='c',
    # ~ fix_aspect=True)

    # plot the picture over the ground
    map.pcolormesh(lon, lat, I, latlon=True,
                   cmap=plt.cm.Greys_r,
                   vmin=np.nanmin(I),
                   vmax=np.nanmax(I))

    # add date / time
    # ~ plt.annotate(str(date_img) + ' UTC',
    # ~ color="cyan", size=16,
    # ~ xy=(1970, 460), xycoords='axes pixels')

    # save figure
    w, h = figResolution(resolution, lon1, lat1, lon2, lat2)
    dpi = 100
    fig.set_size_inches(w / dpi, h / dpi)
    fig.tight_layout()
    plt.savefig('opc_' + datef + '.png', transparent=True,
                bbox_inches='tight', pad_inches=0, dpi=dpi)
    plt.close()


def getgeotiff():
    filein = 'opc_' + datef + '.png'
    fileout = 'opc_' + datef + '.tif'
    exec_gdal = 'gdal_translate -of Gtiff -a_ullr ' + str(lon2) \
                + ' ' + str(lat1) + ' ' + str(lon2) \
                + ' ' + str(lat1) + ' -a_srs "WGS84" ' + filein + ' ' + fileout
    subprocess.Popen(exec_gdal)


if __name__ == '__main__':

    # ~ #root = '/share/work/nicp0003/opc/work/todo/'
    # ~ root = '/home/nicopa01/data/opc/2017/test_res'
    # ~ mfile = '/home/nicopa01/data/opc/2017/image/20170217/opc_20170217_dist.mat'

    root = raw_input('Folder to process: ')
    mfile = raw_input('georect matrix mfile: ')
    lat1 = float(raw_input('latitude of the lower left corner: '))
    lon1 = float(raw_input('longitude of the lower left corner : '))
    lat2 = float(raw_input('latitude of the upper right corner : '))
    lon2 = float(raw_input('longitude of the upper right corner : '))
    resolution = float(raw_input('pixel resolution (in meter): '))

    os.chdir(root)
    for f in os.listdir(root):
        if f.endswith(".jpg"):
            georect_opc(f, mfile)
            print os.path.splitext(f)[0] + ' processed at ' + dt.datetime.now().strftime("%H:%M:%S")
    writeReadMe()