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 def process_paths(outer_coords,inner_coords,min_points=50):
  11     outer_paths = []
  12     inner_paths = []
  13     for n,outer in enumerate(outer_coords):
  14         outer = np.array(outer)
  15         if len(outer[:,0]) > min_points:
  16             path = Path(outer)
  17             outer_paths.append(path)
  18             inners = np.array(inner_coords[n])
  19             if len(inners.shape) > 2:
  20                 in_paths = []
  21                 for inner in inners:
  22                     if len(inner) > min_points:
  23                         path = Path(inner)
  24                         in_paths.append(path)
  25             else:
  26                 if len(inners) > min_points:
  27                     in_paths = Path(inners)
  28                 else:
  29                     in_paths = []
  30             inner_paths.append(in_paths)
  31     return outer_paths,inner_paths
  32
  33 def process_ignitions(igns,bounds,time=None):
  34     """
  35     Process ignitions the same way than satellite data.
  36
  37     :param igns: ([lons],[lats],[dates]) where lons and lats in degrees and dates in ESMF format
  38     :param bounds: coordinate bounding box filtering to
  39     :return ignitions: dictionary with all the information from each ignition similar to satellite data
  40
  41     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  42     Angel Farguell (angel.farguell@gmail.com), 2019-05-29
  43     """
  44
  45     # Time interval
  46     if time:
  47         interval_datetime = map(time_iso2datetime,time)
  48
  49     # prefix of the elements in the dictionary
  50     prefix = "IGN"
  51     # initializing dictionary
  52     ignitions = dict({})
  53     # scan and track dimensions of the observation (in km)
  54     scan = 1.
  55     track = 1.
  56     # confidences
  57     conf_fire = 100.
  58
  59     # for each ignition
  60     for lon, lat, time_iso in zip(igns[0],igns[1],igns[2]):
  61         try:
  62             # take coordinates
  63             lon = np.array(lon)
  64             lat = np.array(lat)
  65             # look if coordinates in bounding box
  66             mask = np.logical_and(np.logical_and(np.logical_and(lon >= bounds[0], lon <= bounds[1]),lat >= bounds[2]),lat <= bounds[3])
  67             if not mask.sum():
  68                 continue
  69             lons = lon[mask]
  70             lats = lat[mask]
  71             # get time elements
  72             time_num = time_iso2num(time_iso)
  73             time_datetime = time_iso2datetime(time_iso)
  74             # skip if not inside interval (only if time is determined)
  75             if time and (time_datetime < interval_datetime[0] or time_datetime > interval_datetime[1]):
  76                 print 'Perimeter from %s skipped, not inside the simulation interval!' % file
  77                 continue
  78             time_data = '_A%04d%03d_%02d%02d_%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
  79                                         time_datetime.hour, time_datetime.minute, time_datetime.second)
  80             acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
  81             acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
  82         except Exception as e:
  83             print 'Error: bad ignition %s specified.' % igns
  84             print 'Exception: %s.' % e
  85             sys.exit(1)
  86
  87         # no nofire detection
  88         lon_nofire = np.array([])
  89         lat_nofire = np.array([])
  90         conf_nofire = np.array([])
  91
  92         # update ignitions dictionary
  93         ignitions.update({prefix + time_data: Dict({'lon': lons, 'lat': lats,
  94                                 'fire': np.array(9*np.ones(lons.shape)), 'conf_fire': np.array(conf_fire*np.ones(lons.shape)),
  95                                 'lon_fire': lons, 'lat_fire': lats, 'lon_nofire': lon_nofire, 'lat_nofire': lat_nofire,
  96                                 'scan_fire': scan*np.ones(lons.shape), 'track_fire': track*np.ones(lons.shape),
  97                                 'conf_nofire' : conf_nofire, 'scan_nofire': scan*np.ones(lon_nofire.shape),
  98                                 'track_nofire': track*np.ones(lon_nofire.shape), 'time_iso': time_iso,
  99                                 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 100     return ignitions
 101
 102
 103 def process_infrared_perimeters(dst,bounds,time=None,maxp=500,ngrid=50,plot=False,gen_polys=False):
 104     """
 105     Process infrared perimeters the same way than satellite data.
 106
 107     :param dst: path to kml perimeter files
 108     :param bounds: coordinate bounding box filtering to
 109     :param time: optional, time interval in ISO
 110     :param maxp: optional, maximum number of points for each perimeter
 111     :param ngrid: optional, number of nodes for the grid of in/out nodes at each direction
 112     :param plot: optional, boolean to plot or not the result at each perimeter iteration
 113     :return perimeters: dictionary with all the information from each perimeter similar to satellite data
 114
 115     Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 116     Angel Farguell (angel.farguell@gmail.com), 2019-05-29
 117     """
 118
 119     # Time interval
 120     if time:
 121         interval_datetime = map(time_iso2datetime,time)
 122
 123     # list of kml files in 'dst' path
 124     files = glob.glob(osp.join(dst, '*.kml'))
 125     # prefix of the elements in the dictionary
 126     prefix = "PER"
 127     # initializing dictionary
 128     perimeters = Dict({})
 129     # scan and track dimensions of the observation (in km)
 130     scan = .5
 131     track = .5
 132     # confidences
 133     conf_fire = 100.
 134     conf_nofire = 70.
 135
 136     # Creating grid where to evaluate in/out of the perimeter
 137     [X,Y] = np.meshgrid(np.linspace(bounds[0],bounds[1],ngrid),np.linspace(bounds[2],bounds[3],ngrid))
 138     XX = np.c_[np.ravel(X),np.ravel(Y)]
 139
 140     # if any file
 141     if files:
 142         # for each file
 143         for file in files:
 144             print 'Retrieving perimeters from %s' % file
 145             try:
 146                 # open the file
 147                 f = open(file,"r")
 148                 # read all the lines of the file
 149                 f_str = ''.join(f.readlines())
 150                 # close the file
 151                 f.close()
 152             except Exception as e:
 153                 print 'Error: exception when opening file %s, %s' % (file,e.message)
 154                 sys.exit(1)
 155             try:
 156                 # Read name and get time elements
 157                 # read name of the file
 158                 name = re.findall(r'<name>(.*?)</name>',f_str,re.DOTALL)[0]
 159                 # read date of the perimeter
 160                 date = re.match(r'.* ([0-9]+)-([0-9]+)-([0-9]+) ([0-9]{2})([0-9]{2})',name).groups()
 161                 date = (date[2],date[0],date[1],date[3],date[4])
 162                 # create ISO time of the date
 163                 time_iso = '%04d-%02d-%02dT%02d:%02d:00' % tuple([ int(d) for d in date ])
 164                 # create numerical time from the ISO time
 165                 time_num = time_iso2num(time_iso)
 166                 # create datetime element from the ISO time
 167                 time_datetime = time_iso2datetime(time_iso)
 168                 # skip if not inside interval (only if time is determined)
 169                 if time and (time_datetime < interval_datetime[0] or time_datetime > interval_datetime[1]):
 170                     print 'Perimeter from %s skipped, not inside the simulation interval!' % file
 171                     continue
 172                 # create time stamp
 173                 time_data = '_A%04d%03d_%02d%02d' % (time_datetime.year, time_datetime.timetuple().tm_yday,
 174                                                     time_datetime.hour, time_datetime.minute)
 175                 # create acquisition date
 176                 acq_date = '%04d-%02d-%02d' % (time_datetime.year, time_datetime.month, time_datetime.day)
 177                 # create acquisition time
 178                 acq_time = '%02d%02d' % (time_datetime.hour, time_datetime.minute)
 179
 180                 # Get the coordinates of all the perimeters
 181                 # regex of the polygons (or perimeters)
 182                 polygons = re.findall(r'<Polygon>(.*?)</Polygon>',f_str,re.DOTALL)
 183                 # for each polygon, outer boundary
 184                 outer_buckets = [re.findall(r'<outerBoundaryIs>(.*?)</outerBoundaryIs>',p,re.DOTALL)[0] for p in polygons]
 185                 # for each outer polygon, regex of the coordinates
 186                 buckets = [re.split(r'\s',re.findall(r'<coordinates>(.*?)</coordinates>',p,re.DOTALL)[0])[1:] for p in outer_buckets]
 187                 # array of arrays with each outer polygon coordinates
 188                 outer_coordinates = [[np.array(re.split(',',b)[0:2]).astype(float) for b in bucket if b is not ''] for bucket in buckets]
 189                 # for each polygon, inner boundary
 190                 inner_buckets = [re.findall(r'<innerBoundaryIs>(.*?)</innerBoundaryIs>',p,re.DOTALL)[0] if re.findall(r'<innerBoundaryIs>(.*?)</innerBoundaryIs>',p,re.DOTALL) else '' for p in polygons]
 191                 # for each inner polygon, regex of the coordinates
 192                 buckets = [re.split(r'\s',re.findall(r'<coordinates>(.*?)</coordinates>',p,re.DOTALL)[0])[1:] if p != '' else '' for p in inner_buckets]
 193                 # array of arrays with each inner polygon coordinates
 194                 inner_coordinates = [[np.array(re.split(',',b)[0:2]).astype(float) for b in bucket if b is not ''] for bucket in buckets]
 195             except Exception as e:
 196                 print 'Error: file %s has not proper structure.' % file
 197                 print 'Exception: %s.' % e
 198                 sys.exit(1)
 199
 200             # Plot perimeter
 201             if plot:
 202                 plt.ion()
 203                 for outer in outer_coordinates:
 204                     if len(outer):
 205                         x = np.array(outer)[:,0]
 206                         y = np.array(outer)[:,1]
 207                         plt.plot(x,y,'gx')
 208                 for inner in inner_coordinates:
 209                     if len(inner):
 210                         x = np.array(outer)[:,0]
 211                         y = np.array(outer)[:,1]
 212                         plt.plot(x,y,'rx')
 213
 214             # Create paths for each polygon (outer and inner)
 215             outer_paths,inner_paths = process_paths(outer_coordinates,inner_coordinates)
 216             if len(outer_paths):
 217                 # compute mask of coordinates inside outer polygons
 218                 outer_mask = np.logical_or.reduce([path.contains_points(XX) for path in outer_paths if path])
 219                 # compute mask of coordinates inside inner polygons
 220                 inner_mask = np.logical_or.reduce([path.contains_points(XX) for path in inner_paths if path])
 221                 # mask inside outer polygons and outside inner polygons
 222                 inmask = outer_mask * ~inner_mask
 223                 # upper and lower bounds arifitial from in/out polygon
 224                 up_arti = XX[inmask]
 225                 lon_fire = np.array([up[0] for up in up_arti])
 226                 lat_fire = np.array([up[1] for up in up_arti])
 227                 low_arti = XX[~inmask]
 228                 lon_nofire = np.array([low[0] for low in low_arti])
 229                 lat_nofire = np.array([low[1] for low in low_arti])
 230
 231                 # take a coarsening of the outer polygons
 232                 coordinates = []
 233                 for k,coord in enumerate(outer_coordinates):
 234                     if len(coord) > maxp:
 235                         coarse = len(coord)/maxp
 236                         if coarse > 0:
 237                             coordinates += [coord[ind] for ind in np.concatenate(([0],range(len(coord))[coarse:-coarse:coarse]))]
 238
 239                 if coordinates:
 240                     # append perimeter nodes
 241                     lon_fire = np.append(lon_fire,np.array(coordinates)[:,0])
 242                     lat_fire = np.append(lat_fire,np.array(coordinates)[:,1])
 243
 244                 # create general arrays
 245                 lon = np.concatenate((lon_nofire,lon_fire))
 246                 lat = np.concatenate((lat_nofire,lat_fire))
 247                 fire = np.concatenate((5*np.ones(lon_nofire.shape),9*np.ones(lon_fire.shape)))
 248
 249                 # mask in bounds
 250                 mask = np.logical_and(np.logical_and(np.logical_and(lon >= bounds[0],lon <= bounds[1]),lat >= bounds[2]),lat <= bounds[3])
 251                 if not mask.sum():
 252                     break
 253                 lons = lon[mask]
 254                 lats = lat[mask]
 255                 fires = fire[mask]
 256
 257                 # plot results
 258                 if plot:
 259                     plt.plot(lons[fires==5],lats[fires==5],'g.')
 260                     plt.plot(lons[fires==9],lats[fires==9],'r.')
 261                     plt.show()
 262                     plt.pause(.5)
 263                     plt.cla()
 264             else:
 265                 lons = np.array([])
 266                 lats = np.array([])
 267                 fires = np.array([])
 268
 269             if gen_polys:
 270                 from shapely.geometry import Polygon
 271                 from shapely.ops import transform
 272                 from functools import partial
 273                 import pyproj
 274
 275                 proj4_wgs84 = '+proj=latlong +ellps=WGS84 +datum=WGS84 +units=degree +no_defs'
 276                 proj4_moll = '+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
 277                 proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'),
 278                                 pyproj.Proj(proj4_moll))
 279                 polys = []
 280                 for n,outer in enumerate(outer_coordinates):
 281                     x = np.array(outer)[:,0]
 282                     y = np.array(outer)[:,1]
 283                     shape = Polygon([(l[0],l[1]) for l in zip(x,y)]).buffer(0)
 284                     inner = np.array(inner_coordinates[n])
 285                     if len(inner):
 286                         if len(inner.shape) > 2:
 287                             for h in inner_coordinates[n]:
 288                                 xh = np.array(h)[:,0]
 289                                 yh = np.array(h)[:,1]
 290                         else:
 291                             xh = np.array(inner)[:,0]
 292                             yh = np.array(inner)[:,1]
 293                         shapeh = Polygon([(l[0],l[1]) for l in zip(xh,yh)]).buffer(0)
 294                         shape.difference(shapeh)
 295                     poly = transform(proj,shape)
 296                     polys.append(poly)
 297
 298             idn = prefix + time_data
 299             if idn in perimeters.keys():
 300                 idn = idn + '_2'
 301             # update perimeters dictionary
 302             perimeters.update({idn: Dict({'file': file, 'lon': lons, 'lat': lats,
 303                             'fire': fires, 'conf_fire': np.array(conf_fire*np.ones(lons[fires==9].shape)),
 304                             'lon_fire': lons[fires==9], 'lat_fire': lats[fires==9], 'lon_nofire': lons[fires==5], 'lat_nofire': lats[fires==5],
 305                             'scan_fire': scan*np.ones(lons[fires==9].shape), 'track_fire': track*np.ones(lons[fires==9].shape),
 306                             'conf_nofire': np.array(conf_nofire*np.ones(lons[fires==5].shape)),
 307                             'scan_nofire': scan*np.ones(lons[fires==5].shape), 'track_nofire': track*np.ones(lons[fires==9].shape),
 308                             'time_iso': time_iso, 'time_num': time_num, 'acq_date': acq_date, 'acq_time': acq_time})})
 309             if gen_polys:
 310                 perimeters[idn].update({'polys': polys})
 311     else:
 312         print 'Warning: No KML files in the path specified'
 313         perimeters = []
 314     return perimeters
 315
 316
 317 if __name__ == "__main__":
 318     # Experiment to do
 319     exp = 4
 320     # Plot perimeters as created
 321     plot = True
 322
 323     # Defining options
 324     def pioneer():
 325         bounds = (-115.97282409667969, -115.28449249267578, 43.808258056640625, 44.302913665771484)
 326         time_iso = ('2016-07-18T00:00:00Z', '2016-08-31T00:00:00Z')
 327         igns = None
 328         perims = './pioneer/perim'
 329         return bounds, time_iso, igns, perims
 330     def patch():
 331         bounds = (-113.85068, -111.89413, 39.677563, 41.156837)
 332         time_iso = ('2013-08-10T00:00:00Z', '2013-08-15T00:00:00Z')
 333         igns = ([-112.676039],[40.339372],['2013-08-10T20:00:00Z'])
 334         perims = './patch/perim'
 335         return bounds, time_iso, igns, perims
 336     def saddleridge():
 337         bounds = (-118.60684204101562, -118.35965728759766, 34.226539611816406, 34.43047332763672)
 338         time_iso = ('2019-10-10T00:00:00Z', '2019-10-15T00:00:00Z')
 339         igns = None
 340         perims = './saddleridge/perim'
 341         return bounds, time_iso, igns, perims
 342     def polecreek():
 343         bounds = (-111.93914, -111.311035, 39.75985, 40.239746)
 344         time_iso = ('2018-09-09T00:00:00Z', '2018-09-23T00:00:00Z')
 345         igns = None
 346         perims = './polecreek/perim'
 347         return bounds, time_iso, igns, perims
 348
 349     # Creating the options
 350     options = {1: pioneer, 2: patch, 3: saddleridge, 4: polecreek}
 351
 352     # Defining the option depending on the experiment
 353     bounds, time_iso, igns, perims = options[exp]()
 354
 355     # Processing infrared perimeters
 356     p = process_infrared_perimeters(perims,bounds,time=time_iso,plot=plot,gen_polys=True)
 357     # Processing ignitions if defined
 358     if igns:
 359         p.update(process_ignitions(igns,bounds,time=time_iso))
 360     # Saving results
 361     sl.save(p,'perimeters')
 362