In this tutorial, we will learn about how to remap satellite data to a uniform grid. This process involves mapping the data to a uniform grid, which makes it simpler to view, manipulate, and store. The tutorial will focus on remapping a VIIRS L2 Fire Mask Granule.

Author: Benjamin Marchant

In [1]:
from pylab import matplotlib
from datetime import date
from os import path
from scipy.spatial.distance import cdist
from scipy.interpolate import griddata
from scipy.interpolate.interpnd import _ndim_coords_from_arrays
from scipy.spatial import cKDTree

import urllib.request
import urllib.request, json 
import pprint
import os
import pandas as pd
import glob
import netCDF4
import random
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import numpy as np
import numpy.ma as ma
import warnings
import cartopy.crs as ccrs
import numpy.ma as ma

from matplotlib.pyplot import figure

warnings.filterwarnings('ignore')

Read VIIRS NASA IBand AF Mask

You can download the file used in this tutorial by clicking on the provided link:

https://drive.google.com/file/d/1hp5vdonRXYNJ2AEgk-RS_rtij0GSSzy6/view?usp=sharing

In [2]:
filename = 'VNP14IMG.A2019215.2106.001.2019216045602.nc'

file_name = "/Volumes/HD2/Datasets/Research//NASA/Suomi-NPP/VIIRS/VNP14IMG/2019/2019_08_03/VNP14IMG.A2019215.2106.001.2019216045602.nc"
In [3]:
time_stamp = '{}.{}'.format(filename.split('.')[1],filename.split('.')[2])

time_stamp
Out[3]:
'A2019215.2106'
In [4]:
f = netCDF4.Dataset(file_name)
In [5]:
def plot_fire_mask(data,filename,show_plot,dpi,nb_ticks):
    
    cmap = cm.get_cmap('RdBu_r', 11)    # PiYG

    color_list = ['#808080']

    for i in range(cmap.N):
        rgba = cmap(i)
        color_list.append(matplotlib.colors.rgb2hex(rgba))
    
    fig = figure(num=None, figsize=(12, 10), dpi=dpi, facecolor='w', edgecolor='k')

    ax = fig.add_subplot(111)

    cmap = color_list
    cmap = mpl.colors.ListedColormap(cmap)

    bounds = [i for i in range(11)]

    norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

    img = plt.imshow(data, cmap=cmap, norm=norm, interpolation='none')

    cbar_bounds = bounds
    cbar_ticks =  [(cbar_bounds[i+1]-cbar_bounds[i])/2.0+cbar_bounds[i] for i in range( len(cbar_bounds) - 1 )] #  [0.5, 1.5, 155.0, 254.0, 255.0]  
    cbar_labels = [i for i in range(10)]

    cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
    cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

    plt.title('FRP MASK')

    plt.grid(c='black',ls='--')

    plt.xticks(np.linspace(0, data.shape[1], nb_ticks))
    plt.yticks(np.linspace(0, data.shape[0], nb_ticks))

    ax.spines['bottom'].set_color('black')
    ax.spines['top'].set_color('black')
    ax.spines['left'].set_color('black')
    ax.spines['right'].set_color('black')

    #plt.savefig('Fire_mask.png', dpi=200, bbox_inches='tight')
    
    if show_plot:
        plt.show()
    
    plt.close()
    
    return None

masked_data = f.variables['fire mask']
data =  ma.getdata(masked_data)

filename = 'test'

plot_fire_mask(data,filename,show_plot=True,dpi=200,nb_ticks=10)
In [6]:
frp_mask_data =  ma.getdata(masked_data)
In [7]:
agg_y1 = 640
agg_y2 = 368 
agg_y3 = 592
In [8]:
idx_bounds = [0,
              agg_y1*2,
              (agg_y1+agg_y2)*2,
              (agg_y1+agg_y2+2*agg_y3)*2,
              (agg_y1+agg_y2+2*agg_y3+agg_y2)*2,
              (agg_y1+agg_y2+2*agg_y3+agg_y2+agg_y1)*2]

