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 clean_polys(paths,bounds,plot=True):
  13     npaths = []
  14     maxd = .05
  15     eps = 1e-5
  16     dx = abs(bounds[1]-bounds[0])
  17     dy = abs(bounds[3]-bounds[2])
  18     xb = np.array(bounds[:2])
  19     yb = np.array(bounds[2:])
  20     for p in paths:
  21         diff = np.abs(p[-1]-p[0])
  22         if (diff>maxd).all() or (abs(diff[0]-dx)<eps or abs(diff[1]-dy)<eps):
  23             x = [c[0] for c in p]
  24             y = [c[1] for c in p]
  25             xm = (x[0]+x[-1])/2
  26             ym = (y[0]+y[-1])/2
  27             xn = x[:]
  28             yn = y[:]
  29             x.append(p[0][0])
  30             y.append(p[0][1])
  31
  32             if abs(diff[0]-dx)<eps:
  33                 ex = xn[-1]
  34                 ey = yb[np.argmin(abs(yb-xn[-1]))]
  35                 p = np.append(p,np.array([[ex,ey]]),axis=0)
  36                 ex = xn[0]
  37                 p = np.append(p,np.array([[ex,ey]]),axis=0)
  38                 #p = np.append(p,np.array([p[0]]),axis=0)
  39                 if False:
  40                     plt.ion()
  41                     plt.plot([xb[0],xb[0],xb[1],xb[1],xb[0]],[yb[0],yb[1],yb[1],yb[0],yb[0]],'k')
  42                     plt.plot(x,y,'r')
  43                     x = [c[0] for c in p]
  44                     y = [c[1] for c in p]
  45                     plt.plot(x,y,'g')
  46                     plt.xlim(xb)
  47                     plt.ylim(yb)
  48                     plt.show()
  49                     plt.pause(.001)
  50                     plt.cla()
  51             elif abs(diff[1]-dy)<eps:
  52                 #ex = xb[np.argmin(abs(xb-xn[-1]))]
  53                 ex = xb[0]
  54                 ey = yn[-1]
  55                 p = np.append(p,np.array([[ex,ey]]),axis=0)
  56                 ey = yn[0]
  57                 p = np.append(p,np.array([[ex,ey]]),axis=0)
  58                 #p = np.append(p,np.array([p[0]]),axis=0)
  59                 if plot:
  60                     #plt.ion()
  61                     plt.plot([xb[0],xb[0],xb[1],xb[1],xb[0]],[yb[0],yb[1],yb[1],yb[0],yb[0]],'k')
  62                     plt.plot(x,y,'r')
  63                     x = [c[0] for c in p]
  64                     y = [c[1] for c in p]
  65                     plt.plot(x,y,'g')
  66                     plt.xlim(xb)
  67                     plt.ylim(yb)
  68                     plt.show()
  69                     #plt.pause(.001)
  70                     #plt.cla()
  71             else:
  72                 ex = xb[np.argmin(abs(xm-xb))]
  73                 ey = yb[np.argmin(abs(ym-yb))]
  74                 p = np.append(p,np.array([[ex,ey]]),axis=0)
  75                 #p = np.append(p,np.array([p[0]]),axis=0)
  76                 if False:
  77                     plt.ion()
  78                     plt.plot([xb[0],xb[0],xb[1],xb[1],xb[0]],[yb[0],yb[1],yb[1],yb[0],yb[0]],'k')
  79                     plt.plot(x,y,'r')
  80                     x = [c[0] for c in p]
  81                     y = [c[1] for c in p]
  82                     plt.plot(xn,yn,'g')
  83                     plt.xlim(xb)
  84                     plt.ylim(yb)
  85                     plt.show()
  86                     plt.pause(.001)
  87                     plt.cla()
  88
  89         else:
  90             npaths.append(p)
  91     return np.array(npaths)
  92
  93 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):
  94     bounds = (xx.min(),xx.max(),yy.min(),yy.max())
  95     fig = plt.figure()
  96     # Computing the levels
  97     # Datetimes for the first and last granule
  98     dt1=datetime.datetime.fromtimestamp(time_num_granules[0])
  99     dt2=datetime.datetime.fromtimestamp(time_num_granules[-1])
 100     # Starting 6 hours before the first granule and finishing 6 hours after the last granule
 101     dti=dt1-datetime.timedelta(hours=contour_dt_init)
 102     dtf=dt2+datetime.timedelta(hours=contour_dt_final)
 103     # Number of periods of 6 hours we need to compute rounded
 104     d=dtf-dti
 105     hours=d.total_seconds()/3600
 106     M=int(np.round((hours+1)/contour_dt_hours ))
 107     # Datetimes where we are going to compute the levels
 108     dts=[dti+datetime.timedelta(hours=k*contour_dt_hours) for k in range(1,M)]
 109     if levels_gran:
 110         levels=time_num_granules
 111     else:
 112         levels=[time.mktime(t.timetuple()) for t in dts]
 113
 114     # Scaling the time components as in detections
 115     time_num=np.array(levels)
 116     time_iso=[time_num2iso(t) for t in time_num]
 117
 118     # Computing and generating the coordinates for the contours
 119     contours = []
 120     if gauss_filter:
 121         # for each level
 122         for level in levels:
 123             # copy the fire arrival time
 124             Z = np.array(zz)
 125             # distinguish either in or out the perimeter
 126             Z[Z < level] = 0
 127             Z[Z >= level] = 1
 128             # smooth the perimeter using a gaussian filter
 129             sigma = 2.
 130             ZZ = gaussian_filter(Z,sigma)
 131             # find the contour line in between
 132             cn = plt.contour(xx,yy,ZZ,levels=0.5)
 133             # contour line
 134             cc = cn.collections[0]
 135             # initialize the path
 136             paths = []
 137             # for each separate section of the contour line
 138             for pp in cc.get_paths():
 139                 # read all the vertices
 140                 paths.append(pp.vertices)
 141             contours.append(paths)
 142             #contours.append(clean_polys(paths,bounds))
 143     else:
 144         # computing all the contours
 145         cn = plt.contour(xx,yy,zz,levels=levels)
 146         # for each contour line
 147         for cc in cn.collections:
 148             # initialize the path
 149             paths = []
 150             # for each separate section of the contour line
 151             for pp in cc.get_paths():
 152                 # read all the vertices
 153                 paths.append(pp.vertices)
 154             contours.append(paths)
 155
 156     # Plotting or not the contour lines
 157     if plot_contours:
 158         print 'contours are collections of line, each line consisting of poins with x and y coordinates'
 159         for c in contours:
 160             for cc in c:
 161                 xx=[x[0] for x in cc]
 162                 yy=[y[1] for y in cc]
 163                 plt.scatter(xx,yy)
 164         plt.show()
 165
 166     if col_repr:
 167         import matplotlib.colors as colors
 168         col = np.flip(np.divide([(230, 25, 75, 150), (245, 130, 48, 150), (255, 255, 25, 150),
 169                                 (210, 245, 60, 150), (60, 180, 75, 150), (70, 240, 240, 150),
 170                                 (0, 0, 128, 150), (145, 30, 180, 150), (240, 50, 230, 150),
 171                                 (128, 128, 128, 150)],255.),0)
 172         cm = colors.LinearSegmentedColormap.from_list('BuRd',col,len(contours))
 173         cols = ['%02x%02x%02x%02x' % tuple(255*np.flip(c)) for c in cm(range(cm.N))]
 174
 175         # Creating an array of dictionaries for each perimeter
 176         conts=[Dict({'text':time_iso[k],
 177             'LineStyle':{
 178                 'color': cols[k],
 179                 'width':'2.5',
 180             },
 181             'PolyStyle':{
 182                 'color':'66000086',
 183                 'colorMode':'random'
 184             },
 185             'time_begin':time_iso[k],
 186             'polygons': contours[k] }) for k in range(len(contours))]
 187     else:
 188         # Creating an array of dictionaries for each perimeter
 189         conts=[Dict({'text':time_iso[k],
 190             'LineStyle':{
 191                 'color':'ff081388',
 192                 'width':'2.5',
 193             },
 194             'PolyStyle':{
 195                 'color':'66000086',
 196                 'colorMode':'random'
 197             },
 198             'time_begin':time_iso[k],
 199             'polygons': contours[k] }) for k in range(len(contours))]
 200
 201     # Creating a dictionary to store the KML file information
 202     data=Dict({'name':'contours.kml',
 203         'folder_name':'Perimeters'})
 204     data.update({'contours': conts})
 205
 206     return data
 207
 208 if __name__ == "__main__":
 209     import os,sys
 210     result_file = 'result.mat'
 211     mgout_file = 'mgout.mat'
 212     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):
 213         print 'Loading the data...'
 214         result=loadmat(result_file)
 215         mgout=loadmat(mgout_file)
 216     else:
 217         print 'Error: file %s or %s not exist or not readable' % (result_file,mgout_file)
 218         sys.exit(1)
 219
 220     # Indexing the coordinates into the same as the multigrid solution
 221     xind=mgout['sm'][0]-1
 222     yind=mgout['sn'][0]-1
 223     x=np.array(result['fxlon'])
 224     xx=x[np.ix_(xind,yind)]
 225     y=np.array(result['fxlat'])
 226     yy=y[np.ix_(xind,yind)]
 227     tscale=mgout['tscale'][0]
 228     time_scale_num=mgout['time_scale_num'][0]
 229     zz=mgout['a']*tscale+time_scale_num[0]
 230
 231     print 'Computing the contours...'
 232     # Granules numeric times
 233     time_num_granules = result['time_num_granules'][0]
 234     data = get_contour_verts(xx, yy, zz, time_num_granules, contour_dt_hours=24, gauss_filter=False)
 235
 236     print 'Creating the KML file...'
 237     # Creating the KML file
 238     contour2kml(data,'perimeters.kml')
 239
 240     print 'perimeters.kml generated'
 241