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