idx_bounds
Out[8]:
[0, 1280, 2016, 4384, 5120, 6400]
In [9]:
data_img_01 = data[:,0:agg_y1*2]
data_img_02 = data[:,agg_y1*2:(agg_y1+agg_y2)*2]
data_img_03 = data[:,(agg_y1+agg_y2)*2:(agg_y1+agg_y2+2*agg_y3)*2]
data_img_04 = data[:,(agg_y1+agg_y2+2*agg_y3)*2:(agg_y1+agg_y2+2*agg_y3+agg_y2)*2]
data_img_05 = data[:,(agg_y1+agg_y2+2*agg_y3+agg_y2)*2:]
In [10]:
imgs_list = [data_img_01,data_img_02,data_img_03,data_img_04,data_img_05]
In [11]:
plot_fire_mask(data_img_01,filename,show_plot=True,dpi=200,nb_ticks=5)
In [12]:
plot_fire_mask(data_img_02,filename,show_plot=True,dpi=200,nb_ticks=3)
In [13]:
plot_fire_mask(data_img_03,filename,show_plot=True,dpi=200,nb_ticks=5)
In [14]:
plot_fire_mask(data_img_04,filename,show_plot=True,dpi=200,nb_ticks=3)
In [15]:
plot_fire_mask(data_img_05,filename,show_plot=True,dpi=200,nb_ticks=5)
In [16]:
f.close()

Get Geolocation Data

To download the VNP03IMG file, click on the provided link. The next 5 steps demonstrate how to obtain the match VNP03IMG file for a specific day.

https://drive.google.com/file/d/1cgZKa5LtP2Ba_OprFPrAEr1x7PnoCYxH/view?usp=sharing

In [17]:
root = '/Volumes/HD2/Datasets/Research'
In [18]:
VNP03IMG_Files = glob.glob('{}/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/*.nc'.format(root))

len(VNP03IMG_Files)
Out[18]:
240
In [19]:
VNP03IMG_Files[:10]
Out[19]:
['/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0000.001.2019215063252.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0006.001.2019215063432.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0012.001.2019215063428.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0018.001.2019215063430.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0024.001.2019215063437.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0030.001.2019215063433.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0036.001.2019215063434.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0042.001.2019215063434.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0048.001.2019215063434.nc',
 '/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.0054.001.2019215063432.nc']
In [20]:
[i for i in VNP03IMG_Files if time_stamp in i]
Out[20]:
['/Volumes/HD2/Datasets/Research/NASA/Suomi-NPP/VIIRS/VNP03IMG/2019/2019_08_03/VNP03IMG.A2019215.2106.001.2019216024358.nc']
In [21]:
VNP03IMG_filename = [i for i in VNP03IMG_Files if time_stamp in i][0]
In [22]:
f = netCDF4.Dataset(VNP03IMG_filename)
In [23]:
geo_data_gp = f.groups['geolocation_data']
In [24]:
df = pd.DataFrame()

for variable in ['latitude','longitude']:
    masked_data =  geo_data_gp.variables[variable]
    data =  ma.getdata(masked_data)
    df[variable] = data.ravel()

granule_min_long = df['longitude'].min()
granule_max_long = df['longitude'].max()

granule_min_lat = df['latitude'].min()
granule_max_lat = df['latitude'].max()

print(granule_min_long,granule_max_long)
print(granule_min_lat,granule_max_lat)

-156.75677490234375 -103.08659362792969
37.07435989379883 63.50246810913086

Remaping (Approach 1)

In [25]:
df = pd.DataFrame()

for variable in ['latitude','longitude']:
    masked_data =  geo_data_gp.variables[variable]
    data =  ma.getdata(masked_data)
    df[variable] = data.ravel()

sds_shape = frp_mask_data.shape
    
