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