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 and infrared_perimeters.py file:
  15 #               *) Find granules overlaping fire domain and time interval.
  16 #               *) Download Active Satellite Data.
  17 #               *) Read and process Active Satellite Data files.
  18 #               *) Process ignitions.
  19 #               *) Read and process infrared perimeters files.
  20 #               *) Save observed data information in 'data' file.
  21 #       2) Methods from interpolation.py, JPSSD.py, and plot_pixels.py files:
  22 #               *) Write KML file 'fire_detections.kml' with fire detection pixels (satellite, ignitions and perimeters).
  23 #               *) Write KMZ file 'googlearth.kmz' with saved images and KML file of observed data.
  24 #       3) Method process_detections from setup.py file:
  25 #               *) Sort all the granules from all the sources in time order.
  26 #               *) Construct upper and lower bounds using a mask to prevent ground after fire.
  27 #               *) Save results in 'results.mat' file.
  28 #       4) Methods preprocess_data_svm and SVM3 from svm.py file:
  29 #               *) Preprocess bounds as an input of Support Vector Machine method.
  30 #               *) Run Support Vector Machine method.
  31 #               *) Save results in svm.mat file.
  32 #       5) Methods from contline.py and contour2kml.py files:
  33 #               *) Construct a smooth contour line representation of the fire arrival time.
  34 #               *) Write the contour lines in a KML file called 'perimeters_svm.kml'.
  35 #
  36 # OUTPUTS
  37 #       - 'data': binary file containing satellite granules information.
  38 #       - 'result.mat': matlab file containing upper and lower bounds (U and L) from satellite data.
  39 #       - 'svm.mat': matlab file containing the solution to the Support Vector Machine execution.
  40 #                                Contains estimation of the fire arrival time in tign_g variable.
  41 #       - 'fire_detections.kml': KML file with fire detection pixels (satellite, ignitions and perimeters).
  42 #       - 'googlearth.kmz': KMZ file with saved images and KML file of observed data.
  43 #       - 'perimeters_svm.kml': KML file with perimeters from estimation of the fire arrival time using SVM.
  44 #
  45 # Developed in Python 2.7.15 :: Anaconda, Inc.
  46 # Angel Farguell (angel.farguell@gmail.com), 2019-04-29
  47 #---------------------------------------------------------------------------------------------------------------------
  48 from JPSSD import read_fire_mesh, retrieve_af_data, sdata2json, json2kml, time_iso2num
  49 from interpolation import sort_dates
  50 from setup import process_detections
  51 from infrared_perimeters import process_ignitions, process_infrared_perimeters
  52 from svm import preprocess_data_svm, SVM3
  53 from mpl_toolkits.basemap import Basemap
  54 from plot_pixels import basemap_scatter_mercator, create_kml
  55 from contline import get_contour_verts
  56 from contour2kml import contour2kml
  57 import saveload as sl
  58 from utils import Dict
  59 from scipy.io import loadmat, savemat
  60 import numpy as np
  61 import datetime as dt
  62 import sys, os, re
  63 from time import time
  64
  65 # plot observed information
  66 plot_observed = False
  67
  68 # if ignitions are known: ([lons],[lats],[dates]) where lons and lats in degrees and dates in ESMF format
  69 # 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'])
  70 igns = None
  71 # if infrared perimeters: path to KML files
  72 # examples: perim_path = './pioneer/perim'
  73 perim_path = ''
  74
  75 satellite_file = 'data'
  76 fire_file = 'fire_detections.kml'
  77 gearth_file = 'googlearth.kmz'
  78 bounds_file = 'result.mat'
  79 svm_file = 'svm.mat'
  80 contour_file = 'perimeters_svm.kml'
  81
  82 def exist(path):
  83         return (os.path.isfile(path) and os.access(path,os.R_OK))
  84
  85 satellite_exists = exist(satellite_file)
  86 fire_exists = exist(fire_file)
  87 gearth_exists = exist(gearth_file)
  88 bounds_exists = exist(bounds_file)
  89
  90 if len(sys.argv) != 4 and (not bounds_exists) and (not satellite_exists):
  91         print 'Error: python %s wrfout start_time days' % sys.argv[0]
  92         print ' * wrfout - string, wrfout file of WRF-SFIRE simulation'
  93         print ' * start_time - string, YYYYMMDDHHMMSS where: '
  94         print '         YYYY - year'
  95         print '         MM - month'
  96         print '         DD - day'
  97         print '         HH - hour'
  98         print '         MM - minute'
  99         print '         SS - second'
 100         print ' * days - float, number of days of simulation (can be less than a day)'
 101         print 'or link an existent file %s or %s' % (satellite_file,bounds_file)
 102         sys.exit(0)
 103
 104 t_init = time()
 105
 106 print ''
 107 if bounds_exists:
 108         print '>> File %s already created! Skipping all satellite processing <<' % bounds_file
 109         print 'Loading from %s...' % bounds_file
 110         result = loadmat(bounds_file)
 111         # Taking necessary variables from result dictionary
 112         scale = result['time_scale_num'][0]
 113         time_num_granules = result['time_num_granules'][0]
 114         time_num_interval = result['time_num'][0]
 115         lon = np.array(result['fxlon']).astype(float)
 116         lat = np.array(result['fxlat']).astype(float)
 117 else:
 118         if satellite_exists:
 119                 print '>> File %s already created! Skipping satellite retrieval <<' % satellite_file
 120                 print 'Loading from %s...' % satellite_file
 121                 data,fxlon,fxlat,time_num = sl.load(satellite_file)
 122                 bbox = [fxlon.min(),fxlon.max(),fxlat.min(),fxlat.max()]
 123         else:
 124                 print '>> Reading the fire mesh <<'
 125                 sys.stdout.flush()
 126                 fxlon,fxlat,bbox,time_esmf = read_fire_mesh(sys.argv[1])
 127                 # converting times to ISO
 128                 dti = dt.datetime.strptime(sys.argv[2],'%Y%m%d%H%M%S')
 129                 time_start_iso = '%d-%02d-%02dT%02d:%02d:%02dZ' % (dti.year,dti.month,dti.day,dti.hour,dti.minute,dti.second)
 130                 dtf = dti+dt.timedelta(days=float(sys.argv[3]))
 131                 time_final_iso = '%d-%02d-%02dT%02d:%02d:%02dZ' % (dtf.year,dtf.month,dtf.day,dtf.hour,dtf.minute,dtf.second)
 132                 time_iso = (time_start_iso,time_final_iso)
 133
 134                 print ''
 135                 print '>> Retrieving satellite data <<'
 136                 sys.stdout.flush()
 137                 data = retrieve_af_data(bbox,time_iso)
 138                 if igns:
 139                         data.update(process_ignitions(igns,bounds=bbox))
 140                 if perim_path:
 141                         data.update(process_infrared_perimeters(perim_path,bounds=bbox))
 142
 143                 if data:
 144                         print ''
 145                         print '>> Saving satellite data file (data) <<'
 146                         sys.stdout.flush()
 147                         time_num = map(time_iso2num,time_iso)
 148                         sl.save((data,fxlon,fxlat,time_num),satellite_file)
 149                         print 'data file saved correctly!'
 150                 else:
 151                         print ''
 152                         print 'ERROR: No data obtained...'
 153                         sys.exit(1)
 154
 155         print ''
 156         if (not fire_exists) or (not gearth_exists and plot_observed):
 157                 print '>> Generating KML of fire and ground detections <<'
 158                 sys.stdout.flush()
 159                 # sort the granules by dates
 160                 sdata=sort_dates(data)
 161         if fire_exists:
 162                 print '>> File %s already created! <<' % fire_file
 163         else:
 164                 # writting fire detections file
 165                 print 'writting KML with fire detections'
 166                 keys = ['latitude','longitude','brightness','scan','track','acq_date','acq_time','satellite','instrument','confidence','bright_t31','frp','scan_angle']
 167                 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']
 168                 prods = {'AF':'Active Fires','FRP':'Fire Radiative Power','TF':'Temporal Fire coloring'}
 169                 # filter out perimeter, ignition, and forecast information (too many pixels)
 170                 regex = re.compile(r'^((?!(PER_A|IGN_A|FOR_A)).)*$')
 171                 nsdata = [d for d in sdata if regex.match(d[0])]
 172                 # compute number of elements for each granule
 173                 N = [len(d[1]['lat_fire']) if 'lat_fire' in d[1] else 0 for d in nsdata]
 174                 # transform dictionary notation to json notation
 175                 json = sdata2json(nsdata,keys,dkeys,N)
 176                 # write KML file from json notation
 177                 json2kml(json,fire_file,bbox,prods)
 178         if gearth_exists or not plot_observed:
 179                 print ''
 180                 print '>> File %s already created! <<' % gearth_file
 181         else:
 182                 # creating KMZ overlay of each information
 183                 # create the Basemap to plot into
 184                 bmap = Basemap(projection='merc',llcrnrlat=bbox[2], urcrnrlat=bbox[3], llcrnrlon=bbox[0], urcrnrlon=bbox[1])
 185                 # initialize array
 186                 kmld = []
 187                 # for each observed information
 188                 for idx, g in enumerate(sdata):
 189                         # plot a scatter basemap
 190                         raster_png_data,corner_coords = basemap_scatter_mercator(g[1],bbox,bmap)
 191                         # compute bounds
 192                         bounds = (corner_coords[0][0],corner_coords[1][0],corner_coords[0][1],corner_coords[2][1])
 193                         # create png name
 194                         pngfile = g[0]+'.png'
 195                         # create timestamp for KML
 196                         timestamp = g[1].acq_date + 'T' + g[1].acq_time[0:2] + ':' + g[1].acq_time[2:4] + 'Z'
 197                         # write PNG file
 198                         with open(pngfile, 'w') as f:
 199                                 f.write(raster_png_data)
 200                         print '> File %s saved.' % g[0]
 201                         # append dictionary information for the KML creation
 202                         kmld.append(Dict({'name': g[0], 'png_file': pngfile, 'bounds': bbox, 'time': timestamp}))
 203                 # create KML
 204                 create_kml(kmld,'./doc.kml')
 205                 # create KMZ with all the PNGs included
 206                 os.system('zip -r %s doc.kml *_A*_*.png' % gearth_file)
 207                 print 'Created file %s' % gearth_file
 208                 # eliminate images and KML after creation of KMZ
 209                 os.system('rm doc.kml *_A*_*.png')
 210
 211         print ''
 212         print '>> Processing satellite data <<'
 213         sys.stdout.flush()
 214         result = process_detections(data,fxlon,fxlat,time_num)
 215         # Taking necessary variables from result dictionary
 216         scale = result['time_scale_num']
 217         time_num_granules = result['time_num_granules']
 218         time_num_interval = result['time_num']
 219         lon = np.array(result['fxlon']).astype(float)
 220         lat = np.array(result['fxlat']).astype(float)
 221
 222 U = np.array(result['U']).astype(float)
 223 L = np.array(result['L']).astype(float)
 224 T = np.array(result['T']).astype(float)
 225
 226 if 'C' in result.keys():
 227         conf = np.array(result['C'])
 228 else:
 229         conf = None
 230
 231 print ''
 232 print '>> Preprocessing the data <<'
 233 sys.stdout.flush()
 234 X,y,c = preprocess_data_svm(lon,lat,U,L,T,scale,time_num_granules,C=conf)
 235
 236 print ''
 237 print '>> Running Support Vector Machine <<'
 238 sys.stdout.flush()
 239 C = 500.
 240 kgam = 10.
 241 F = SVM3(X,y,C=C,kgam=kgam,fire_grid=(lon,lat))
 242
 243 print ''
 244 print '>> Saving the results <<'
 245 sys.stdout.flush()
 246 tscale = 24*3600 # scale from seconds to days
 247 # Fire arrival time in seconds from the begining of the simulation
 248 tign_g = np.array(F[2])*float(tscale)+scale[0]-time_num_interval[0]
 249 # Creating the dictionary with the results
 250 svm = {'dxlon': lon, 'dxlat': lat, 'U': U/tscale, 'L': L/tscale,
 251                 'fxlon': F[0], 'fxlat': F[1], 'Z': F[2],
 252                 'tign_g': tign_g,
 253                 'tscale': tscale, 'time_num_granules': time_num_granules,
 254                 'time_scale_num': scale, 'time_num': time_num_interval}
 255 # Save resulting file
 256 savemat(svm_file, mdict=svm)
 257 print 'The results are saved in svm.mat file'
 258
 259 print ''
 260 print '>> Computing contour lines of the fire arrival time <<'
 261 print 'Computing the contours...'
 262 try:
 263         # Granules numeric times
 264         Z = F[2]*tscale+scale[0]
 265         # Creating contour lines
 266         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)
 267         print 'Creating the KML file...'
 268         # Creating the KML file
 269         contour2kml(contour_data,contour_file)
 270         print 'The resulting contour lines are saved in perimeters_svm.kml file'
 271 except:
 272         print 'Warning: contour creation problem'
 273         print 'Run: python contlinesvm.py'
 274
 275 print ''
 276 print '>> DONE <<'
 277 t_final = time()
 278 print 'Elapsed time for all the process: %ss.' % str(abs(t_final-t_init))
 279 sys.exit()