df['fire mask'] = frp_mask_data.ravel()
In [26]:
df
Out[26]:
latitude longitude fire mask
0 42.438839 -104.463242 1
1 42.438988 -104.473495 1
2 42.439144 -104.484024 1
3 42.439297 -104.494576 1
4 42.439445 -104.504715 1
... ... ... ...
41574395 55.158630 -156.715332 1
41574396 55.153946 -156.725677 1
41574397 55.149254 -156.736023 1
41574398 55.144554 -156.746399 1
41574399 55.139851 -156.756775 1

41574400 rows × 3 columns

In [27]:
df_no_bowtie = df[ df['fire mask'] > 1 ]

df_no_bowtie
Out[27]:
latitude longitude fire mask
2016 41.965446 -116.653458 4
2017 41.964771 -116.660751 5
2018 41.964092 -116.668030 5
2019 41.963448 -116.674904 5
2020 41.962830 -116.681458 4
... ... ... ...
41572379 60.060963 -142.336624 5
41572380 60.058632 -142.346039 3
41572381 60.056309 -142.355453 3
41572382 60.053982 -142.364883 3
41572383 60.051651 -142.374313 3

36221696 rows × 3 columns

In [28]:
df_no_bowtie[['longitude','latitude']].to_numpy()
Out[28]:
array([[-116.65346 ,   41.965446],
       [-116.66075 ,   41.96477 ],
       [-116.66803 ,   41.964092],
       ...,
       [-142.35545 ,   60.05631 ],
       [-142.36488 ,   60.05398 ],
       [-142.37431 ,   60.05165 ]], dtype=float32)
In [29]:
xy1_igrid = df_no_bowtie[['longitude','latitude']].to_numpy()

xy1_igrid.shape
Out[29]:
(36221696, 2)
In [30]:
df_no_bowtie[['fire mask']].to_numpy()
Out[30]:
array([[4],
       [5],
       [5],
       ...,
       [3],
       [3],
       [3]], dtype=uint8)
In [31]:
z_igrid_01 = df_no_bowtie[['fire mask']].to_numpy()

z_igrid_01.shape
Out[31]:
(36221696, 1)
In [32]:
z_igrid_01 = z_igrid_01.astype('float64')

z_igrid_01
Out[32]:
array([[4.],
       [5.],
       [5.],
       ...,
       [3.],
       [3.],
       [3.]])
In [33]:
%%time

xi, yi = np.mgrid[granule_min_long:granule_max_long:1000j, granule_min_lat:granule_max_lat:1000j]

z_01 = griddata(xy1_igrid, z_igrid_01, (xi, yi), method='nearest')

#----------#

THRESHOLD = 0.05

tree = cKDTree(xy1_igrid)
arr_x = _ndim_coords_from_arrays((xi, yi))
dists, indexes = tree.query(arr_x)

z_01[dists > THRESHOLD] = np.nan

#----------#

whereAreNaNs = np.isnan(z_01);
z_01[whereAreNaNs] = 0.;

CPU times: user 43.8 s, sys: 931 ms, total: 44.7 s
Wall time: 44.8 s
In [34]:
min_long = granule_min_long
max_long = granule_max_long

min_lat = granule_min_lat
max_lat = granule_max_lat
In [35]:
plt.figure(figsize=(12,12),dpi=200, facecolor='w')

proj = ccrs.PlateCarree()

offset = 0.0

ease_extent = [min_long-offset, 
               max_long+offset, 
               min_lat-offset, 
               max_lat+offset]

ax = plt.axes(projection=proj)

ax.set_extent(ease_extent, crs=proj) 


swe_extent = [min_long, max_long, min_lat, max_lat]


cmap = cm.get_cmap('RdBu_r', 11)    # PiYG

color_list = ['#808080']

for i in range(cmap.N):
    rgba = cmap(i)
    color_list.append(matplotlib.colors.rgb2hex(rgba))


cmap = color_list
cmap = mpl.colors.ListedColormap(cmap)


bounds = [i for i in range(11)]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)


im = ax.imshow(np.rot90(np.fliplr(z_01)), extent=swe_extent, transform=proj, origin='lower', aspect=1., cmap=cmap, norm=norm)

