infrared_perimeters.py

   1 from matplotlib.path import Path
   2 import matplotlib.pyplot as plt
   3 import numpy as np
   4 import os.path as osp
   5 from JPSSD import time_iso2num, time_iso2datetime
   6 from utils import Dict
   7 import saveload as sl
   8 import re, glob, sys, os
   9
  10
  11 def process_ignitions(igns,bounds,time=None):
  12     """
  13     Process ignitions the same way than satellite data.
  14
  15     :param igns: ([lons],[lats],[dates]) where lons and lats in degrees and dates in ESMF format
  16     :param bounds: coordinate bounding box filtering to
  17     :return ignitions: dictionary with all the information from each ignition similar to satellite data
  18
  19     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  20     Angel Farguell (angel.farguell@gmail.com), 2019-05-29
  21     """
  22
  23     # Time interval
  24     if time:
  25         interval_datetime = map(time_iso2datetime,time)
  26
  27     # prefix of the elements in the dictionary
  28     prefix = "IGN"
  29     # initializing dictionary
  30     ignitions = dict({})
  31     # scan and track dimensions of the observation (in km)
  32     scan = 1.
  33     track = 1.
  34     # confidences
  35     conf_fire = 100
  36
  37     # for each ignition
  38     for lon, lat, time_iso in zip(igns[0],igns[1],igns[2]):
  39         try:
  40             # take coordinates
  41             lon = np.array(lon)
  42             lat = np.array(lat)
  43             # look if coordinates in bounding box
  44             mask = np.logical_and(np.logical_and(np.logical_and(lon >= bounds[0], lon <= bounds[1]),lat >= bounds[2]),lat <= bounds[3])
  45             if not mask.sum():
  46                 continue
  47             lons = lon[mask]
  48             lats = lat[mask]
  49             # get time elements
  50             time_num = time_iso2num(time_iso)
  51             time_datetime = time_iso2datetime(time_iso)
  52             # skip if not inside interval (only if time is determined)
  53             if time and (time_datetime < interval_datetime[0] or time_datetime > interval_datetime[1]):
  54                 print 'Perimeter from %s skipped, not inside the simulation interval!' % file
  55                 continue
  56             time_data = '_A%04d%03d_%02d%02d_%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
  57                                         time_datetime.hour, time_datetime.minute, time_datetime.second)
  58             acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
  59             acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
  60         except Exception as e:
  61             print 'Error: bad ignition %s specified.' % igns
  62             print 'Exception: %s.' % e
  63             sys.exit(1)
  64
  65         # no nofire detection
  66         lon_nofire = np.array([])
  67         lat_nofire = np.array([])
  68         conf_nofire = np.array([])
  69
  70         # update ignitions dictionary
  71         ignitions.update({prefix + time_data: Dict({'lon': lons, 'lat': lats,
  72                                 'fire': np.array(9*np.ones(lons.shape)), 'conf_fire': np.array(conf_fire*np.ones(lons.shape)),
  73                                 'lon_fire': lons, 'lat_fire': lats, 'lon_nofire': lon_nofire, 'lat_nofire': lat_nofire,
  74                                 'scan_fire': scan*np.ones(lons.shape), 'track_fire': track*np.ones(lons.shape),
  75                                 'conf_nofire' : conf_nofire, 'scan_nofire': scan*np.ones(lon_nofire.shape),
  76                                 'track_nofire': track*np.ones(lon_nofire.shape), 'time_iso': time_iso,
  77                                 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
  78     return ignitions
  79
  80
  81 def process_infrared_perimeters(dst,bounds,time=None,maxp=1000,ngrid=100,plot=False):
  82     """
  83     Process infrared perimeters the same way than satellite data.
  84
  85     :param dst: path to kml perimeter files
  86     :param bounds: coordinate bounding box filtering to
  87     :param time: optional, time interval in ISO
  88     :param maxp: optional, maximum number of points for each perimeter
  89     :param ngrid: optional, number of nodes for the grid of in/out nodes at each direction
  90     :param plot: optional, boolean to plot or not the result at each perimeter iteration
  91     :return perimeters: dictionary with all the information from each perimeter similar to satellite data
  92
  93     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  94     Angel Farguell (angel.farguell@gmail.com), 2019-05-29
  95     """
  96
  97     # Time interval
  98     if time:
  99         interval_datetime = map(time_iso2datetime,time)
 100
 101     # list of kml files in 'dst' path
 102     files = glob.glob(osp.join(dst, '*.kml'))
 103     # prefix of the elements in the dictionary
 104     prefix = "PER"
 105     # initializing dictionary
 106     perimeters = Dict({})
 107     # scan and track dimensions of the observation (in km)
 108     scan = .5
 109     track = .5
 110     # confidences
 111     conf_fire = 100
 112     conf_nofire = 100
 113
 114     # Creating grid where to evaluate in/out of the perimeter
 115     [X,Y] = np.meshgrid(np.linspace(bounds[0],bounds[1],ngrid),np.linspace(bounds[2],bounds[3],ngrid))
 116     XX = np.c_[(X.ravel(),Y.ravel())]
 117
 118     # if any file
 119     if files:
 120         # for each file
 121         for file in files:
 122             print 'Retrieving perimeters from %s' % file
 123             try:
 124                 # open the file
 125                 f = open(file,"r")
 126                 # read all the lines of the file
 127                 f_str = ''.join(f.readlines())
 128                 # close the file
 129                 f.close()
 130             except Exception as e:
 131                 print 'Error: exception when opening file %s, %s' % (file,e.message)
 132                 sys.exit(1)
 133             try:
 134                 # Read name and get time elements
 135                 # read name of the file
 136                 name = re.findall(r'<name>(.*?)</name>',f_str,re.DOTALL)[0]
 137                 # read date of the perimeter
 138                 date = re.match(r'.*([0-9]+)-([0-9]+)-([0-9]+) ([0-9]{2})([0-9]{2})',name).groups()
 139                 date = (date[2],date[0],date[1],date[3],date[4])
 140                 # create ISO time of the date
 141                 time_iso = '%04d-%02d-%02dT%02d:%02d:00' % tuple([ int(d) for d in date ])
 142                 # create numerical time from the ISO time
 143                 time_num = time_iso2num(time_iso)
 144                 # create datetime element from the ISO time
 145                 time_datetime = time_iso2datetime(time_iso)
 146                 # skip if not inside interval (only if time is determined)
 147                 if time and (time_datetime < interval_datetime[0] or time_datetime > interval_datetime[1]):
 148                     print 'Perimeter from %s skipped, not inside the simulation interval!' % file
 149                     continue
 150                 # create time stamp
 151                 time_data = '_A%04d%03d_%02d%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
 152                                                     time_datetime.hour, time_datetime.minute)
 153                 # create acquisition date
 154                 acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
 155                 # create acquisition time
 156                 acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
 157
 158                 # Get the coordinates of all the perimeters
 159                 # regex of the polygons (or perimeters)
 160                 polygons = re.findall(r'<Polygon>(.*?)</Polygon>',f_str,re.DOTALL)
 161                 # for each polygon, regex of the coordinates
 162                 buckets = [re.split('\r\n\s+',re.findall(r'<coordinates>(.*?)</coordinates>',p,re.DOTALL)[0])[1:] for p in polygons]
 163                 # array of arrays with each polygon coordinates
 164                 coordinates = [[np.array(re.split(',',b)[0:2]).astype(float) for b in bucket] for bucket in buckets]
 165             except Exception as e:
 166                 print 'Error: file %s has not proper structure.' % file
 167                 print 'Exception: %s.' % e
 168                 sys.exit(1)
 169
 170             # Plot perimeter
 171             if plot:
 172                 plt.ion()
 173                 plt.plot([coord[0] for coordinate in coordinates for coord in coordinate],[coord[1] for coordinate in coordinates for coord in coordinate],'bx')
 174
 175             # Create upper and lower bound coordinates depending on in/out polygons
 176             # compute path elements for each polygon
 177             paths = [Path(coord) for coord in coordinates]
 178             # compute mask of coordinates inside polygon for each polygon
 179             masks = [path.contains_points(XX) for path in paths]
 180             # logical or for all the masks
 181             inmask = np.logical_or.reduce(masks)
 182             # upper and lower bounds arifitial from in/out polygon
 183             up_arti = XX[inmask]
 184             lon_fire = np.array([up[0] for up in up_arti])
 185             lat_fire = np.array([up[1] for up in up_arti])
 186             low_arti = XX[~inmask]
 187             lon_nofire = np.array([low[0] for low in low_arti])
 188             lat_nofire = np.array([low[1] for low in low_arti])
 189
 190             # take a coarsening of the perimeters
 191             for k,coord in enumerate(coordinates):
 192                 if len(coord) > maxp:
 193                     coarse = len(coord)/maxp
 194                     if coarse > 0:
 195                         coordinates[k] = [coord[ind] for ind in np.concatenate(([0],range(len(coord))[coarse:-coarse:coarse]))]
 196
 197             # append perimeter nodes
 198             lon_fire = np.append(lon_fire,np.array([coord[0] for coordinate in coordinates for coord in coordinate]))
 199             lat_fire = np.append(lat_fire,np.array([coord[1] for coordinate in coordinates for coord in coordinate]))
 200
 201             # create general arrays
 202             lon = np.concatenate((lon_nofire,lon_fire))
 203             lat = np.concatenate((lat_nofire,lat_fire))
 204             fire = np.concatenate((5*np.ones(lon_nofire.shape),9*np.ones(lon_fire.shape)))
 205
 206             # mask in bounds
 207             mask = np.logical_and(np.logical_and(np.logical_and(lon >= bounds[0],lon <= bounds[1]),lat >= bounds[2]),lat <= bounds[3])
 208             if not mask.sum():
 209                 break
 210             lons = lon[mask]
 211             lats = lat[mask]
 212             fires = fire[mask]
 213
 214             # plot results
 215             if plot:
 216                 plt.plot(lons[fires==5],lats[fires==5],'g.')
 217                 plt.plot(lons[fires==9],lats[fires==9],'r.')
 218                 plt.show()
 219                 plt.pause(.5)
 220                 plt.cla()
 221
 222             # update perimeters dictionary
 223             perimeters.update({prefix + time_data: Dict({'file': file, 'lon': lons, 'lat': lats,
 224                             'fire': fires, 'conf_fire': np.array(conf_fire*np.ones(lons[fires==9].shape)),
 225                             'lon_fire': lons[fires==9], 'lat_fire': lats[fires==9], 'lon_nofire': lons[fires==5], 'lat_nofire': lats[fires==5],
 226                             'scan_fire': scan*np.ones(lons[fires==9].shape), 'track_fire': track*np.ones(lons[fires==9].shape),
 227                             'conf_nofire': np.array(conf_nofire*np.ones(lons[fires==5].shape)),
 228                             'scan_nofire': scan*np.ones(lons[fires==5].shape), 'track_nofire': track*np.ones(lons[fires==9].shape),
 229                             'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 230     else:
 231         print 'Warning: No KML files in the path specified'
 232         perimeters = []
 233     return perimeters
 234
 235
 236 if __name__ == "__main__":
 237     # Experiment to do
 238     exp = 1
 239     # Plot perimeters as created
 240     plot = True
 241
 242     # Defining options
 243     def pioneer():
 244         bounds = (-115.97282409667969, -115.28449249267578, 43.808258056640625, 44.302913665771484)
 245         time_iso = ('2016-07-18T00:00:00Z', '2016-08-31T00:00:00Z')
 246         igns = None
 247         perims = './pioneer/perim'
 248         return bounds, time_iso, igns, perims
 249     def patch():
 250         bounds = (-113.85068, -111.89413, 39.677563, 41.156837)
 251         time_iso = ('2013-08-10T00:00:00Z', '2013-08-15T00:00:00Z')
 252         igns = ([-112.676039],[40.339372],['2013-08-10T20:00:00Z'])
 253         perims = './patch/perim'
 254         return bounds, time_iso, igns, perims
 255
 256     # Creating the options
 257     options = {1: pioneer, 2: patch}
 258
 259     # Defining the option depending on the experiment
 260     bounds, time_iso, igns, perims = options[exp]()
 261
 262     # Processing infrared perimeters
 263     p = process_infrared_perimeters(perims,bounds,time=time_iso,plot=plot)
 264     # Processing ignitions if defined
 265     if igns:
 266         p.update(process_ignitions(igns,bounds,time=time_iso))
 267     # Saving results
 268     sl.save(p,'perimeters')
 269