forecast.py

   1 import matplotlib.pyplot as plt
   2 from mpl_toolkits.mplot3d import axes3d
   3 import matplotlib.colors as colors
   4 import numpy as np
   5 from datetime import timedelta
   6 from JPSSD import time_iso2num, time_iso2datetime, time_datetime2iso
   7 from utils import Dict
   8 import re, glob, sys, os
   9
  10 def process_tign_g(lon,lat,tign_g,ctime,scale,time_num,epsilon=5,plot=False):
  11     """
  12     Process forecast from lon, lat, and tign_g as 3D data points
  13
  14     :param lon: array of longitudes
  15     :param lat: array of latitudes
  16     :param tign_g: array of fire arrival time
  17     :param ctime: time char in wrfout of the last fire arrival time
  18     :param scale: numerical time scale in seconds of start and end of simulation
  19     :param time_num: numerical time interval in seconds of start and end of simulation
  20     :param epsilon: optional, epsilon in seconds to define the separation of artificial data
  21     :param plot: optional, plotting results
  22
  23     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  24     Angel Farguell (angel.farguell@gmail.com) and James Haley, 2019-06-05
  25     """
  26
  27     # confidences
  28     conf_ground = 25.
  29     conf_fire = 25.
  30
  31     # ctime transformations
  32     ctime_iso = ctime.replace('_','T')
  33     ctime_datetime = time_iso2datetime(ctime_iso)
  34     tscale = 3600.*24.
  35
  36     # tign_g as data points smaller than maximum tign_g
  37     t1d = np.ravel(tign_g)
  38     etime = t1d.max()
  39     mask = t1d < etime
  40     x1d = np.ravel(lon)[mask]
  41     y1d = np.ravel(lat)[mask]
  42     t1d = t1d[mask]
  43
  44     # compute time as days from start of the simulation
  45     tt = (np.array([time_iso2num(time_datetime2iso(
  46                         ctime_datetime-timedelta(seconds=float(etime-T))))
  47                             for T in t1d])-scale[0])/tscale
  48
  49     # define forecast artificial data points
  50     fground = np.c_[x1d.ravel(),y1d.ravel(),tt.ravel() - epsilon/tscale]
  51     ffire = np.c_[x1d.ravel(),y1d.ravel(),tt.ravel() + epsilon/tscale]
  52
  53     # coarsening
  54     maxsize = 5e3
  55     coarsening = np.int(1+len(ffire)/maxsize)
  56     fground = fground[::coarsening,:]
  57     ffire = ffire[::coarsening,:]
  58
  59     # define output data
  60     X = np.concatenate((fground,ffire))
  61     y = np.concatenate((-np.ones(len(fground)),np.ones(len(ffire))))
  62     c = np.concatenate((conf_ground*np.ones(len(fground)),conf_fire*np.ones(len(ffire))))
  63
  64     # plot results
  65     if plot:
  66         from contline import get_contour_verts
  67         from contour2kml import contour2kml
  68         tt = np.array([time_iso2num(time_datetime2iso(
  69                         ctime_datetime-timedelta(seconds=float(etime-T))))
  70                             for T in np.ravel(tign_g)])
  71         data = get_contour_verts(lon, lat, np.reshape(tt,tign_g.shape), time_num, contour_dt_hours=6, contour_dt_init=24, contour_dt_final=12)
  72         contour2kml(data,'perimeters_forecast.kml')
  73         col = [(0, .5, 0), (.5, 0, 0)]
  74         cm_GR = colors.LinearSegmentedColormap.from_list('GrRd',col,N=2)
  75         fig = plt.figure()
  76         ax = fig.gca(projection='3d')
  77         ax.scatter(X[:, 0], X[:, 1], X[:, 2],
  78                     c=y, cmap=cm_GR, s=1, alpha=.5,
  79                     vmin=y.min(), vmax=y.max())
  80         ax.set_xlabel("Longitude")
  81         ax.set_ylabel("Latitude")
  82         ax.set_zlabel("Fire arrival time [days]")
  83         plt.show()
  84
  85     return X,y,c
  86
  87
  88 def process_tign_g_slices(lon,lat,tign_g,bounds,ctime,dx,dy,wrfout_file='',dt_for=3600.,plot=False):
  89     """
  90     Process forecast as slices on time from lon, lat, and tign_g
  91
  92     :param lon: array of longitudes
  93     :param lat: array of latitudes
  94     :param tign_g: array of fire arrival time
  95     :param bounds: coordinate bounding box filtering to
  96     :param ctime: time char in wrfout of the last fire arrival time
  97     :param dx,dy: data resolution in meters
  98     :param dt_for: optional, temporal resolution to get the fire arrival time from in seconds
  99     :param wrfout_file: optional, to get the name of the file in the dictionary
 100
 101     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 102     Angel Farguell (angel.farguell@gmail.com) and James Haley, 2019-06-05
 103     """
 104
 105     if plot:
 106         fig = plt.figure()
 107
 108     # prefix of the elements in the dictionary
 109     prefix = "FOR"
 110     # confidences
 111     conf_fire = 50
 112     conf_nofire = 10
 113     # margin percentage
 114     margin = .5
 115     # scan and track dimensions of the observation (in km)
 116     scan = dx/1000.
 117     track = dy/1000.
 118     # initializing dictionary
 119     forecast = Dict({})
 120
 121     # ctime transformations
 122     ctime_iso = ctime.replace('_','T')
 123     ctime_datetime = time_iso2datetime(ctime_iso)
 124     # mask coordinates to bounding box
 125     mask = np.logical_and(np.logical_and(np.logical_and(lon >= bounds[0], lon <= bounds[1]), lat >= bounds[2]), lat <= bounds[3])
 126     # create a square subset of fire arrival time less than the maximum
 127     mtign = np.logical_and(mask, tign_g < tign_g.max())
 128     mlon = lon[mtign]
 129     mlat = lat[mtign]
 130     mlenlon = margin*(mlon.max()-mlon.min())
 131     mlenlat = margin*(mlat.max()-mlat.min())
 132     sbounds = (mlon.min()-mlenlon, mlon.max()+mlenlon, mlat.min()-mlenlat, mlat.max()+mlenlat)
 133     smask = np.logical_and(np.logical_and(np.logical_and(lon >= sbounds[0], lon <= sbounds[1]), lat >= sbounds[2]), lat <= sbounds[3])
 134
 135     # times to get fire arrival time from
 136     TT = np.arange(tign_g.min()-3*dt_for,tign_g.max(),dt_for)[0:-1]
 137     for T in TT:
 138         # fire arrival time to datetime
 139         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)
 140         # create ISO time of the date
 141         time_iso = time_datetime2iso(time_datetime)
 142         # create numerical time from the ISO time
 143         time_num = time_iso2num(time_iso)
 144         # create time stamp
 145         time_data = '_A%04d%03d_%02d%02d_%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
 146                                                 time_datetime.hour, time_datetime.minute, time_datetime.second)
 147         # create acquisition date
 148         acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
 149         # create acquisition time
 150         acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
 151
 152         print 'Retrieving forecast from %s' % time_iso
 153
 154         # masks
 155         fire = tign_g <= T
 156         nofire = np.logical_and(tign_g > T, np.logical_or(tign_g != tign_g.max(), smask))
 157         # coordinates masked
 158         lon_fire = lon[np.logical_and(mask,fire)]
 159         lat_fire = lat[np.logical_and(mask,fire)]
 160         lon_nofire = lon[np.logical_and(mask,nofire)]
 161         lat_nofire = lat[np.logical_and(mask,nofire)]
 162         # create general arrays
 163         lons = np.concatenate((lon_nofire,lon_fire))
 164         lats = np.concatenate((lat_nofire,lat_fire))
 165         fires = np.concatenate((5*np.ones(lon_nofire.shape),9*np.ones(lon_fire.shape)))
 166
 167         # plot results
 168         if plot:
 169             plt.ion()
 170             plt.cla()
 171             plt.plot(lons[fires==5],lats[fires==5],'g.')
 172             plt.plot(lons[fires==9],lats[fires==9],'r.')
 173             plt.pause(.0001)
 174             plt.show()
 175
 176         # update perimeters dictionary
 177         forecast.update({prefix + time_data: Dict({'file': wrfout_file, 'time_tign_g': T, 'lon': lons, 'lat': lats,
 178                             'fire': fires, 'conf_fire': np.array(conf_fire*np.ones(lons[fires==9].shape)),
 179                             'lon_fire': lons[fires==9], 'lat_fire': lats[fires==9], 'lon_nofire': lons[fires==5], 'lat_nofire': lats[fires==5],
 180                             'scan_fire': scan*np.ones(lons[fires==9].shape), 'track_fire': track*np.ones(lons[fires==9].shape),
 181                             'conf_nofire': np.array(conf_nofire*np.ones(lons[fires==5].shape)),
 182                             'scan_nofire': scan*np.ones(lons[fires==5].shape), 'track_nofire': track*np.ones(lons[fires==9].shape),
 183                             'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 184
 185     return forecast
 186
 187 def read_forecast_wrfout(wrfout_file):
 188     """
 189     Read forecast from a wrfout.
 190
 191     :param wrfout_file: path to wrfout file
 192
 193     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 194     Angel Farguell (angel.farguell@gmail.com), 2019-06-05
 195     """
 196
 197     # read file
 198     import netCDF4 as nc
 199     try:
 200         data = nc.Dataset(wrfout_file)
 201     except Exception as e:
 202         print 'Warning: No netcdf file %s in the path' % wrfout_file
 203         return []
 204     # current time
 205     ctime = ''.join(data['Times'][0])
 206     # getting rid of strip
 207     atmlenx = len(data.dimensions['west_east'])
 208     atmleny = len(data.dimensions['south_north'])
 209     staglenx = len(data.dimensions['west_east_stag'])
 210     stagleny = len(data.dimensions['south_north_stag'])
 211     dimlenx = len(data.dimensions['west_east_subgrid'])
 212     dimleny = len(data.dimensions['south_north_subgrid'])
 213     ratiox = dimlenx/max(staglenx,atmlenx+1)
 214     ratioy = dimleny/max(stagleny,atmleny+1)
 215     lenx = dimlenx-ratiox
 216     leny = dimleny-ratioy
 217     # coordinates
 218     lon = data['FXLONG'][0][0:lenx,0:leny]
 219     lat = data['FXLAT'][0][0:lenx,0:leny]
 220     # fire arrival time
 221     tign_g = data['TIGN_G'][0][0:lenx,0:leny]
 222     # resolutions
 223     DX = data.getncattr('DX')
 224     DY = data.getncattr('DY')
 225     sx = data.dimensions['west_east_subgrid'].size/data.dimensions['west_east_stag'].size
 226     sy = data.dimensions['south_north_subgrid'].size/data.dimensions['south_north_stag'].size
 227     dx = DX/sx
 228     dy = DY/sy
 229     # close netcdf file
 230     data.close()
 231     return lon,lat,tign_g,ctime,dx,dy
 232
 233 def process_forecast_slides_wrfout(wrfout_file,bounds,plot=False):
 234     """
 235     Process forecast from a wrfout.
 236
 237     :param wrfout_file: path to wrfout file
 238     :param bounds: coordinate bounding box filtering to
 239
 240     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 241     Angel Farguell (angel.farguell@gmail.com), 2019-06-05
 242     """
 243
 244     # read forecast from wrfout
 245     lon,lat,tign_g,ctime,dx,dy = read_forecast_wrfout(wrfout_file)
 246     # create forecast
 247     forecast = process_tign_g_slices(lon,lat,tign_g,bounds,ctime,dx,dy,wrfout_file=wrfout_file,plot=plot)
 248
 249     return forecast
 250
 251 def process_forecast_wrfout(wrfout_file,scale,time_num,epsilon=5,plot=False):
 252     """
 253     Process forecast from a wrfout.
 254
 255     :param wrfout_file: path to wrfout file
 256     :param scale: numerical time in seconds of start and end of simulation
 257
 258     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 259     Angel Farguell (angel.farguell@gmail.com), 2019-06-05
 260     """
 261
 262     # read forecast from wrfout
 263     lon,lat,tign_g,ctime,dx,dy = read_forecast_wrfout(wrfout_file)
 264     # create forecast
 265     X,y,c = process_tign_g(lon,lat,tign_g,ctime,scale,time_num,epsilon=epsilon,plot=plot)
 266     print 'shape X_f: ', X.shape
 267     print 'shape y_f: ', y.shape
 268     print 'shape c_f: ', c.shape
 269     print 'len fire: ', len(X[y==1])
 270     print 'len ground: ', len(X[y==-1])
 271
 272     return X,y,c
 273
 274 if __name__ == "__main__":
 275     import saveload as sl
 276     real = True
 277     cloud = True
 278
 279     if real:
 280         plot = True
 281         if cloud:
 282             time_num = [1376114400.0, 1376546400.0]
 283             scale = [1375898400.0, 1377410400.0]
 284             dst = './patch/wrfout_patch'
 285             X,y,c = process_forecast_wrfout(dst,scale,time_num,plot=plot)
 286             sl.save(dict({'X': X, 'y': y, 'c': c}),'forecast')
 287         else:
 288             bounds = (-113.85068, -111.89413, 39.677563, 41.156837)
 289             dst = './patch/wrfout_patch'
 290             f = process_forecast_slides_wrfout(dst,bounds,plot=plot)
 291             sl.save(f,'forecast')
 292     else:
 293         from infrared_perimeters import process_ignitions
 294         from setup import process_detections
 295         dst = 'ideal_test'
 296         plot = False
 297         ideal = sl.load(dst)
 298         kk = 4
 299         data = process_tign_g_slices(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)
 300         if 'point' in ideal.keys():
 301             p = [[ideal['point'][0]],[ideal['point'][1]],[ideal['point'][2]]]
 302             data.update(process_ignitions(p,ideal['bounds']))
 303         elif 'points' in ideal.keys():
 304             p = [[point[0] for point in ideal['points']],[point[1] for point in ideal['points']],[point[2] for point in ideal['points']]]
 305             data.update(process_ignitions(p,ideal['bounds']))
 306         etime = time_iso2num(ideal['ctime'].replace('_','T'))
 307         time_num_int = (etime-ideal['tign_g'].max(),etime)
 308         sl.save((data,ideal['lon'],ideal['lat'],time_num_int),'data')
 309         process_detections(data,ideal['lon'],ideal['lat'],time_num_int)
 310