cbar_bounds = bounds
cbar_ticks =  [(cbar_bounds[i+1]-cbar_bounds[i])/2.0+cbar_bounds[i] for i in range( len(cbar_bounds) - 1 )] #  [0.5, 1.5, 155.0, 254.0, 255.0]  
cbar_labels = [i for i in range(10)]

cbar = plt.colorbar(im, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks,fraction=0.020)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

g1 = ax.gridlines(draw_labels=True, color='gray', linestyle='--')

ax.coastlines()

#plt.title('FIREX-AQ MASTER \n Fire Radiative Power (MW) \n {}-{}-{} {}:{} \n FRP TOT {:.2f}MW'.format(year,month,day,hour,minute,FRP_tot),fontsize=14)

g1.top_labels = False
g1.right_labels = False

plt.tight_layout()

plt.title('')

plt.savefig('NASA_VIIRS_AF_Remapping.png', dpi=1000, bbox_inches='tight')

plt.show()

plt.close()

Remaping (Approach 2)

I am sharing my initial attempt at VIIRS remapping, which involves going through various sections of a VIIRS granule. While this approach is not as efficient as the one explained above, I have included it as an additional example.

In [36]:
df = pd.DataFrame()

for variable in ['latitude','longitude']:
    masked_data =  geo_data_gp.variables[variable]
    data =  ma.getdata(masked_data)
    data = data[:,(agg_y1+agg_y2)*2:(agg_y1+agg_y2+2*agg_y3)*2]
    df[variable] = data.ravel()

sds_shape = data_img_03.shape
    
df['fire mask'] = data_img_03.ravel()
    
df
Out[36]:
latitude longitude fire mask
0 41.965446 -116.653458 4
1 41.964771 -116.660751 5
2 41.964092 -116.668030 5
3 41.963448 -116.674904 5
4 41.962830 -116.681458 4
... ... ... ...
15382523 60.060963 -142.336624 5
15382524 60.058632 -142.346039 3
15382525 60.056309 -142.355453 3
15382526 60.053982 -142.364883 3
15382527 60.051651 -142.374313 3

15382528 rows × 3 columns

In [37]:
#min_long = df['longitude'].min()
#max_long = df['longitude'].max()

#min_lat = df['latitude'].min()
#max_lat = df['latitude'].max()

min_long = granule_min_long
max_long = granule_max_long

min_lat = granule_min_lat
max_lat = granule_max_lat

#print(min_long,max_long)
#print(min_lat,max_lat)
In [38]:
proj = ccrs.PlateCarree()

lat_long_grid = proj.transform_points(                 
                    x = df['longitude'].to_numpy().reshape(sds_shape),
                    y = df['latitude'].to_numpy().reshape(sds_shape),
                    src_crs = proj)

x_igrid = lat_long_grid[:,:,0] ## long
y_igrid = lat_long_grid[:,:,1] ## lat

z_igrid_01 = np.zeros(sds_shape)

z_igrid_01[:,:] = df['fire mask'].to_numpy().reshape(sds_shape)

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()

z_igrid_01 = z_igrid_01.ravel()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T

print( xy1_igrid.shape )
print( z_igrid_01.shape )

(15382528, 2)
(15382528,)
In [39]:
xy1_igrid
Out[39]:
array([[-116.65345764,   41.96544647],
       [-116.66075134,   41.96477127],
       [-116.66802979,   41.96409225],
       ...,
       [-142.35545349,   60.05630875],
       [-142.36488342,   60.05398178],
       [-142.37431335,   60.051651  ]])
In [40]:
z_igrid_01
Out[40]:
array([4., 5., 5., ..., 3., 3., 3.])
In [41]:
proj = ccrs.PlateCarree()

lat_long_grid = proj.transform_points(                 
                    x = df['longitude'].to_numpy().reshape(sds_shape),
                    y = df['latitude'].to_numpy().reshape(sds_shape),
                    src_crs = proj)

x_igrid = lat_long_grid[:,:,0] ## long
y_igrid = lat_long_grid[:,:,1] ## lat

