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