How to aggregate daily NASA MODIS active fires using python ?
Author: Benjamin Marchant
License: CC BY 4.0
In [10]:
from calendar import monthrange
from matplotlib.pyplot import figure
from pylab import matplotlib
from pyhdf.SD import SD, SDC
from pyhdf.HDF import *
from pyhdf.VS import *
import numpy as np
import numpy.ma as ma
import pandas as pd
import glob
import re
import matplotlib.pyplot as plt
import matplotlib as mpl
import cartopy.crs as ccrs
import matplotlib.cm as cm
import warnings
warnings.filterwarnings('ignore')
In [11]:
root = '/Volumes/HD2/Datasets/Research'
In [12]:
year = 2019
month = 8
day = 1
In [13]:
platform = 'AQUA'
instrument = 'MODIS'
product = 'MYD14'
MODIS_Files = glob.glob('{}/NASA/{}/{}/{}/{:04d}/{:04d}_{:02d}_{:02d}/*.hdf'.format(
root,
platform,
instrument,
product,
year,
year,
month,
day))
print( len(MODIS_Files) )
In [21]:
%%time
df = pd.DataFrame()
for idx,MODIS_File in enumerate(MODIS_Files):
df_new = pd.DataFrame()
f = SD(MODIS_File, SDC.READ)
sds_obj = f.select('fire mask') # select sds
fm_data = sds_obj.get() # get sds data
if fm_data[ fm_data > 6 ].shape[0] > 0:
for key in ['FP_latitude',
'FP_longitude',
'FP_ViewZenAng',
'FP_power']:
sds_obj = f.select(key) # select sds
data = sds_obj.get() # get sds data
df_new[key] = data.ravel()
df = pd.concat([df,df_new], axis=0)
f.end()
In [22]:
df
Out[22]:
In [23]:
plt.figure(figsize=(16,9))
proj = ccrs.PlateCarree()
ease_extent = [-180., 180., 90., -90.]
ax = plt.axes(projection=proj)
ax.set_extent(ease_extent, crs=proj)
ax.gridlines(color='gray', linestyle='--')
ax.coastlines()
longs = df['FP_longitude']
lats = df['FP_latitude']
plt.scatter( longs, lats,
color='red', linewidth=2, marker='o', s=2,
transform=ccrs.PlateCarree(),
)
#plt.title('({}) VIIRS AF IBands {}-{}-{} / FP_PersistentAnomalyCategory = {} ({})'.format(platform,year,month,day,PersistentAnomalyCategory,FP_PersistentAnomalyCategory_Legend[PersistentAnomalyCategory]),fontsize=16)
plt.tight_layout()
#plt.savefig('./outputs/{}_AF_Iband_PAC_{}_2022_07_08.png'.format(platform,PersistentAnomalyCategory), dpi=150, facecolor='white')
plt.show()
plt.close()
In [ ]:
In [24]:
resolution = 0.5
In [25]:
def quantitative_to_categorical_converter(x,col1,col2,resolution):
lat_idx = int( (x[col1]+90.0) / resolution )
long_idx = int( (x[col2]+180) / resolution )
return int( 360 / resolution ) * lat_idx + long_idx
In [26]:
%%time
df['Index'] = df.apply(quantitative_to_categorical_converter,
axis=1,
col1 = 'FP_latitude',
col2 = 'FP_longitude',
resolution=resolution)
df[ ['Index', 'FP_latitude', 'FP_longitude'] ]
Out[26]:
In [ ]:
In [27]:
dfg = df[['Index','FP_power']].groupby(['Index']).sum()
dfg
Out[27]:
In [33]:
dfg = dfg.reset_index()
dfg
Out[33]:
In [ ]:
In [28]:
data = np.zeros(( int(180 / resolution) , int(360 / resolution) ))
shape = data.shape
print(data.shape)
print(data.size)
In [29]:
df_map = pd.DataFrame()
xv, yv = np.meshgrid(np.arange( -90,90, resolution), np.arange(-180,180,resolution))
df_map['latitude'] = xv.ravel()
df_map['longitude'] = yv.ravel()
df_map
Out[29]:
In [30]:
%%time
df_map['Index'] = df_map.apply(quantitative_to_categorical_converter,
axis=1,
col1 = 'latitude',
col2 = 'longitude',
resolution=resolution)
In [31]:
df_map
Out[31]:
In [34]:
df_merged = pd.merge(dfg,df_map, on=['Index'], how='right')
df_merged
Out[34]:
In [35]:
df_merged[ df_merged['FP_power'] > 0 ]
Out[35]:
In [ ]:
In [36]:
data = df_merged['FP_power'].to_numpy()
#data = data.reshape(180, 360)
data = data.reshape(shape[1],shape[0])
In [37]:
data = np.nan_to_num(data)
In [ ]:
In [57]:
cmap = cm.get_cmap('hot', 100)
#cmap = cm.get_cmap('jet', 100)
color_list = ['#808080']
#color_list = ['#FFFFFF']
for i in range(cmap.N):
rgba = cmap(i)
# rgb2hex accepts rgb or rgba
#print(matplotlib.colors.rgb2hex(rgba))
color_list.append(matplotlib.colors.rgb2hex(rgba))
cmap = color_list
cmap = mpl.colors.ListedColormap(cmap)
In [70]:
plt.figure(figsize=(16,9))
proj = ccrs.PlateCarree()
ease_extent = [-180., 180., 90., -90.]
ax = plt.axes(projection=proj)
ax.set_extent(ease_extent, crs=proj)
ax.gridlines(color='gray', linestyle='--')
ax.coastlines()
m = ax.imshow(data.T, origin='lower', extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap=cmap,vmin=0, vmax=100.0)
plt.title('AQUA MODIS Fire Radiative Power (2019-08-01)',fontsize=16)
plt.tight_layout()
plt.colorbar(m,ax=ax, fraction=0.02)
plt.savefig("./outputs/Fire_Radiative_Power_Map.png", bbox_inches='tight', dpi=200, facecolor='white')
plt.show()
plt.close()
In [ ]:
In [ ]:
In [ ]: