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.

"""

# 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

"""
write a ReadMe File with georectification informations
"""
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')

# 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
LON = M_data['LON']
LAT = M_data['LAT']

# masking upper and lower part of image

# apply mask over image and LAT/LON

# 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,
# ~ 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))

# ~ 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,
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")