#print(x_igrid.shape)

#----------#

z_igrid_01 = np.zeros(sds_shape)

z_igrid_01[:,:] = df['fire mask'].to_numpy().reshape(sds_shape)

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()

z_igrid_01 = z_igrid_01.ravel()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T
xi, yi = np.mgrid[min_long:max_long:1000j, min_lat:max_lat:1000j]

z_01 = griddata(xy1_igrid, z_igrid_01, (xi, yi), method='nearest')

#----------#

THRESHOLD = 0.05

tree = cKDTree(xy1_igrid)
arr_x = _ndim_coords_from_arrays((xi, yi))
dists, indexes = tree.query(arr_x)

z_01[dists > THRESHOLD] = np.nan

#----------#

whereAreNaNs = np.isnan(z_01);
z_01[whereAreNaNs] = 0.;

plt.figure(figsize=(12,12),dpi=200)

proj = ccrs.PlateCarree()

offset = 0.0

ease_extent = [min_long-offset, 
               max_long+offset, 
               min_lat-offset, 
               max_lat+offset]

ax = plt.axes(projection=proj)

ax.set_extent(ease_extent, crs=proj) 


swe_extent = [min_long, max_long, min_lat, max_lat]


cmap = cm.get_cmap('RdBu_r', 11)    # PiYG

color_list = ['#808080']

for i in range(cmap.N):
    rgba = cmap(i)
    color_list.append(matplotlib.colors.rgb2hex(rgba))


cmap = color_list
cmap = mpl.colors.ListedColormap(cmap)


bounds = [i for i in range(11)]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)


im = ax.imshow(np.rot90(np.fliplr(z_01)), extent=swe_extent, transform=proj, origin='lower', aspect=1., cmap=cmap, norm=norm)

cbar_bounds = bounds
cbar_ticks =  [(cbar_bounds[i+1]-cbar_bounds[i])/2.0+cbar_bounds[i] for i in range( len(cbar_bounds) - 1 )] #  [0.5, 1.5, 155.0, 254.0, 255.0]  
cbar_labels = [i for i in range(10)]

cbar = plt.colorbar(im, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks,fraction=0.026)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

g1 = ax.gridlines(draw_labels=True, color='gray', linestyle='--')

ax.coastlines()

#plt.title('FIREX-AQ MASTER \n Fire Radiative Power (MW) \n {}-{}-{} {}:{} \n FRP TOT {:.2f}MW'.format(year,month,day,hour,minute,FRP_tot),fontsize=14)

g1.top_labels = False
g1.right_labels = False

plt.tight_layout()

plt.title('')

#plt.savefig('./outputs/FIREX_AQ_MASTER_FRP_Regridding_{}{}{}_{}{}.png'.format(year,month,day,hour,minute), dpi=100, bbox_inches='tight')

plt.show()

plt.close()
In [42]:
data_img_03_z = z_01.copy()
In [ ]:
In [43]:
df = pd.DataFrame()

for variable in ['latitude','longitude']:
    masked_data =  geo_data_gp.variables[variable]
    data =  ma.getdata(masked_data)
    data = data[:,0:agg_y1*2]
    df[variable] = data.ravel()

sds_shape = data_img_01.shape
    
df['fire mask'] = data_img_01.ravel()
    
df
Out[43]:
latitude longitude fire mask
0 42.438839 -104.463242 1
1 42.438988 -104.473495 1
2 42.439144 -104.484024 1
3 42.439297 -104.494576 1
4 42.439445 -104.504715 1
... ... ... ...
8314875 63.444656 -116.048149 1
8314876 63.444458 -116.054916 1
8314877 63.444275 -116.061554 1
8314878 63.444092 -116.068161 1
8314879 63.443905 -116.074722 1

8314880 rows × 3 columns

In [44]:
df['fire mask'].value_counts()
Out[44]:
5    3219105
4    2529832
1    2078720
3     487190
8         19
9         11
7          3
Name: fire mask, dtype: int64
In [45]:
df = df[ df['fire mask'] != 1 ]

