contline.py

   1 # following https://stackoverflow.com/questions/18304722/python-find-contour-lines-from-matplotlib-pyplot-contour
   2 import matplotlib.pyplot as plt
   3 import numpy as np
   4 from scipy.io import loadmat
   5 from utils import *
   6 import datetime
   7 import time
   8 from JPSSD import time_num2iso
   9 from contour2kml import contour2kml
  10 from scipy.ndimage import gaussian_filter
  11 import os
  12 import sys
  13
  14 def get_contour_verts(xx, yy, zz, time_num_granules, contour_dt_hours=6, contour_dt_init=6, contour_dt_final=6, gauss_filter=True, plot_contours=False):
  15     # Computing the levels
  16     # Datetimes for the first and last granule
  17     dt1=datetime.datetime.fromtimestamp(time_num_granules[0])
  18     dt2=datetime.datetime.fromtimestamp(time_num_granules[-1])
  19     # Starting 6 hours before the first granule and finishing 6 hours after the last granule
  20     dti=dt1-datetime.timedelta(hours=contour_dt_init)
  21     dtf=dt2+datetime.timedelta(hours=contour_dt_final)
  22     # Number of periods of 6 hours we need to compute rounded
  23     d=dtf-dti
  24     hours=d.total_seconds()/3600
  25     M=int(np.round((hours+1)/contour_dt_hours ))
  26     # Datetimes where we are going to compute the levels
  27     dts=[dti+datetime.timedelta(hours=k*contour_dt_hours) for k in range(1,M)]
  28     levels=[time.mktime(t.timetuple()) for t in dts]
  29
  30     # Scaling the time components as in detections
  31     time_num=np.array(levels)
  32     time_iso=[time_num2iso(t) for t in time_num]
  33
  34     # Computing and generating the coordinates for the contours
  35     contours = []
  36     if gauss_filter:
  37         # for each level
  38         for level in levels:
  39             # copy the fire arrival time
  40             Z = np.array(zz)
  41             # distinguish either in or out the perimeter
  42             Z[Z < level] = 0
  43             Z[Z >= level] = 1
  44             # smooth the perimeter using a gaussian filter
  45             sigma = 2.
  46             ZZ = gaussian_filter(Z,sigma)
  47             # find the contour line in between
  48             cn = plt.contour(xx,yy,ZZ,levels=0.5)
  49             # contour line
  50             cc = cn.collections[0]
  51             # initialize the path
  52             paths = []
  53             # for each separate section of the contour line
  54             for pp in cc.get_paths():
  55                 # read all the vertices
  56                 paths.append(pp.vertices)
  57             contours.append(paths)
  58     else:
  59         # computing all the contours
  60         cn = plt.contour(xx,yy,zz,levels=levels)
  61         # for each contour line
  62         for cc in cn.collections:
  63             # initialize the path
  64             paths = []
  65             # for each separate section of the contour line
  66             for pp in cc.get_paths():
  67                 # read all the vertices
  68                 paths.append(pp.vertices)
  69             contours.append(paths)
  70
  71     # Plotting or not the contour lines
  72     if plot_contours:
  73         print 'contours are collections of line, each line consisting of poins with x and y coordinates'
  74         for c in contours:
  75             for cc in c:
  76                 xx=[x[0] for x in cc]
  77                 yy=[y[1] for y in cc]
  78                 plt.scatter(xx,yy)
  79         plt.show()
  80
  81     # Creating an array of dictionaries for each perimeter
  82     conts=[Dict({'text':time_iso[k],
  83         'LineStyle':{
  84             'color':'ff081388',
  85             'width':'2.5',
  86         },
  87         'PolyStyle':{
  88             'color':'66000086',
  89             'colorMode':'random'
  90         },
  91         'time_begin':time_iso[k],
  92         'polygons': contours[k] }) for k in range(0,len(contours))]
  93
  94     # Creating a dictionary to store the KML file information
  95     data=Dict({'name':'contours.kml',
  96         'folder_name':'Perimeters'})
  97     data.update({'contours': conts})
  98
  99     return data
 100
 101 if __name__ == "__main__":
 102     result_file = 'result.mat'
 103     mgout_file = 'mgout.mat'
 104     if os.path.isfile(result_file) and os.access(result_file,os.R_OK) and os.path.isfile(mgout_file) and os.access(mgout_file,os.R_OK):
 105         print 'Loading the data...'
 106         result=loadmat(result_file)
 107         mgout=loadmat(mgout_file)
 108     else:
 109         print 'Error: file %s or %s not exist or not readable' % (result_file,mgout_file)
 110         sys.exit(1)
 111
 112     # Indexing the coordinates into the same as the multigrid solution
 113     xind=mgout['sm'][0]-1
 114     yind=mgout['sn'][0]-1
 115     x=np.array(result['fxlon'])
 116     xx=x[np.ix_(xind,yind)]
 117     y=np.array(result['fxlat'])
 118     yy=y[np.ix_(xind,yind)]
 119     tscale=mgout['tscale'][0]
 120     time_scale_num=mgout['time_scale_num'][0]
 121     zz=mgout['a']*tscale+time_scale_num[0]
 122
 123     print 'Computing the contours...'
 124     # Granules numeric times
 125     time_num_granules = result['time_num_granules'][0]
 126     data = get_contour_verts(xx, yy, zz, time_num_granules, contour_dt_hours=24, gauss_filter=False)
 127
 128     print 'Creating the KML file...'
 129     # Creating the KML file
 130     contour2kml(data,'perimeters.kml')
 131
 132     print 'perimeters.kml generated'
 133