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):
  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     # prefix of the elements in the dictionary
  24     prefix = "IGN"
  25     # initializing dictionary
  26     ignitions = dict({})
  27     # scan and track dimensions of the observation (in km)
  28     scan = 1.
  29     track = 1.
  30
  31     # for each ignition
  32     for lon, lat, time_iso in zip(igns[0],igns[1],igns[2]):
  33         try:
  34             # take coordinates
  35             lon = np.array(lon)
  36             lat = np.array(lat)
  37             # look if coordinates in bounding box
  38             mask = np.logical_and(np.logical_and(np.logical_and(lon>bounds[0],lon<bounds[1]),lat>bounds[2]),lat<bounds[3])
  39             if not mask.sum():
  40                 break
  41             lons = lon[mask]
  42             lats = lat[mask]
  43             # get time elements
  44             time_num = time_iso2num(time_iso)
  45             time_datetime = time_iso2datetime(time_iso)
  46             time_data = '_A%04d%03d_%02d%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
  47                                         time_datetime.hour, time_datetime.minute)
  48             acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
  49             acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
  50         except Exception as e:
  51             print 'Error: bad ignition %s specified.' % igns
  52             print 'Exception: %s.' % e
  53             sys.exit(1)
  54
  55         # no nofire detection
  56         lon_nofire = np.array([])
  57         lat_nofire = np.array([])
  58
  59         # update ignitions dictionary
  60         ignitions.update({prefix + time_data: Dict({'lon': lons, 'lat': lats,
  61                                 'fire': np.array(9*np.ones(lons.shape)), 'conf_fire': np.array(100*np.ones(lons.shape)),
  62                                 'lon_fire': lons, 'lat_fire': lats, 'lon_nofire': lon_nofire, 'lat_nofire': lat_nofire,
  63                                 'scan_fire': scan*np.ones(lons.shape), 'track_fire': track*np.ones(lons.shape),
  64                                 'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
  65     return ignitions
  66
  67
  68 def process_infrared_perimeters(dst,bounds,maxp=1000,ngrid=50,plot=False):
  69     """
  70     Process infrared perimeters the same way than satellite data.
  71
  72     :param dst: path to kml perimeter files
  73     :param bounds: coordinate bounding box filtering to
  74     :param maxp: optional, maximum number of points for each perimeter
  75     :param ngrid: optional, number of nodes for the grid of in/out nodes at each direction
  76     :param plot: optional, boolean to plot or not the result at each perimeter iteration
  77     :return perimeters: dictionary with all the information from each perimeter similar to satellite data
  78
  79     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  80     Angel Farguell (angel.farguell@gmail.com), 2019-05-29
  81     """
  82
  83     # list of kml files in 'dst' path
  84     files = glob.glob(osp.join(dst, '*.kml'))
  85     # prefix of the elements in the dictionary
  86     prefix = "PER"
  87     # initializing dictionary
  88     perimeters = Dict({})
  89     # scan and track dimensions of the observation (in km)
  90     scan = .01
  91     track = .01
  92
  93     # Creating grid where to evaluate in/out of the perimeter
  94     [X,Y] = np.meshgrid(np.linspace(bounds[0],bounds[1],ngrid),np.linspace(bounds[2],bounds[3],ngrid))
  95     XX = np.c_[(X.ravel(),Y.ravel())]
  96
  97     # if any file
  98     if files:
  99         # for each file
 100         for file in files:
 101             print 'Retrieving perimeters from %s' % file
 102             try:
 103                 # open the file
 104                 f = open(file,"r")
 105                 # read all the lines of the file
 106                 f_str = ''.join(f.readlines())
 107                 # close the file
 108                 f.close()
 109             except Exception as e:
 110                 print 'Error: exception when opening file %s, %s' % (file,e.message)
 111                 sys.exit(1)
 112             try:
 113                 # Read name and get time elements
 114                 # read name of the file
 115                 name = re.findall(r'<name>(.*?)</name>',f_str,re.DOTALL)[0]
 116                 # read date of the perimeter
 117                 date = re.match(r'.*([0-9]{2}-[0-9]{2}-[0-9]{4} [0-9]{4})',name).groups()[0]
 118                 # create ISO time of the date
 119                 time_iso = date[6:10]+'-'+date[0:2]+'-'+date[3:5]+'T'+date[11:13]+':'+date[13:15]+':00'
 120                 # create numerical time from the ISO time
 121                 time_num = time_iso2num(time_iso)
 122                 # create datetime element from the ISO time
 123                 time_datetime = time_iso2datetime(time_iso)
 124                 # create time stamp
 125                 time_data = '_A%04d%03d_%02d%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
 126                                                     time_datetime.hour, time_datetime.minute)
 127                 # create acquisition date
 128                 acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
 129                 # create acquisition time
 130                 acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
 131
 132                 # Get the coordinates of all the perimeters
 133                 # regex of the polygons (or perimeters)
 134                 polygons = re.findall(r'<Polygon>(.*?)</Polygon>',f_str,re.DOTALL)
 135                 # for each polygon, regex of the coordinates
 136                 buckets = [re.split('\r\n\s+',re.findall(r'<coordinates>(.*?)</coordinates>',p,re.DOTALL)[0])[1:] for p in polygons]
 137                 # array of arrays with each polygon coordinates
 138                 coordinates = [[np.array(re.split(',',b)[0:2]).astype(float) for b in bucket] for bucket in buckets]
 139             except Exception as e:
 140                 print 'Error: file %s has not proper structure.' % file
 141                 print 'Exception: %s.' % e
 142                 sys.exit(1)
 143
 144             # Plot perimeter
 145             if plot:
 146                 plt.ion()
 147                 plt.plot([coord[0] for coordinate in coordinates for coord in coordinate],[coord[1] for coordinate in coordinates for coord in coordinate],'bx')
 148
 149             # Create upper and lower bound coordinates depending on in/out polygons
 150             # compute path elements for each polygon
 151             paths = [Path(coord) for coord in coordinates]
 152             # compute mask of coordinates inside polygon for each polygon
 153             masks = [path.contains_points(XX) for path in paths]
 154             # logical or for all the masks
 155             inmask = np.logical_or.reduce(masks)
 156             # upper and lower bounds arifitial from in/out polygon
 157             up_arti = XX[inmask]
 158             lon_fire = np.array([up[0] for up in up_arti])
 159             lat_fire = np.array([up[1] for up in up_arti])
 160             low_arti = XX[~inmask]
 161             lon_nofire = np.array([low[0] for low in low_arti])
 162             lat_nofire = np.array([low[1] for low in low_arti])
 163
 164             # take a coarsening of the perimeters
 165             for k,coord in enumerate(coordinates):
 166                 if len(coord) > maxp:
 167                     coarse = len(coord)/maxp
 168                     if coarse > 0:
 169                         coordinates[k] = [coord[ind] for ind in np.concatenate(([0],range(len(coord))[coarse:-coarse:coarse]))]
 170
 171             # append perimeter nodes
 172             lon_fire = np.append(lon_fire,np.array([coord[0] for coordinate in coordinates for coord in coordinate]))
 173             lat_fire = np.append(lat_fire,np.array([coord[1] for coordinate in coordinates for coord in coordinate]))
 174
 175             # create general arrays
 176             lon = np.concatenate((lon_nofire,lon_fire))
 177             lat = np.concatenate((lat_nofire,lat_fire))
 178             fire = np.concatenate((5*np.ones(lon_nofire.shape),9*np.ones(lon_fire.shape)))
 179
 180             # mask in bounds
 181             mask = np.logical_and(np.logical_and(np.logical_and(lon>bounds[0],lon<bounds[1]),lat>bounds[2]),lat<bounds[3])
 182             if not mask.sum():
 183                 break
 184             lons = lon[mask]
 185             lats = lat[mask]
 186             fires = fire[mask]
 187
 188             # plot results
 189             if plot:
 190                 plt.plot(lons[fires==5],lats[fires==5],'g.')
 191                 plt.plot(lons[fires==9],lats[fires==9],'r.')
 192                 plt.show()
 193                 plt.pause(.001)
 194                 plt.cla()
 195
 196             # update perimeters dictionary
 197             perimeters.update({prefix + time_data: Dict({'file': file, 'lon': lons, 'lat': lats,
 198                             'fire': fires, 'conf_fire': np.array(100*np.ones(lons[fires==9].shape)),
 199                             'lon_fire': lons[fires==9], 'lat_fire': lons[fires==9], 'lon_nofire': lats[fires==5], 'lat_nofire': lats[fires==5],
 200                             'scan_fire': scan*np.ones(lons.shape), 'track_fire': track*np.ones(lons.shape),
 201                             'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 202     else:
 203         print 'Warning: No KML files in the path specified'
 204         perimeters = []
 205     return perimeters
 206
 207
 208 if __name__ == "__main__":
 209     plot = True
 210     bounds = (-115.97282409667969, -115.28449249267578, 43.808258056640625, 44.302913665771484)
 211     dst = './pioneer/perim'
 212
 213     p = process_infrared_perimeters(dst,bounds,plot=plot)
 214     print p