df.shape
Out[45]:
(6236160, 3)
In [46]:
print( data_img_01.shape )
print( data_img_01[ data_img_01 == 1 ].shape ) 

(6496, 1280)
(2078720,)
In [47]:
data_img_01[ data_img_01 == 1 ]
Out[47]:
masked_array(data=[1, 1, 1, ..., 1, 1, 1],
             mask=False,
       fill_value=999999,
            dtype=uint8)
In [48]:
data_img_01[ data_img_01[:,-1] != 1 ].shape
Out[48]:
(4872, 1280)
In [49]:
data_img_01[ data_img_01[:,-1] != 1 ]
Out[49]:
masked_array(
  data=[[5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        [5, 5, 5, ..., 5, 5, 5],
        ...,
        [4, 4, 4, ..., 4, 4, 4],
        [4, 4, 4, ..., 4, 4, 4],
        [4, 4, 4, ..., 4, 4, 4]],
  mask=False,
  fill_value=999999,
  dtype=uint8)
In [50]:
data_img_01_no_bowtie = ma.getdata(data_img_01[ data_img_01[:,-1] != 1 ])

data_img_01_no_bowtie
Out[50]:
array([[5, 5, 5, ..., 5, 5, 5],
       [5, 5, 5, ..., 5, 5, 5],
       [5, 5, 5, ..., 5, 5, 5],
       ...,
       [4, 4, 4, ..., 4, 4, 4],
       [4, 4, 4, ..., 4, 4, 4],
       [4, 4, 4, ..., 4, 4, 4]], dtype=uint8)
In [51]:
sds_shape = data_img_01_no_bowtie.shape

sds_shape
Out[51]:
(4872, 1280)
In [52]:
4872 * 1280
Out[52]:
6236160
In [53]:
plot_fire_mask(data_img_01_no_bowtie,filename,show_plot=True,dpi=200,nb_ticks=5)
In [54]:
#min_long = df['longitude'].min()
#max_long = df['longitude'].max()

#min_lat = df['latitude'].min()
#max_lat = df['latitude'].max()

#print(min_long,max_long)
#print(min_lat,max_lat)
In [55]:
proj = ccrs.PlateCarree()

lat_long_grid = proj.transform_points(                 
                    x = df['longitude'].to_numpy().reshape(sds_shape),
                    y = df['latitude'].to_numpy().reshape(sds_shape),
                    src_crs = proj)

x_igrid = lat_long_grid[:,:,0] ## long
y_igrid = lat_long_grid[:,:,1] ## lat

#print(x_igrid.shape)

#----------#

z_igrid_01 = np.zeros(sds_shape)

z_igrid_01[:,:] = df['fire mask'].to_numpy().reshape(sds_shape)

x1_igrid = x_igrid.ravel()
y1_igrid = y_igrid.ravel()

z_igrid_01 = z_igrid_01.ravel()

xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T
xi, yi = np.mgrid[granule_min_long:granule_max_long:1000j, granule_min_lat:granule_max_lat:1000j]

z_01 = griddata(xy1_igrid, z_igrid_01, (xi, yi), method='nearest')

#----------#

THRESHOLD = 0.05

tree = cKDTree(xy1_igrid)
arr_x = _ndim_coords_from_arrays((xi, yi))
dists, indexes = tree.query(arr_x)

z_01[dists > THRESHOLD] = np.nan

#----------#

whereAreNaNs = np.isnan(z_01);
z_01[whereAreNaNs] = 0.;

plt.figure(figsize=(12,12),dpi=200)

proj = ccrs.PlateCarree()

offset = 0.0

ease_extent = [min_long-offset, 
               max_long+offset, 
               min_lat-offset, 
               max_lat+offset]

ax = plt.axes(projection=proj)

ax.set_extent(ease_extent, crs=proj) 


swe_extent = [min_long, max_long, min_lat, max_lat]


cmap = cm.get_cmap('RdBu_r', 11)    # PiYG

color_list = ['#808080']

for i in range(cmap.N):
    rgba = cmap(i)
    color_list.append(matplotlib.colors.rgb2hex(rgba))


cmap = color_list
cmap = mpl.colors.ListedColormap(cmap)


bounds = [i for i in range(11)]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)


