forecast.py

   1 import matplotlib.pyplot as plt
   2 import numpy as np
   3 from datetime import timedelta
   4 from JPSSD import time_iso2num, time_iso2datetime, time_datetime2iso
   5 from utils import Dict
   6 import re, glob, sys, os
   7
   8
   9 def process_tign_g(lon,lat,tign_g,bounds,ctime,dx,dy,wrfout_file='',dt_for=3600.,plot=False):
  10     """
  11     Process forecast from lon, lat, and tign_g
  12
  13     :param lon: array of longitudes
  14     :param lat: array of latitudes
  15     :param tign_g: array of fire arrival time
  16     :param bounds: coordinate bounding box filtering to
  17     :param ctime: time char in wrfout of the last fire arrival time
  18     :param dx,dy: data resolution in meters
  19     :param dt_for: optional, temporal resolution to get the fire arrival time from in seconds
  20     :param wrfout_file: optional, to get the name of the file in the dictionary
  21
  22     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  23     Angel Farguell (angel.farguell@gmail.com) and James Haley, 2019-06-05
  24     """
  25
  26     if plot:
  27         fig = plt.figure()
  28
  29     # prefix of the elements in the dictionary
  30     prefix = "FOR"
  31     # confidences
  32     conf_fire = 50
  33     conf_nofire = 10
  34     # margin percentage
  35     margin = .5
  36     # scan and track dimensions of the observation (in km)
  37     scan = dx/1000.
  38     track = dy/1000.
  39     # initializing dictionary
  40     forecast = Dict({})
  41
  42     # ctime transformations
  43     ctime_iso = ctime.replace('_','T')
  44     ctime_datetime = time_iso2datetime(ctime_iso)
  45     # mask coordinates to bounding box
  46     mask = np.logical_and(np.logical_and(np.logical_and(lon >= bounds[0], lon <= bounds[1]), lat >= bounds[2]), lat <= bounds[3])
  47     # create a square subset of fire arrival time less than the maximum
  48     mtign = np.logical_and(mask, tign_g < tign_g.max())
  49     mlon = lon[mtign]
  50     mlat = lat[mtign]
  51     mlenlon = margin*(mlon.max()-mlon.min())
  52     mlenlat = margin*(mlat.max()-mlat.min())
  53     sbounds = (mlon.min()-mlenlon, mlon.max()+mlenlon, mlat.min()-mlenlat, mlat.max()+mlenlat)
  54     smask = np.logical_and(np.logical_and(np.logical_and(lon >= sbounds[0], lon <= sbounds[1]), lat >= sbounds[2]), lat <= sbounds[3])
  55
  56     # times to get fire arrival time from
  57     TT = np.arange(tign_g.min()-3*dt_for,tign_g.max(),dt_for)[0:-1]
  58     for T in TT:
  59         # fire arrival time to datetime
  60         time_datetime = ctime_datetime-timedelta(seconds=float(tign_g.max()-T)) # there is an error of about 10 seconds (wrfout not ending at exact time of simulation)
  61         # create ISO time of the date
  62         time_iso = time_datetime2iso(time_datetime)
  63         # create numerical time from the ISO time
  64         time_num = time_iso2num(time_iso)
  65         # create time stamp
  66         time_data = '_A%04d%03d_%02d%02d_%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
  67                                                 time_datetime.hour, time_datetime.minute, time_datetime.second)
  68         # create acquisition date
  69         acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
  70         # create acquisition time
  71         acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
  72
  73         print 'Retrieving forecast from %s' % time_iso
  74
  75         # masks
  76         fire = tign_g <= T
  77         nofire = np.logical_and(tign_g > T, np.logical_or(tign_g != tign_g.max(), smask))
  78         # coordinates masked
  79         lon_fire = lon[np.logical_and(mask,fire)]
  80         lat_fire = lat[np.logical_and(mask,fire)]
  81         lon_nofire = lon[np.logical_and(mask,nofire)]
  82         lat_nofire = lat[np.logical_and(mask,nofire)]
  83         # create general arrays
  84         lons = np.concatenate((lon_nofire,lon_fire))
  85         lats = np.concatenate((lat_nofire,lat_fire))
  86         fires = np.concatenate((5*np.ones(lon_nofire.shape),9*np.ones(lon_fire.shape)))
  87
  88         # plot results
  89         if plot:
  90             plt.ion()
  91             plt.cla()
  92             plt.plot(lons[fires==5],lats[fires==5],'g.')
  93             plt.plot(lons[fires==9],lats[fires==9],'r.')
  94             plt.pause(.0001)
  95             plt.show()
  96
  97         # update perimeters dictionary
  98         forecast.update({prefix + time_data: Dict({'file': wrfout_file, 'time_tign_g': T, 'lon': lons, 'lat': lats,
  99                             'fire': fires, 'conf_fire': np.array(conf_fire*np.ones(lons[fires==9].shape)),
 100                             'lon_fire': lons[fires==9], 'lat_fire': lats[fires==9], 'lon_nofire': lons[fires==5], 'lat_nofire': lats[fires==5],
 101                             'scan_fire': scan*np.ones(lons[fires==9].shape), 'track_fire': track*np.ones(lons[fires==9].shape),
 102                             'conf_nofire': np.array(conf_nofire*np.ones(lons[fires==5].shape)),
 103                             'scan_nofire': scan*np.ones(lons[fires==5].shape), 'track_nofire': track*np.ones(lons[fires==9].shape),
 104                             'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 105
 106     return forecast
 107
 108 def process_forecast_wrfout(wrfout_file,bounds,plot=False):
 109     """
 110     Process forecast from a wrfout.
 111
 112     :param dst: path to wrfout file
 113     :param bounds: coordinate bounding box filtering to
 114
 115     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 116     Angel Farguell (angel.farguell@gmail.com), 2019-06-05
 117     """
 118
 119     # read file
 120     import netCDF4 as nc
 121     try:
 122         data = nc.Dataset(wrfout_file)
 123     except Exception as e:
 124         print 'Warning: No netcdf file %s in the path' % wrfout_file
 125         return []
 126     # current time
 127     ctime = ''.join(data['Times'][-1])
 128     # getting rid of strip
 129     atmlenx = len(data.dimensions['west_east'])
 130     atmleny = len(data.dimensions['south_north'])
 131     staglenx = len(data.dimensions['west_east_stag'])
 132     stagleny = len(data.dimensions['south_north_stag'])
 133     dimlenx = len(data.dimensions['west_east_subgrid'])
 134     dimleny = len(data.dimensions['south_north_subgrid'])
 135     ratiox = dimlenx/max(staglenx,atmlenx+1)
 136     ratioy = dimleny/max(stagleny,atmleny+1)
 137     lenx = dimlenx-ratiox
 138     leny = dimleny-ratioy
 139     # coordinates
 140     lon = data['FXLONG'][0][0:lenx,0:leny]
 141     lat = data['FXLAT'][0][0:lenx,0:leny]
 142     # fire arrival time
 143     tign_g = data['TIGN_G'][0][0:lenx,0:leny]
 144     # resolutions
 145     DX = data.getncattr('DX')
 146     DY = data.getncattr('DY')
 147     sx = data.dimensions['west_east_subgrid'].size/data.dimensions['west_east_stag'].size
 148     sy = data.dimensions['south_north_subgrid'].size/data.dimensions['south_north_stag'].size
 149     dx = DX/sx
 150     dy = DY/sy
 151     # create forecast
 152     forecast = process_tign_g(lon,lat,tign_g,bounds,ctime,dx,dy,wrfout_file=wrfout_file,plot=plot)
 153     # close netcdf file
 154     data.close()
 155
 156     return forecast
 157
 158
 159 if __name__ == "__main__":
 160     import saveload as sl
 161     real = False
 162
 163     if real:
 164         plot = True
 165         bounds = (-113.85068, -111.89413, 39.677563, 41.156837)
 166         dst = './patch/wrfout_patch'
 167         f = process_forecast_wrfout(dst,bounds,plot=plot)
 168         sl.save(f,'forecast')
 169     else:
 170         from infrared_perimeters import process_ignitions
 171         from setup import process_detections
 172         dst = 'ideal_test'
 173         plot = False
 174         ideal = sl.load(dst)
 175         kk = 4
 176         data = process_tign_g(ideal['lon'][::kk,::kk],ideal['lat'][::kk,::kk],ideal['tign_g'][::kk,::kk],ideal['bounds'],ideal['ctime'],ideal['dx'],ideal['dy'],wrfout_file='ideal',dt_for=ideal['dt'],plot=plot)
 177         if 'point' in ideal.keys():
 178             p = [[ideal['point'][0]],[ideal['point'][1]],[ideal['point'][2]]]
 179             data.update(process_ignitions(p,ideal['bounds']))
 180         elif 'points' in ideal.keys():
 181             p = [[point[0] for point in ideal['points']],[point[1] for point in ideal['points']],[point[2] for point in ideal['points']]]
 182             data.update(process_ignitions(p,ideal['bounds']))
 183         etime = time_iso2num(ideal['ctime'].replace('_','T'))
 184         time_num_int = (etime-ideal['tign_g'].max(),etime)
 185         sl.save((data,ideal['lon'],ideal['lat'],time_num_int),'data')
 186         process_detections(data,ideal['lon'],ideal['lat'],time_num_int)
 187