plot_pixels.py

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