im = ax.imshow(np.rot90(np.fliplr(z_01)), extent=swe_extent, transform=proj, origin='lower', aspect=1., cmap=cmap, norm=norm)

cbar_bounds = bounds
cbar_ticks =  [(cbar_bounds[i+1]-cbar_bounds[i])/2.0+cbar_bounds[i] for i in range( len(cbar_bounds) - 1 )] #  [0.5, 1.5, 155.0, 254.0, 255.0]  
cbar_labels = [i for i in range(10)]

cbar = plt.colorbar(im, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks,fraction=0.020)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

g1 = ax.gridlines(draw_labels=True, color='gray', linestyle='--')

ax.coastlines()

#plt.title('FIREX-AQ MASTER \n Fire Radiative Power (MW) \n {}-{}-{} {}:{} \n FRP TOT {:.2f}MW'.format(year,month,day,hour,minute,FRP_tot),fontsize=14)

g1.top_labels = False
g1.right_labels = False

plt.tight_layout()

plt.title('')

#plt.savefig('./outputs/FIREX_AQ_MASTER_FRP_Regridding_{}{}{}_{}{}.png'.format(year,month,day,hour,minute), dpi=100, bbox_inches='tight')

plt.show()

plt.close()
In [56]:
data_img_01_z = z_01.copy()
In [ ]:
In [57]:
plt.figure(figsize=(12,12),dpi=200)

proj = ccrs.PlateCarree()

offset = 0.0

ease_extent = [min_long-offset, 
               max_long+offset, 
               min_lat-offset, 
               max_lat+offset]

ax = plt.axes(projection=proj)

ax.set_extent(ease_extent, crs=proj) 


swe_extent = [min_long, max_long, min_lat, max_lat]


cmap = cm.get_cmap('RdBu_r', 11)    # PiYG

color_list = ['#808080']

for i in range(cmap.N):
    rgba = cmap(i)
    color_list.append(matplotlib.colors.rgb2hex(rgba))


cmap = color_list
cmap = mpl.colors.ListedColormap(cmap)


bounds = [i for i in range(11)]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)


im = ax.imshow(np.rot90(np.fliplr(data_img_01_z+data_img_03_z)), extent=swe_extent, transform=proj, origin='lower', aspect=1., cmap=cmap, norm=norm)

cbar_bounds = bounds
cbar_ticks =  [(cbar_bounds[i+1]-cbar_bounds[i])/2.0+cbar_bounds[i] for i in range( len(cbar_bounds) - 1 )] #  [0.5, 1.5, 155.0, 254.0, 255.0]  
cbar_labels = [i for i in range(10)]

cbar = plt.colorbar(im, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks,fraction=0.020)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

g1 = ax.gridlines(draw_labels=True, color='gray', linestyle='--')

ax.coastlines()

#plt.title('FIREX-AQ MASTER \n Fire Radiative Power (MW) \n {}-{}-{} {}:{} \n FRP TOT {:.2f}MW'.format(year,month,day,hour,minute,FRP_tot),fontsize=14)

g1.top_labels = False
g1.right_labels = False

plt.tight_layout()

plt.title('')

#plt.savefig('./outputs/FIREX_AQ_MASTER_FRP_Regridding_{}{}{}_{}{}.png'.format(year,month,day,hour,minute), dpi=100, bbox_inches='tight')

plt.show()

plt.close()
In [ ]:
In [58]:
res_z = []

