plot_pixels.py

   1 from utils import *
   2 import numpy as np
   3 import matplotlib as mpl
   4 mpl.use('Agg')
   5 import matplotlib.pyplot as plt
   6 from matplotlib.patches import Rectangle
   7 from matplotlib.collections import PatchCollection
   8 from matplotlib.transforms import Affine2D
   9 from mpl_toolkits.mplot3d import Axes3D
  10 from matplotlib.colors import LinearSegmentedColormap
  11 from mpl_toolkits.basemap import Basemap
  12 import saveload as sl
  13 from interpolation import sort_dates
  14 import StringIO
  15 import sys
  16 import os
  17
  18 def create_pixels(lons,lats,widths,heights,alphas,color,label):
  19         """
  20         Plot of pixels using centers (lons,lats), with sizes (widhts,heights), angle alphas, and color color.
  21
  22         :param lons: array of longitude coordinates at the center of the pixels
  23         :param lats: array of latitude coordinates at the center of the pixels
  24         :param widths: array of widhts of the pixels
  25         :param heights: array of heights of the pixels
  26         :param alphas: array of angles of the pixels
  27         :param color: tuple with RGB color values
  28         :return col: collection of rectangles with the pixels
  29
  30         Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  31         Angel Farguell (angel.farguell@gmail.com) 2018-12-28
  32         """
  33
  34         # Create pixels
  35         pixels=[]
  36         for x, y, h, w, a in zip(lons, lats, heights, widths, alphas):
  37                 xpos = x - w/2 # The x position will be half the width from the center
  38                 ypos = y - h/2 # same for the y position, but with height
  39                 rect = Rectangle( (xpos, ypos), w, h, a) # Create a rectangle
  40                 pixels.append(rect) # Add the rectangle patch to our list
  41
  42         # Create a collection from the rectangles
  43         col = PatchCollection(pixels)
  44         # set the alpha for all rectangles
  45         col.set_alpha(0.5)
  46         # Set the colors using the colormap
  47         col.set_facecolor( [color]*len(lons) )
  48         # No lines
  49         col.set_linewidth( 0 )
  50
  51         return col
  52
  53
  54 def basemap_scatter_mercator(g,bounds,bmap,only_fire=False):
  55         """
  56         Scatter plot of fire and ground pixels in a png file using a basemap with mercator projection
  57
  58         :param g: granule dictionary from read_*_files functions in JPSSD.py
  59         :param bounds: array with the coordinates of the bounding box of the case
  60         :param map: basemap mercator map where to plot the center of the pixels
  61         :return png: string with png plot
  62         :return float_bounds: coordinates of the borders of the png file
  63
  64         Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
  65         Angel Farguell (angel.farguell@gmail.com) 2018-04-05
  66         """
  67
  68         size = 50
  69         # Satellite pixels
  70         flon = np.ravel(g.lon)
  71         flat = np.ravel(g.lat)
  72         mask = np.ravel(g.fire)
  73
  74         fil = np.logical_and(np.logical_and(np.logical_and(flon >= bounds[0], flon <= bounds[1]), flat >= bounds[2]), flat <= bounds[3])
  75
  76         if only_fire:
  77                 categories = (np.array(mask[fil] == 7), np.array(mask[fil] == 8),
  78                                         np.array(mask[fil] == 9))
  79                 alphas = (.5,.6,.7)
  80                 labels = ('Fire low','Fire nominal','Fire high')
  81                 colors = [(1,1,0),(1,.65,0),(.5,0,0)]
  82         else:
  83                 categories = (np.array(mask[fil] == 3), np.array(mask[fil] == 5),
  84                                         np.array(mask[fil] == 7), np.array(mask[fil] == 8),
  85                                         np.array(mask[fil] == 9))
  86                 alphas = (.4,.4,.5,.6,.7)
  87                 labels = ('Water','Ground','Fire low','Fire nominal','Fire high')
  88                 colors = [(0,0,.5),(0,.5,0),(1,1,0),(1,.65,0),(.5,0,0)]
  89
  90         lon = []
  91         lat = []
  92         val = []
  93         for i,cat in enumerate(categories):
  94                 lon.append(flon[fil][cat])
  95                 lat.append(flat[fil][cat])
  96                 val.append(np.ones(lon[i].shape)*i)
  97
  98         N = len(categories)
  99         fig = plt.figure(frameon=False,figsize=(12,8))
 100         plt.axis('off')
 101         cmap = LinearSegmentedColormap.from_list('fire_detections', colors, N=N)
 102         for i in range(N):
 103                 bmap.scatter(lon[i],lat[i],size,c=val[i],latlon=True,marker='.',cmap=cmap,vmin=-.5,vmax=N-.5,alpha=alphas[i],linewidths=0)
 104
 105         str_io = StringIO.StringIO()
 106         plt.savefig(str_io,bbox_inches='tight',format='png',pad_inches=0,transparent=True)
 107
 108         #colorbar
 109         orientation = 'vertical'
 110         size_in = 2
 111         cb_label = ''
 112
 113         kwargs = { 'norm': mpl.colors.Normalize(-.5,N-.5),
 114                            'orientation': orientation,
 115                            'spacing': 'proportional',
 116                            'ticks': range(0,N),
 117                            'cmap': cmap}
 118
 119         # build figure according to requested orientation
 120         hsize, wsize = (size_in,size_in*0.5) if orientation == 'vertical' else (size_in*0.5,size_in)
 121         fig = plt.figure(figsize=(wsize,hsize))
 122
 123         # proportions that work with axis title (horizontal not tested)
 124         ax = fig.add_axes([.5,.03,.12,.8]) if orientation=='vertical' else fig.add_axes([0.03,.4,.8,.12])
 125
 126         # construct the colorbar and modify properties
 127         cb = mpl.colorbar.ColorbarBase(ax,**kwargs)
 128         cb.ax.tick_params(length=0)
 129         cb.ax.set_yticklabels(labels)
 130         cb.set_label(cb_label,color='0',fontsize=8,labelpad=-50)
 131
 132         # move ticks to left side
 133         ax.yaxis.set_ticks_position('left')
 134         for tick_lbl in ax.get_yticklabels():
 135                 tick_lbl.set_color('0')
 136                 tick_lbl.set_fontsize(5)
 137
 138         #plt.show()
 139         plt.close()
 140
 141         png = str_io.getvalue()
 142
 143         numpy_bounds = [ (bounds[0],bounds[2]),(bounds[1],bounds[2]),(bounds[1],bounds[3]),(bounds[0],bounds[3]) ]
 144         float_bounds = [ (float(x), float(y)) for x,y in numpy_bounds ]
 145         return png, float_bounds
 146
 147
 148 def center_pixels_plot(g,bounds):
 149         """
 150         Center pixels plot: generates a plot of the center of the pixels.
 151
 152         :param g: granule dictionary from read_*_files functions in JPSSD.py
 153         :param bounds: array with the coordinates of the bounding box of the case
 154         :return: 2D plot of the center of the pixels
 155
 156         Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 157         Angel Farguell (angel.farguell@gmail.com) 2018-12-28
 158         """
 159
 160         fig=plt.figure()
 161         ax=fig.add_subplot(111)
 162
 163         # Size of the center of the pixels
 164         size=80
 165
 166         # Ground pixels
 167         lons=np.array(g.lon_nofire)
 168         lats=np.array(g.lat_nofire)
 169         color=(0,0.59765625,0)
 170         plt.scatter(lons,lats,size,marker='.',color=color,edgecolors='k')
 171
 172         # Fire pixels
 173         lons=np.array(g.lon_fire)
 174         lats=np.array(g.lat_fire)
 175         color=(0.59765625,0,0)
 176         plt.scatter(lons,lats,size,marker='.',color=color,edgecolors='k')
 177
 178         ax.set_xlim(bounds[0],bounds[1])
 179         ax.set_ylim(bounds[2],bounds[3])
 180
 181
 182 def pixels_plot(g,bounds):
 183         """
 184         Regular pixels plot: generates a plot of the pixels using a regular grid by distances.
 185
 186         :param g: granule dictionary from read_*_files functions in JPSSD.py
 187         :param bounds: array with the coordinates of the bounding box of the case
 188         :return: 2D plot of the pixels
 189
 190         Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 191         Angel Farguell (angel.farguell@gmail.com) 2018-12-28
 192         """
 193
 194         lon=g.lon
 195         lat=g.lat
 196         fire=g.fire
 197
 198         lonh0=lon[:,0:-1]
 199         lonh1=lon[:,1:]
 200         lonv0=lon[0:-1,:]
 201         lonv1=lon[1:,:]
 202         dlonh=lonh0-lonh1
 203         dlonh=np.hstack((dlonh,dlonh[:,[-1]]))
 204         dlonv=lonv0-lonv1
 205         dlonv=np.vstack((dlonv,dlonv[-1]))
 206
 207         lath0=lat[:,0:-1]
 208         lath1=lat[:,1:]
 209         latv0=lat[0:-1,:]
 210         latv1=lat[1:,:]
 211         dlath=lath0-lath1
 212         dlath=np.hstack((dlath,dlath[:,[-1]]))
 213         dlatv=latv0-latv1
 214         dlatv=np.vstack((dlatv,dlatv[-1]))
 215
 216         dh=np.sqrt(dlonh**2+dlath**2)
 217         dv=np.sqrt(dlonv**2+dlatv**2)
 218
 219         # Mean between the distances
 220         w=dh
 221         w[:,1:-1]=np.hstack([(dh[:,[k]]+dh[:,[k+1]])/2 for k in range(dh.shape[1]-2)])
 222         h=dv
 223         h[1:-1,:]=np.vstack([(dv[k]+dv[k+1])/2 for k in range(dv.shape[0]-2)])
 224
 225         '''
 226         # Maximum between the distances
 227         w=dh
 228         w[:,1:-1]=np.hstack([np.maximum(dh[:,[k]],dh[:,[k+1]]) for k in range(dh.shape[1]-2)])
 229         h=dv
 230         h[1:-1,:]=np.vstack([np.maximum(dv[k],dv[k+1]) for k in range(dv.shape[0]-2)])
 231         '''
 232
 233         angle=np.arctan(dlath/dlonh)*180/np.pi
 234
 235         plot=False
 236         if plot:
 237                 fig = plt.figure()
 238                 ax = fig.gca(projection='3d')
 239                 surf = ax.plot_surface(lon,lat,w)
 240                 fig = plt.figure()
 241                 ax = fig.gca(projection='3d')
 242                 surf = ax.plot_surface(lon,lat,h)
 243                 fig = plt.figure()
 244                 ax = fig.gca(projection='3d')
 245                 surf = ax.plot_surface(lon,lat,angle)
 246
 247         flon=np.ravel(lon)
 248         flat=np.ravel(lat)
 249         fil=np.logical_and(np.logical_and(np.logical_and(flon >= bounds[0], flon <= bounds[1]), flat >= bounds[2]), flat <= bounds[3])
 250         lon=flon[fil]
 251         lat=flat[fil]
 252
 253         ff=np.ravel(fire)
 254         ffg=ff[fil]
 255         nofi=np.logical_or(ffg == 3, ffg == 5)
 256         fi=np.array(ffg > 6)
 257
 258         width=np.ravel(w)[fil]
 259         height=np.ravel(h)[fil]
 260         alpha=np.ravel(angle)[fil]
 261
 262         fig = plt.figure()
 263         ax = fig.add_subplot(111)
 264
 265         # Ground pixels
 266         lons=lon[nofi]
 267         lats=lat[nofi]
 268         widths=width[nofi]
 269         heights=height[nofi]
 270         alphas=alpha[nofi]
 271         color=(0,0.59765625,0)
 272         col_nofire=create_pixels(lons,lats,widths,heights,alphas,color,'Ground')
 273         ax.add_collection(col_nofire)
 274
 275         # Fire pixels
 276         lons=lon[fi]
 277         lats=lat[fi]
 278         widths=width[fi]
 279         heights=height[fi]
 280         alphas=alpha[fi]
 281         color=(0.59765625,0,0)
 282         col_fire=create_pixels(lons,lats,widths,heights,alphas,color,'Fire')
 283         ax.add_collection(col_fire)
 284
 285         # Set axis
 286         ax.set_xlim(bounds[0],bounds[1])
 287         ax.set_ylim(bounds[2],bounds[3])
 288
 289
 290 def create_kml(kml_data,kml_path):
 291         """
 292         Creation of KML file with png files from basemap_scatter_mercator function
 293
 294         :param kml_data: dictionary with granules to plot information: name, png_file, bounds, time
 295         :param kml_path: path to kml to generate
 296
 297         Developed in Python 2.7.15 :: Anaconda 4.5.10, on MACINTOSH.
 298         Angel Farguell (angel.farguell@gmail.com) 2018-04-05
 299         """
 300
 301         with open(kml_path,'w') as kml:
 302                 kml.write("""<?xml version="1.0" encoding="UTF-8"?>\n""")
 303                 kml.write("""<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">\n""")
 304                 kml.write("""<Document><name>Observed Data</name>""")
 305                 kml.write("""
 306                 <Folder>
 307                         <name>Overlays</name>""")
 308                 for idx,data in enumerate(kml_data):
 309                         name = data['name']
 310                         time = data['time']
 311                         png = data['png_file']
 312                         bounds = data['bounds']
 313                         coord = (str(bounds[0]),str(bounds[2]),
 314                                 str(bounds[1]),str(bounds[2]),
 315                                 str(bounds[1]),str(bounds[3]),
 316                                 str(bounds[0]),str(bounds[3]),
 317                                 str(bounds[0]),str(bounds[2]))
 318                         kml.write("""
 319                         <GroundOverlay>
 320                                 <name>%s</name>
 321                                 <Icon><href>%s</href></Icon>
 322                                 <gx:TimeStamp><when>%s</when></gx:TimeStamp>""" % (
 323                                 name, png, time))
 324                         kml.write("""
 325                                 <gx:LatLonQuad>
 326                                         <coordinates>
 327                                                 %s,%s,0
 328                                                 %s,%s,0
 329                                                 %s,%s,0
 330                                                 %s,%s,0
 331                                                 %s,%s,0
 332                                         </coordinates>
 333                                 </gx:LatLonQuad>""" % coord)
 334                         kml.write("""
 335                         </GroundOverlay>""")
 336                 kml.write("""
 337                 </Folder>""")
 338                 kml.write("""</Document>\n</kml>\n""")
 339
 340         print 'Created file %s' % kml_path
 341
 342
 343 def perror():
 344         print "Error: python %s pixel_plot_type" % sys.argv[0]
 345         print "  - Center pixels in basemap: 0"
 346         print "  - Center pixels plot: 1"
 347         print "  - Regular pixels plot: 2"
 348
 349 if __name__ == "__main__":
 350         if len(sys.argv)!=2:
 351                 perror()
 352                 sys.exit(1)
 353         else:
 354                 if sys.argv[1] not in ['0','1','2']:
 355                         perror()
 356                         print "Invalid value given: %s" % sys.argv[1]
 357                         sys.exit(1)
 358
 359         try:
 360                 # Upload data
 361                 print 'Loading data...'
 362                 data,fxlon,fxlat,time_num=sl.load('data')
 363         except Exception:
 364                 print "Error: No data file in the directory. First run a case using case.py."
 365                 sys.exit(1)
 366
 367         granules=sort_dates(data)
 368         bounds=[fxlon.min(),fxlon.max(),fxlat.min(),fxlat.max()]
 369
 370         print 'Plotting data...'
 371         # Plot pixels
 372         if sys.argv[1] is '0':
 373                 m = Basemap(projection='merc',llcrnrlat=bounds[2], urcrnrlat=bounds[3], llcrnrlon=bounds[0], urcrnrlon=bounds[1])
 374                 kmld = []
 375                 for idx, g in enumerate(granules):
 376                         raster_png_data,corner_coords = basemap_scatter_mercator(g[1],bounds,m)
 377                         bounds = (corner_coords[0][0],corner_coords[1][0],corner_coords[0][1],corner_coords[2][1])
 378                         pngfile = g[0]+'.png'
 379                         timestamp = g[1].acq_date + 'T' + g[1].acq_time[0:2] + ':' + g[1].acq_time[2:4] + 'Z'
 380                         with open(pngfile, 'w') as f:
 381                                 f.write(raster_png_data)
 382                         print '> File %s saved.' % g[0]
 383                         kmld.append(Dict({'name': g[0], 'png_file': pngfile, 'bounds': bounds, 'time': timestamp}))
 384                 create_kml(kmld,'./doc.kml')
 385                 os.system('zip -r granules.kmz doc.kml *.png')
 386                 print 'Created file granules.kmz'
 387                 os.system('rm doc.kml *_A*_*.png')
 388         elif sys.argv[1] is '1':
 389                 for g in granules:
 390                         center_pixels_plot(g[1],bounds)
 391                         plt.title('Granule %s' % g[0])
 392                 plt.show()
 393         else:
 394                 for g in granules:
 395                         pixels_plot(g[1],bounds)
 396                         plt.title('Granule %s' % g[0])
 397                 plt.show()