forecast.py

   1 import matplotlib.pyplot as plt
   2 import numpy as np
   3 import os.path as osp
   4 import netCDF4 as nc
   5 from datetime import timedelta
   6 from JPSSD import time_iso2num, time_iso2datetime, time_datetime2iso
   7 from utils import Dict
   8 import saveload as sl
   9 import re, glob, sys, os
  10
  11
  12 def process_forecast(wrfout_file,bounds,plot=False):
  13     """
  14     Process infrared perimeters the same way than satellite data.
  15
  16     :param dst: path to kml perimeter files
  17     :param bounds: coordinate bounding box filtering to
  18
  19     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  20     Angel Farguell (angel.farguell@gmail.com), 2019-06-05
  21     """
  22
  23     # prefix of the elements in the dictionary
  24     prefix = "FOR"
  25     # confidences
  26     conf_fire = 70
  27     conf_nofire = 70
  28     # margin percentage
  29     margin = .1
  30     # initializing dictionary
  31     forecast = Dict({})
  32
  33     # read file
  34     try:
  35         data = nc.Dataset(wrfout_file)
  36     except Exception as e:
  37         print 'Warning: No netcdf file %s in the path' % wrfout_file
  38         return []
  39     # current time
  40     ctime = ''.join(data['Times'][-1])
  41     ctime_iso = ctime.replace('_','T')
  42     ctime_datetime = time_iso2datetime(ctime_iso)
  43     # getting rid of strip
  44     atmlenx = len(data.dimensions['west_east'])
  45     atmleny = len(data.dimensions['south_north'])
  46     staglenx = len(data.dimensions['west_east_stag'])
  47     stagleny = len(data.dimensions['south_north_stag'])
  48     dimlenx = len(data.dimensions['west_east_subgrid'])
  49     dimleny = len(data.dimensions['south_north_subgrid'])
  50     ratiox = dimlenx/max(staglenx,atmlenx+1)
  51     ratioy = dimleny/max(stagleny,atmleny+1)
  52     lenx = dimlenx-ratiox
  53     leny = dimleny-ratioy
  54     # coordinates
  55     lon = data['FXLONG'][0][0:lenx,0:leny]
  56     lat = data['FXLAT'][0][0:lenx,0:leny]
  57     # mask coordinates to bounding box
  58     mask = np.logical_and(np.logical_and(np.logical_and(lon>bounds[0],lon<bounds[1]),lat>bounds[2]),lat<bounds[3])
  59     # fire arrival time
  60     tign_g = data['TIGN_G'][0][0:lenx,0:leny]
  61     # create a square subset of fire arrival time less than the maximum
  62     mtign = np.logical_and(mask,tign_g < tign_g.max())
  63     mlon = lon[mtign]
  64     mlat = lat[mtign]
  65     mlen = margin*(mlon.max()-mlon.min())
  66     sbounds = (mlon.min()-mlen, mlon.max()+mlen, mlat.min()-mlen, mlat.max()+mlen)
  67     smask = np.logical_and(np.logical_and(np.logical_and(lon>sbounds[0],lon<sbounds[1]),lat>sbounds[2]),lat<sbounds[3])
  68     # resolutions
  69     dx = data.getncattr('DX')
  70     dy = data.getncattr('DY')
  71     # scan and track dimensions of the observation (in km)
  72     scan = dx/1000.
  73     track = dy/1000.
  74     # close netcdf file
  75     data.close()
  76
  77     # times to get fire arrival time from
  78     dt_forecast = 1800 # in seconds
  79     TT = np.arange(tign_g.min(),tign_g.max(),dt_forecast)[0:-1]
  80     for T in TT:
  81         # fire arrival time to datetime
  82         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)
  83         # create ISO time of the date
  84         time_iso = time_datetime2iso(time_datetime)
  85         # create numerical time from the ISO time
  86         time_num = time_iso2num(time_iso)
  87         # create time stamp
  88         time_data = '_A%04d%03d_%02d%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
  89                                                 time_datetime.hour, time_datetime.minute)
  90         # create acquisition date
  91         acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
  92         # create acquisition time
  93         acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
  94
  95         print 'Retrieving forecast from %s' % time_iso
  96
  97         # masks
  98         fire = tign_g <= T
  99         nofire = np.logical_and(tign_g > T, np.logical_or(tign_g != tign_g.max(), smask))
 100         # coordinates masked
 101         lon_fire = lon[np.logical_and(mask,fire)]
 102         lat_fire = lat[np.logical_and(mask,fire)]
 103         lon_nofire = lon[np.logical_and(mask,nofire)]
 104         lat_nofire = lat[np.logical_and(mask,nofire)]
 105         # create general arrays
 106         lons = np.concatenate((lon_nofire,lon_fire))
 107         lats = np.concatenate((lat_nofire,lat_fire))
 108         fires = np.concatenate((5*np.ones(lon_nofire.shape),9*np.ones(lon_fire.shape)))
 109
 110         # plot results
 111         if plot:
 112             plt.ion()
 113             plt.plot(lons[fires==5],lats[fires==5],'g.')
 114             plt.plot(lons[fires==9],lats[fires==9],'r.')
 115             plt.show()
 116             plt.pause(.001)
 117             plt.cla()
 118
 119         # update perimeters dictionary
 120         forecast.update({prefix + time_data: Dict({'file': wrfout_file, 'time_tign_g': T, 'lon': lons, 'lat': lats,
 121                             'fire': fires, 'conf_fire': np.array(conf_fire*np.ones(lons[fires==9].shape)),
 122                             'lon_fire': lons[fires==9], 'lat_fire': lats[fires==9], 'lon_nofire': lons[fires==5], 'lat_nofire': lats[fires==5],
 123                             'scan_fire': scan*np.ones(lons[fires==9].shape), 'track_fire': track*np.ones(lons[fires==9].shape),
 124                             'scan_nofire': scan*np.ones(lons[fires==5].shape), 'track_nofire': track*np.ones(lons[fires==9].shape),
 125                             'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 126
 127     return forecast
 128
 129
 130 if __name__ == "__main__":
 131     plot = True
 132     bounds = (-113.85068, -111.89413, 39.677563, 41.156837)
 133     dst = './patch/wrfout_patch'
 134
 135     f = process_forecast(dst,bounds,plot=plot)
 136     sl.save(f,'forecast')