for i in range(5):
    
    data_img = imgs_list[i]

    df = pd.DataFrame()

    for variable in ['latitude','longitude']:
        masked_data =  geo_data_gp.variables[variable]
        data =  ma.getdata(masked_data)
        data = data[:,idx_bounds[i]:idx_bounds[i+1]]
        df[variable] = data.ravel()

    sds_shape = data_img.shape

    df['fire mask'] = data_img.ravel()

    df = df[ df['fire mask'] != 1 ]

    data_img_no_bowtie = ma.getdata(data_img[ data_img[:,-1] != 1 ])

    sds_shape = data_img_no_bowtie.shape
    
    proj = ccrs.PlateCarree()

    lat_long_grid = proj.transform_points(                 
                        x = df['longitude'].to_numpy().reshape(sds_shape),
                        y = df['latitude'].to_numpy().reshape(sds_shape),
                        src_crs = proj)

    x_igrid = lat_long_grid[:,:,0] ## long
    y_igrid = lat_long_grid[:,:,1] ## lat

    #print(x_igrid.shape)

    #----------#

    z_igrid_01 = np.zeros(sds_shape)

    z_igrid_01[:,:] = df['fire mask'].to_numpy().reshape(sds_shape)

    x1_igrid = x_igrid.ravel()
    y1_igrid = y_igrid.ravel()

    z_igrid_01 = z_igrid_01.ravel()

    xy1_igrid = np.vstack((x1_igrid, y1_igrid)).T
    xi, yi = np.mgrid[granule_min_long:granule_max_long:1000j, granule_min_lat:granule_max_lat:1000j]

    z_01 = griddata(xy1_igrid, z_igrid_01, (xi, yi), method='nearest')

    #----------#

    THRESHOLD = 0.05

    tree = cKDTree(xy1_igrid)
    arr_x = _ndim_coords_from_arrays((xi, yi))
    dists, indexes = tree.query(arr_x)

    z_01[dists > THRESHOLD] = np.nan

    #----------#

    whereAreNaNs = np.isnan(z_01);
    z_01[whereAreNaNs] = 0.;

    plt.figure(figsize=(12,12),dpi=200)

    proj = ccrs.PlateCarree()

    offset = 0.0

    ease_extent = [min_long-offset, 
                   max_long+offset, 
                   min_lat-offset, 
                   max_lat+offset]

    ax = plt.axes(projection=proj)

    ax.set_extent(ease_extent, crs=proj) 


    swe_extent = [min_long, max_long, min_lat, max_lat]


    cmap = cm.get_cmap('RdBu_r', 11)    # PiYG

    color_list = ['#808080']

    for i in range(cmap.N):
        rgba = cmap(i)
        color_list.append(matplotlib.colors.rgb2hex(rgba))


    cmap = color_list
    cmap = mpl.colors.ListedColormap(cmap)


    bounds = [i for i in range(11)]

    norm = mpl.colors.BoundaryNorm(bounds, cmap.N)


    im = ax.imshow(np.rot90(np.fliplr(z_01)), extent=swe_extent, transform=proj, origin='lower', aspect=1., cmap=cmap, norm=norm)

    cbar_bounds = bounds
    cbar_ticks =  [(cbar_bounds[i+1]-cbar_bounds[i])/2.0+cbar_bounds[i] for i in range( len(cbar_bounds) - 1 )] #  [0.5, 1.5, 155.0, 254.0, 255.0]  
    cbar_labels = [i for i in range(10)]

    cbar = plt.colorbar(im, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks,fraction=0.020)
    cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

    g1 = ax.gridlines(draw_labels=True, color='gray', linestyle='--')

    ax.coastlines()

    #plt.title('FIREX-AQ MASTER \n Fire Radiative Power (MW) \n {}-{}-{} {}:{} \n FRP TOT {:.2f}MW'.format(year,month,day,hour,minute,FRP_tot),fontsize=14)

    g1.top_labels = False
    g1.right_labels = False

    plt.tight_layout()

    plt.title('')

    #plt.savefig('./outputs/FIREX_AQ_MASTER_FRP_Regridding_{}{}{}_{}{}.png'.format(year,month,day,hour,minute), dpi=100, bbox_inches='tight')

    plt.show()

    plt.close()
    
    print('Done')
    
    data_img_z = z_01.copy()
    
    res_z.append(data_img_z)

Done

Done

Done