synthetic.py

   1 import numpy as np
   2 import matplotlib.pyplot as plt
   3 from mpl_toolkits.mplot3d import axes3d
   4 from forecast import process_tign_g
   5 from JPSSD import time_iso2num
   6 from setup import process_detections
   7 from infrared_perimeters import process_ignitions
   8 import saveload as sl
   9
  10 # Realistic bounds
  11 nx = ny = 201
  12 bounds = (-113.85068, -111.89413, 39.677563, 41.156837)
  13
  14 # Creation of synthetic coordinates
  15 lon,lat = np.meshgrid(np.linspace(bounds[0],bounds[1],nx),
  16                     np.linspace(bounds[2],bounds[3],ny))
  17 ilon = .5*(bounds[0]+bounds[1])
  18 ilat = .5*(bounds[2]+bounds[3])
  19
  20 # Scaling fire arrival time to be from 0 to 1000
  21 lenlon = abs(bounds[1]-bounds[0])*.5
  22 lenlat = abs(bounds[3]-bounds[2])*.5
  23 clon = 5e5/(lenlon**2)
  24 clat = 5e5/(lenlat**2)
  25
  26 # Creation of cone synthetic fire arrival time
  27 tign_g = np.minimum(600,10+np.sqrt(clon*(lon-ilon)**2+clat*(lat-ilat)**2))
  28 print 'min tign_g %gs, max tign_g %gs' % (tign_g.min(),tign_g.max())
  29 fig = plt.figure()
  30 ax = fig.gca(projection='3d')
  31 ax.contour(lon,lat,tign_g,20)
  32 ax.scatter(-113.25, ilat, 150, color='r')
  33 plt.show()
  34
  35 # Some more set-ups needed
  36 dx = dy = 1000.
  37 dt = 10
  38 ctime = '2013-08-14_00:00:00'
  39
  40 # Forecast creation
  41 data = process_tign_g(lon,lat,tign_g,bounds,ctime,dx,dy,wrfout_file='ideal',dt_for=dt,plot=True)
  42 #points = [[-113.25], [ilat], ['2013-08-13T23:52:30']]
  43 points = [[-113.25, -113.25, -113.25], [ilat], ['2013-08-13T23:52:30']]
  44 data.update(process_ignitions(points, bounds))
  45
  46 # Process detections
  47 isotime = ctime.replace('_','T')
  48 ftime = time_iso2num(isotime)
  49 time_num = (ftime-tign_g.max(),ftime)
  50 sl.save((data,lon,lat,time_num),'data')
  51 result = process_detections(data,lon,lat,time_num)