process.py

   1 # process.py
   2 #
   3 # DESCRIPTION
   4 # Driver python code to estimate fire arrival time using Active Fire Satellite Data
   5 #
   6 # INPUTS
   7 # In the existence of a 'data' satellite granules file and/or 'results.mat' bounds file, any input is necessary.
   8 # Otherwise:
   9 #       wrfout - path to a simulation wrfout file (containing FXLON and FXLAT coordinates).
  10 #       start_time - date string with format: YYYYMMDDHHMMSS
  11 #       days - length of the simulation in decimal days
  12 #
  13 # OVERFLOW
  14 #       1) Methods from JPSSD.py file:
  15 #               *) Find granules overlaping fire domain and time interval.
  16 #               *) Download Active Satellite Data.
  17 #               *) Read Active Satellite Data files.
  18 #               *) Save satellite granule information in 'data' file.
  19 #       2) Methods from interpolation.py and JPSSD.py files:
  20 #               *) Write KML file 'fire_detections.kml' with fire detection pixels.
  21 #               *) Write KML file 'nofire.kml' with saved ground detection pixels.
  22 #       3) Method process_detections from setup.py file:
  23 #               *) Sort all the granules from all the sources in time order.
  24 #               *) Construct upper and lower bounds using a mask to prevent ground after fire.
  25 #               *) Save results in 'results.mat' file.
  26 #       4) Methods preprocess_data_svm and SVM3 from svm.py file:
  27 #               *) Preprocess bounds as an input of Support Vector Machine method.
  28 #               *) Run Support Vector Machine method.
  29 #               *) Save results in svm.mat file.
  30 #       5) Methods from contline.py and contour2kml.py files:
  31 #               *) Construct a smooth contour line representation of the fire arrival time.
  32 #               *) Write the contour lines in a KML file called 'perimeters_svm.kml'.
  33 #
  34 # OUTPUTS
  35 #       - 'data': binary file containing satellite granules information.
  36 #       - 'result.mat': matlab file containing upper and lower bounds (U and L) from satellite data.
  37 #       - 'svm.mat': matlab file containing the solution to the Support Vector Machine execution.
  38 #                                Contains estimation of the fire arrival time in tign_g variable.
  39 #       - 'fire_detections.kml': KML file with fire satellite detection pixels.
  40 #       - 'nofire.kml': KML file with saved ground satellite detection pixels.
  41 #       - 'perimeters_svm.kml': KML file with perimeters from estimation of the fire arrival time using SVM.
  42 #
  43 # Developed in Python 2.7.15 :: Anaconda, Inc.
  44 # Angel Farguell (angel.farguell@gmail.com), 2019-04-29
  45 #---------------------------------------------------------------------------------------------------------------------
  46 from JPSSD import read_fire_mesh, retrieve_af_data, sdata2json, json2kml, time_iso2num
  47 from interpolation import sort_dates
  48 from setup import process_detections
  49 from infrared_perimeters import process_ignitions, process_infrared_perimeters
  50 from svm import preprocess_data_svm, SVM3
  51 from contline import get_contour_verts
  52 from contour2kml import contour2kml
  53 import saveload as sl
  54 from scipy.io import loadmat, savemat
  55 import numpy as np
  56 import datetime as dt
  57 import sys
  58 import os
  59 from time import time
  60
  61 # if ignitions are known: ([lons],[lats],[dates]) where lons and lats in degrees and dates in ESMF format
  62 # examples: igns = ([100],[45],['2015-05-15T20:09:00']) or igns = ([100,105],[45,39],['2015-05-15T20:09:00','2015-05-15T23:09:00'])
  63 igns = None
  64 # if infrared perimeters: path to KML files
  65 # examples: perim_path = './pioneer/perim'
  66 perim_path = ''
  67
  68 satellite_file = 'data'
  69 fire_file = 'fire_detections.kml'
  70 ground_file = 'nofire.kml'
  71 bounds_file = 'result.mat'
  72 svm_file = 'svm.mat'
  73 contour_file = 'perimeters_svm.kml'
  74
  75 def exist(path):
  76         return (os.path.isfile(path) and os.access(path,os.R_OK))
  77
  78 satellite_exists = exist(satellite_file)
  79 fire_exists = exist(fire_file)
  80 ground_exists = exist(ground_file)
  81 bounds_exists = exist(bounds_file)
  82
  83 if len(sys.argv) != 4 and (not bounds_exists) and (not satellite_exists):
  84         print 'Error: python %s wrfout start_time days' % sys.argv[0]
  85         print ' * wrfout - string, wrfout file of WRF-SFIRE simulation'
  86         print ' * start_time - string, YYYYMMDDHHMMSS where: '
  87         print '         YYYY - year'
  88         print '         MM - month'
  89         print '         DD - day'
  90         print '         HH - hour'
  91         print '         MM - minute'
  92         print '         SS - second'
  93         print ' * days - float, number of days of simulation (can be less than a day)'
  94         print 'or link an existent file %s or %s' % (satellite_file,bounds_file)
  95         sys.exit(0)
  96
  97 t_init = time()
  98
  99 print ''
 100 if bounds_exists:
 101         print '>> File %s already created! Skipping all satellite processing <<' % bounds_file
 102         print 'Loading from %s...' % bounds_file
 103         result = loadmat(bounds_file)
 104         # Taking necessary variables from result dictionary
 105         scale = result['time_scale_num'][0]
 106         time_num_granules = result['time_num_granules'][0]
 107         time_num_interval = result['time_num'][0]
 108         lon = np.array(result['fxlon']).astype(float)
 109         lat = np.array(result['fxlat']).astype(float)
 110 else:
 111         if satellite_exists:
 112                 print '>> File %s already created! Skipping satellite retrieval <<' % satellite_file
 113                 print 'Loading from %s...' % satellite_file
 114                 data,fxlon,fxlat,time_num = sl.load(satellite_file)
 115                 bbox = [fxlon.min(),fxlon.max(),fxlat.min(),fxlat.max()]
 116         else:
 117                 print '>> Reading the fire mesh <<'
 118                 sys.stdout.flush()
 119                 fxlon,fxlat,bbox,time_esmf = read_fire_mesh(sys.argv[1])
 120                 # converting times to ISO
 121                 dti = dt.datetime.strptime(sys.argv[2],'%Y%m%d%H%M%S')
 122                 time_start_iso = '%d-%02d-%02dT%02d:%02d:%02dZ' % (dti.year,dti.month,dti.day,dti.hour,dti.minute,dti.second)
 123                 dtf = dti+dt.timedelta(days=float(sys.argv[3]))
 124                 time_final_iso = '%d-%02d-%02dT%02d:%02d:%02dZ' % (dtf.year,dtf.month,dtf.day,dtf.hour,dtf.minute,dtf.second)
 125                 time_iso = (time_start_iso,time_final_iso)
 126
 127                 print ''
 128                 print '>> Retrieving satellite data <<'
 129                 sys.stdout.flush()
 130                 data = retrieve_af_data(bbox,time_iso)
 131                 if igns:
 132                         data.update(process_ignitions(igns,bounds=bbox))
 133                 if perim_path:
 134                         data.update(process_infrared_perimeters(perim_path,bounds=bbox))
 135
 136                 if data:
 137                         print ''
 138                         print '>> Saving satellite data file (data) <<'
 139                         sys.stdout.flush()
 140                         time_num = map(time_iso2num,time_iso)
 141                         sl.save((data,fxlon,fxlat,time_num),satellite_file)
 142                         print 'data file saved correctly!'
 143                 else:
 144                         print ''
 145                         print 'ERROR: No data obtained...'
 146                         sys.exit(1)
 147
 148         print ''
 149         if (not fire_exists) or (not ground_exists):
 150                 print '>> Generating KML of fire and ground detections <<'
 151                 sys.stdout.flush()
 152                 # sort the granules by dates
 153                 sdata=sort_dates(data)
 154         if fire_exists:
 155                 print '>> File %s already created! <<' % fire_file
 156         else:
 157                 # writting fire detections file
 158                 print 'writting KML with fire detections'
 159                 keys = ['latitude','longitude','brightness','scan','track','acq_date','acq_time','satellite','instrument','confidence','bright_t31','frp','scan_angle']
 160                 dkeys = ['lat_fire','lon_fire','brig_fire','scan_fire','track_fire','acq_date','acq_time','sat_fire','instrument','conf_fire','t31_fire','frp_fire','scan_angle_fire']
 161                 prods = {'AF':'Active Fires','FRP':'Fire Radiative Power'}
 162                 N = [len(d[1]['lat_fire']) if 'lat_fire' in d[1] else 0 for d in sdata]
 163                 json = sdata2json(sdata,keys,dkeys,N)
 164                 json2kml(json,fire_file,bbox,prods)
 165         if ground_exists:
 166                 print ''
 167                 print '>> File %s already created! <<' % ground_file
 168         else:
 169                 # writting ground detections file
 170                 print 'writting KML with ground'
 171                 keys = ['latitude','longitude','scan','track','acq_date','acq_time','satellite','instrument','scan_angle']
 172                 dkeys = ['lat_nofire','lon_nofire','scan_nofire','track_nofire','acq_date','acq_time','sat_fire','instrument','scan_angle_nofire']
 173                 prods = {'NF':'No Fire'}
 174                 N = [len(d[1]['lat_nofire']) if 'lat_nofire' in d[1] else 0 for d in sdata]
 175                 json = sdata2json(sdata,keys,dkeys,N)
 176                 json2kml(json,ground_file,bbox,prods)
 177
 178         print ''
 179         print '>> Processing satellite data <<'
 180         sys.stdout.flush()
 181         result = process_detections(data,fxlon,fxlat,time_num)
 182         # Taking necessary variables from result dictionary
 183         scale = result['time_scale_num']
 184         time_num_granules = result['time_num_granules']
 185         time_num_interval = result['time_num']
 186         lon = np.array(result['fxlon']).astype(float)
 187         lat = np.array(result['fxlat']).astype(float)
 188
 189 U = np.array(result['U']).astype(float)
 190 L = np.array(result['L']).astype(float)
 191 T = np.array(result['T']).astype(float)
 192
 193 if 'C' in result.keys():
 194         conf = np.array(result['C'])
 195 else:
 196         conf = None
 197
 198 print ''
 199 print '>> Preprocessing the data <<'
 200 sys.stdout.flush()
 201 X,y,c = preprocess_data_svm(lon,lat,U,L,T,scale,time_num_granules,C=conf)
 202
 203 print ''
 204 print '>> Running Support Vector Machine <<'
 205 sys.stdout.flush()
 206 C = 10.
 207 kgam = 10.
 208 F = SVM3(X,y,C=C,kgam=kgam,fire_grid=(lon,lat))
 209
 210 print ''
 211 print '>> Saving the results <<'
 212 sys.stdout.flush()
 213 tscale = 24*3600 # scale from seconds to days
 214 # Fire arrival time in seconds from the begining of the simulation
 215 tign_g = F[2]*float(tscale)+scale[0]-time_num_interval[0]
 216 # Creating the dictionary with the results
 217 svm = {'dxlon': lon, 'dxlat': lat, 'U': U/tscale, 'L': L/tscale,
 218         'fxlon': F[0], 'fxlat': F[1], 'Z': F[2],
 219         'tign_g': tign_g,
 220         'tscale': tscale, 'time_num_granules': time_num_granules,
 221         'time_scale_num': scale, 'time_num': time_num_interval}
 222 # Save resulting file
 223 savemat(svm_file, mdict=svm)
 224 print 'The results are saved in svm.mat file'
 225
 226 print ''
 227 print '>> Computing contour lines of the fire arrival time <<'
 228 print 'Computing the contours...'
 229 # Fire arrival time in seconds from an old date
 230 Z = F[2]*tscale+scale[0]
 231 # Granules numeric times
 232 contour_data = get_contour_verts(F[0], F[1], Z, time_num_granules, contour_dt_hours=6, contour_dt_init=6, contour_dt_final=6)
 233 print 'Creating the KML file...'
 234 # Creating the KML file
 235 contour2kml(contour_data,contour_file)
 236 print 'The resulting contour lines are saved in perimeters_svm.kml file'
 237
 238 print ''
 239 print '>> DONE <<'
 240 t_final = time()
 241 print 'Elapsed time for all the process: %ss.' % str(abs(t_final-t_init))
 242 sys.exit()