<jupyter_start><jupyter_text>Radar Nowcast ChallengeThis notebook will describe the SEVIR radar nowcasting challenge.  Nowcasts (https://en.wikipedia.org/wiki/Nowcasting_(meteorology)) are short-term forecast of weather variables typically measured by weather radar or satellite.   Nowcasts are different from traditional weather forecasts (like those you see on the news, or on your phones weather app) in that they are based (mostly) on statistical extrapolations of recent data, rather than full physics-based numerical weather prediction (NWP) models.   The advantages of using statistical extrapolation techniques is that they generally run much faster than NWP models (seconds compared to hours), and because of this, nowcasts are able to  leverage the most recently observed data and use it for forecasting.  Because of this extremely low latency, nowcast out-perform traditional weather models in terms of accuracy and precision for short look-aheads (typically 1-2 hours at most). Below is an example that a nowcast generated using the baseline model created in this notebook.  The first half of the movie is observed radar VIL.  The second half of the movie is the forecast generated by a model trained in this notebook.   The challenge is to create the best possible nowcast algorithm given the data in SEVIR.![imgs/mrms_animation.gif](imgs/mrms_animation.gif) Contents:* [Python environment](env)* [Problem Statement](problem)* [Training & testing dataset](datasets)    * [Visualizing training samples](vistraining)* [Model Training](training)* [Visualize Result](visualize)* [Forecast Scoring](scoring)* [Applying to full-sized images](fillimgs)  Python Environment Running this notebook requires a python environment that includes the modules described here.  You can setup the environment with anaconda via```conda create --name sevir_challengesconda activate sevir_challengesconda install --file requirements.txt```  Problem Statement The goal of this challenge is to predict 12 future frames of radar images in 5 minute steps (`vil` in this case) given the previous 13 images of `vil` (also in 5 minute steps) in addition to other modalities in SEVIR. ![imgs/challenge_overview_diagram.png](imgs/challenge_overview_diagram.png)  Training and testing datasets   MIT Supercloud If you have access to MIT supercloud, the data for this challenge can be accessed directly:  AWS Download Sample training and testing datasets are provided on the SEVIR S3 bucket.  These can be obtained by running the following cell.  **This downloads and decompresses files that will take up roughly 40 GB so make sure you have sufficient disk space**<jupyter_code>data_path="/scratch/li.baol/SEVIR"
# Target locations of sample training & testing data
DEST_TRAIN_FILE= os.path.join(data_path,'data/processed/nowcast_training_000.h5')
DEST_TRAIN_META=os.path.join(data_path, 'data/processed/nowcast_training_000_META.csv')
DEST_TEST_FILE= os.path.join(data_path, 'data/processed/nowcast_testing_000.h5')
DEST_TEST_META= os.path.join(data_path, 'data/processed/nowcast_testing_000_META.csv')
# THIS DOWNLOADS APPROXIMATELY 40 GB DATASETS (AFTER DECOMPRESSION)
import boto3
from botocore.handlers import disable_signing
import tarfile
resource = boto3.resource('s3')
resource.meta.client.meta.events.register('choose-signer.s3.*', disable_signing)
bucket=resource.Bucket('sevir')

print('Dowloading sample training data')
if not os.path.exists(DEST_TRAIN_FILE):
    bucket.download_file('data/processed/nowcast_training_000.h5.tar.gz',DEST_TRAIN_FILE+'.tar.gz')
    bucket.download_file('data/processed/nowcast_training_000_META.csv',DEST_TRAIN_META)
    with tarfile.open(DEST_TRAIN_FILE+'.tar.gz') as tfile:
        tfile.extract('data/processed/nowcast_training_000.h5',data_path)
else:
    print('Train file %s already exists' % DEST_TRAIN_FILE)
print('Dowloading sample testing data')
if not os.path.exists(DEST_TEST_FILE):
    bucket.download_file('data/processed/nowcast_testing_000.h5.tar.gz',DEST_TEST_FILE+'.tar.gz')
    bucket.download_file('data/processed/nowcast_testing_000_META.csv',DEST_TEST_META)
    with tarfile.open(DEST_TEST_FILE+'.tar.gz') as tfile:
        tfile.extract('data/processed/nowcast_testing_000.h5',data_path)
else:
    print('Test file %s already exists' % DEST_TEST_FILE)<jupyter_output>Dowloading sample training data
Train file /scratch/li.baol/SEVIR/data/processed/nowcast_training_000.h5 already exists
Dowloading sample testing data
Test file /scratch/li.baol/SEVIR/data/processed/nowcast_testing_000.h5 already exists<jupyter_text>See the contents of the training data file using `h5ls`:<jupyter_code>#!h5ls /scratch/li.baol/SEVIR/data/processed/nowcast_training_000.h5<jupyter_output><empty_output><jupyter_text>This file contains four modalities you can use as input to a nowcast model (keys start with `IN_` -- `vil`,`ir_ir107`,`ir069`, or `lght`).  The challenge is to use these inputs to predict the output vil `OUT_vil`.     If you have downloaded the full SEVIR dataset, you can extract even larger datasets of this format, or even modify the dataset for your puprose.  We assume SEVIR is located in the directory `$SEVIR_ROOT`.  To make your own dataset, run the script `make_dataset.py`:```python make_dataset.py --input_types vil --output_types vil --sevir_data $SEVIR_ROOT/data/ --sevir_catalog $SEVIR_ROOT/CATALOG.csv --output_location ../data/processed/``` The rest of this notebook assumes you are using the sample training and test files located in the directory `../data/processed/` relative to this notebook.  Adjust as needed.  Exploring the training dataset First load some samples from the training data and assocaited metadata for georeferencing and plotting.  If the full training set is too large to fit into memory, adjust `N_TRAIN` and `N_TEST` below.  This will use `TRAIN_VAL_FRAC * N_TRAIN` training samples for a validation set.NOTE:  For this baseline we will only use `IN_vil` frames as inputs.  However the other modalities may be useful for nowcasting.<jupyter_code># Control how many samples are read.   Set to -1 to read all 5000 samples.
N_TRAIN=-1
TRAIN_VAL_FRAC=0.8
#set_trace()
N_TEST=-1
# Loading data takes a few minutes
with h5py.File(DEST_TRAIN_FILE,'r') as hf:
    Nr = N_TRAIN if N_TRAIN>=0 else hf['IN_vil'].shape[0]
    X_train = hf['IN_vil'][:Nr]
    Y_train = hf['OUT_vil'][:Nr]
    training_meta = pd.read_csv(DEST_TRAIN_META).iloc[:Nr]
    X_train,X_val=np.split(X_train,[int(TRAIN_VAL_FRAC*Nr)])
    Y_train,Y_val=np.split(Y_train,[int(TRAIN_VAL_FRAC*Nr)])
    training_meta,val_meta=np.split(training_meta,[int(TRAIN_VAL_FRAC*Nr)])
#set_trace()       
with h5py.File(DEST_TEST_FILE,'r') as hf:
    Nr = N_TEST if N_TEST>=0 else hf['IN_vil'].shape[0]
    X_test = hf['IN_vil'][:Nr]
    Y_test = hf['OUT_vil'][:Nr]
    testing_meta=pd.read_csv(DEST_TEST_META).iloc[:Nr]<jupyter_output><empty_output><jupyter_text>Get a sample batch of image sequences.  Note the shape`[N,L,L,T]` where `N` is the batch size, `L` is the size of the image patch, and `T` is the number of time frames in the video (1 time step = 5 minutes).  The VIL images in SEVIR are scaled to the range [0-255], and can be decoded back to orginal units following the [SEVIR Tutorial](https://nbviewer.jupyter.org/github/MIT-AI-Accelerator/eie-sevir/blob/master/examples/SEVIR_Tutorial.ipynb).<jupyter_code>batch_size,batch_num=8,0
bs,be=batch_size*batch_num,batch_size*(batch_num+1)
X,Y,meta = X_train[bs:be],Y_train[bs:be],training_meta.iloc[bs:be]
print('Input X:',X.shape)
print('Output Y:',Y.shape)
print('META:',meta.shape)<jupyter_output>Input X: (8, 384, 384, 13)
Output Y: (8, 384, 384, 12)
META: (8, 14)<jupyter_text>We'll start by visualizing some of the data.  This includes some colormaps for visualization.  The following plots a few frames with the custom colormap provided for the `vil` type:<jupyter_code>batch_index=1
cmap,norm,vmin,vmax=get_cmap('vil')
fig,axs=plt.subplots(1,4,figsize=(12,4))
for i in range(4):
    axs[i].imshow(X[batch_index,:,:,4*i],origin='lower',cmap=cmap,norm=norm,vmin=vmin,vmax=vmax) # every 4th frame
    axs[i].set_axis_off()<jupyter_output><empty_output><jupyter_text>The variable `meta` contains additional information about the event pictured:<jupyter_code>meta.iloc[batch_index]<jupyter_output><empty_output><jupyter_text>Since this event has an associated `event_id`, more information can be obtained from NCEI storm event database associated to the year of the event (2019 in this case).  Look for the files named `StormEvents_details` for the particular years of interest here: [https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles](https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles)If you download relevant files, additional information for the event can be obtained by searching these files for the `episode_id` or `event_id` found above. The metadata can be used to georeference the patches shown above.  If cartopy is installed in your environment, map features can be added to these plots:<jupyter_code># REQUIRES CARTOPY (and an internet connection if this is your first time using cartopy)
proj,img_extent = make_ccrs(meta.iloc[batch_index])
xll,xur=img_extent[0],img_extent[1]
yll,yur=img_extent[2],img_extent[3]
fig,axs=plt.subplots(1,4,figsize=(12,4), subplot_kw={'projection':proj})
for i in range(4):
    axs[i].set_xlim((xll,xur))
    axs[i].set_ylim((yll,yur))
    axs[i].imshow(X[batch_index,:,:,4*i],interpolation='nearest',
                  origin='lower',transform=proj,extent=[xll,xur,yll,yur],
                  cmap=cmap,norm=norm,vmin=vmin,vmax=vmax) # every 4th frame
    axs[i].add_feature(cfeature.STATES)<jupyter_output><empty_output><jupyter_text>The function `make_animation` creates a short movie loop:<jupyter_code>%%capture off 
anim = make_animation(X[batch_index],meta.iloc[batch_index],title='Inputs')
anim.save('imgs/input_animation.gif', writer='imagemagick', fps=6)
anim = make_animation(Y[batch_index],meta.iloc[batch_index],title='Target')
anim.save('imgs/output_animation.gif', writer='imagemagick', fps=6)
Image('imgs/input_animation.gif')
Image('imgs/output_animation.gif')<jupyter_output><empty_output><jupyter_text>Create baseline modelIn this section we will initialize a baseline model and train it. Initialize U-Net modelThe baseline model used in this notebook is the U-Net model based on the model described [here](https://proceedings.neurips.cc//paper/2020/hash/fa78a16157fed00d7a80515818432169-Abstract.html).This model has the following structure, which takes a sequence of VIL images, and outputs another sequence that covers the next hour of time:![imgs/unet.png](imgs/unet.png)Code for the model is defined in `unet_benchmark.py`.  The following cell defines training hyperparamters we'll use to start.  These can be adjusted as needed.<jupyter_code># Add more as needed
params={
    'start_neurons'   :16,      # Controls size of hidden layers in CNN, higher = more complexity 
    'activation'      :'relu',  # Activation used throughout the U-Net,  see https://www.tensorflow.org/api_docs/python/tf/keras/activations
    'loss'            :'mae',   # Either 'mae' or 'mse', or others as https://www.tensorflow.org/api_docs/python/tf/keras/losses
    'loss_weights'    :0.021,    # Scale for loss.  Recommend squaring this if using MSE
    'opt'             :tf.keras.optimizers.Adam,  # optimizer, see https://www.tensorflow.org/api_docs/python/tf/keras/optimizers
    'learning_rate'   :0.001,   # Learning rate for optimizer
    'num_epochs'      :10,       # Number of epochs to train for
    'batch_size'      :8        # Size of batches during training
}
unet = create_model(start_neurons=params['start_neurons'],activation=params['activation']) 
unet.summary()<jupyter_output><empty_output><jupyter_text>For a loss function, we'll start with mean absolute error scaled by an estimate of the standard deviation computed over the training data ($\sigma^{-1} \approx 0.021$)\begin{equation}L(Y,\hat{Y})=\frac{1}{\sigma}\sum_{i} ||Y-\hat{Y}||_1\end{equation}<jupyter_code>opt=params['opt'](learning_rate=params['learning_rate'])
unet.compile(optimizer=opt, loss=params['loss'],loss_weights=[params['loss_weights']])<jupyter_output><empty_output><jupyter_text>Train Model We'll set up two callbacks.  The first saves checkpoints whenever `val_loss` is minimized.  The other sets up tensorboard.**NOTE:  Model checkpoints will be saved under the `experiments` directory in a newly created time-stamped directory.  A link named `experiments/latest` is also created that points to the most recent experiment.**<jupyter_code># Training 10 epochs takes around 10-20 minutes on GPU
num_epochs=params['num_epochs']
batch_size=params['batch_size']
exprmt_dir=make_log_dir('experiments')

callbacks=[
    tf.keras.callbacks.ModelCheckpoint(exprmt_dir+'/nowcast-unet-{epoch:04d}-{val_loss:04f}.hdf5', 
                    monitor='val_loss',save_best_only=True),
    tf.keras.callbacks.TensorBoard(log_dir=exprmt_dir+'/tboardlogs')
]

history = unet.fit(x=X_train, y=Y_train,
                  batch_size=batch_size,
                  epochs=num_epochs,
                   callbacks=callbacks,
                  validation_data=(X_val, Y_val))<jupyter_output><empty_output><jupyter_text>Plot model training performance for the first 10 epochs<jupyter_code>plt.plot(history.history['loss'],label='Train loss')
plt.plot(history.history['val_loss'],label='Val loss')
plt.legend()
# save model for later use
unet.save('unet_10epochs.h5')
# Reload previously saved model
# Make sure to use `custom_objects` for custom loss functions

# NOTE: The best performing model (on val set) across all epochs was saved under experiments/latest, 
# so you can also load that instead

unet = tf.keras.models.load_model('unet_10epochs.h5',custom_objects={})<jupyter_output><empty_output><jupyter_text>Visualize model on test samples Start by visualizing the resulting model on test samples<jupyter_code>batch_size,batch_num=8,4
bs,be=batch_size*batch_num,batch_size*(batch_num+1)
x_test,y_test,meta = X_test[bs:be],Y_test[bs:be],testing_meta.iloc[bs:be]
%%capture off 
bidx=5
anim = make_animation(x_test[bidx],meta.iloc[bidx],title='Inputs')
anim.save('imgs/input_test_animation.gif', writer='imagemagick', fps=6)
anim = make_animation(y_test[bidx],meta.iloc[bidx],title='Target')
anim.save('imgs/output_test_animation.gif', writer='imagemagick', fps=6)
Image('imgs/input_test_animation.gif')
Image('imgs/output_test_animation.gif')<jupyter_output><empty_output><jupyter_text>Now compute the prediction<jupyter_code>%%capture off 
y_pred = unet.predict(x_test)
anim = make_animation(y_pred[bidx],meta.iloc[bidx],title='Prediction')
anim.save('imgs/pred_test_animation.gif', writer='imagemagick', fps=6)
Image('imgs/pred_test_animation.gif')<jupyter_output><empty_output><jupyter_text>What do you think?  You can adjust `batch_num` and `bidx` to view other cases from the test set.The good:* Seems to be moving the weather in the correct direction, and retaining some of the structureThe bad:* A lot of the detail gets washed out, and for longer leads, the more intense sections of the storm wash away.  Forecast Scoring Let's apply some more quantitative verifications to establish baseline performance.  To define metrics used, let $T_i(t)$ and $F_i(t)$ represent target and forecast images, respectively, for test sample $i$ at lead times $t=5,10,...,60$ minutes.  The first metric to consider is mean-absolute error (MAE): \begin{equation}MAE(t) = \frac{1}{L^2 N_{test}}\sum_{i=1}^{N_{test}} ||T_i(t)-F_i(t)||_1 \end{equation}MAE is nice because it is simple, however it often fails to capture forecast skill in predicting fine detail in storm, such as severe storm cores.  Therefore, we'll also apply more standard forecast verification metrics (see [https://www.cawcr.gov.au/projects/verification/](https://www.cawcr.gov.au/projects/verification/)).Given a threshold $\tau$, we binarize both the target and forecast images, and label each pixel as a "hit" ($H$) if both target and forecast are $\geq \tau$, a "miss" ($M$) if $T\geq\tau$ and $F<\tau$, a "false alarm" ($FA$) if $T<\tau$ and $F\geq \tau$, and a correct rejection if both $T$ and $F$ are less than $\tau$. These pixel scores are rolled up into Critical Success Index (CSI) (which is similar to Intersection over Union, which is another common name for this metric) and is computed as \begin{equation}CSI(t;\tau) = \frac{\H}{\H + \M + \FA}\end{equation}A CSI of 1 represents a perfect pixel-to-pixel forecast.  For scoring VIL nowcasts, we'll average this over three thresholds that are used in practice [16,74,133] representing low, medium, and high strom intensity.  We'll refer to this average as "mCSI", computed for each lead time $t$.\begin{equation}mCSI(t) = \frac{1}{3}(CSI(t;16) + CSI(t;74) + CSI(t;133))\end{equation}A "final" test scores can be computed by averaging $MAE$ and $mCSI$ over all lead times in 0-60 minutes.This repo contains some code to compute CSI for the test set.<jupyter_code># Run unet over the test set
y_pred = unet.predict(X_test,batch_size=4)
# Compute metrics over the test set separately for each lead (takes a few minutes)
from tqdm import tqdm
from src.metrics import critical_success_index
CSI = lambda yt,yp,tau: critical_success_index(yt.astype(np.float32),
                                               yp.astype(np.float32),np.array([tau],dtype=np.float32)).numpy()

THRESHOLDS=np.array([16,74,133],dtype=np.float32)
nT=Y_test.shape[3]
mae_t=[]
csi_t = []
for t in tqdm(range(nT)):
    mae_t.append(np.mean(tf.keras.losses.MAE(Y_test[:,:,:,t:t+1],y_pred[:,:,:,t:t+1])))
    csi_t.append(np.mean([CSI(Y_test[:,:,:,t:t+1], y_pred[:,:,:,t:t+1],tau) for tau in THRESHOLDS]))<jupyter_output><empty_output><jupyter_text>We'll also include another standard baseline model in nowcasting, the *persistence model*.  This model represents the "do nothing" forecast where the last image in the input sequence is just repeated for every future frame.<jupyter_code># Also compute scores of persistence model
nT=Y_test.shape[3]
pers_mae_t=[]
pers_csi_t = []
for t in tqdm(range(nT)):
    pers_mae_t.append(np.mean(tf.keras.losses.MAE(Y_test[:,:,:,t:t+1],X_test[:,:,:,-1:].astype(np.float32))))
    pers_csi_t.append(np.mean([CSI(Y_test[:,:,:,t:t+1], X_test[:,:,:,-1:],tau) for tau in THRESHOLDS]))
t=np.arange(5,65,5)
fig,ax=plt.subplots(1,2,figsize=(15,5))
ax[0].plot(t,mae_t,label='U-Net Baseline')
ax[0].plot(t,pers_mae_t,label='Persistence')
ax[0].set_xlabel('Lead time (minutes)')
ax[0].set_ylabel('Nowcast MAE (lower is better)')
ax[0].set_title('MAE')
ax[0].legend()

ax[1].plot(t,csi_t,label='U-Net Baseline')
ax[1].plot(t,pers_csi_t,label='Persistence')
ax[1].set_xlabel('Lead time (minutes)')
ax[1].set_ylabel('Nowcast mCSI (higher is better)')
ax[1].set_title('mCSI')
ax[1].legend()
plt.show()

print('Average MAE of U-Net on test set: ',np.mean(mae_t))
print('Average mCSI of U-Net on test set: ',np.mean(csi_t))<jupyter_output><empty_output><jupyter_text>The baseline model beats persistence (phew!) but it can still be better!   Can you create a model that improves upon the above test scores?  (Optional) Applying to US-sized images The model trained above is fully convolutional and can be apllied to full-sized radar mosaics.  Sample US radar mosaics of VIL can be obtained from the MRMS system [https://mrms.ncep.noaa.gov/data/2D/VIL/](https://mrms.ncep.noaa.gov/data/2D/VIL/).  The past 24 hours of files are avaialble online for download.  The model *may* also work for other radar fields, e.g. reflectivity, so long as they are first normalized to the range [0-255] first.MRMS data is not a part of this repo so you'll need to download and prepared this data yourself.  Before feeding this data to the CNN, some preparation is required first.  First, the MRMS files are available in `grib2` format and need to be converted to a numpy array.  Second, the units need to be converted to the [0-255] scale.  Code for doing these steps can be found in [this gist](https://gist.github.com/markveillette/048d7e6b6be35d26bac37378bd66ae01), which provides a script for doing this preprocessing.<jupyter_code># MRMS updates every 2 minutes.  We will approximate 5 minute updates by switching between 4 and 6 minute 
# steps. 

# UPDATE WITH YOUR OWN LOCATION AND FILES
mrms_root='../data/MRMS'
files=[
    'MRMS_VIL_00.50_20201125-160039.grib2.npy','MRMS_VIL_00.50_20201125-160436.grib2.npy',
    'MRMS_VIL_00.50_20201125-161039.grib2.npy','MRMS_VIL_00.50_20201125-161442.grib2.npy',
    'MRMS_VIL_00.50_20201125-162039.grib2.npy','MRMS_VIL_00.50_20201125-162434.grib2.npy',
    'MRMS_VIL_00.50_20201125-163037.grib2.npy','MRMS_VIL_00.50_20201125-163436.grib2.npy',
    'MRMS_VIL_00.50_20201125-164041.grib2.npy','MRMS_VIL_00.50_20201125-164438.grib2.npy',
    'MRMS_VIL_00.50_20201125-165039.grib2.npy','MRMS_VIL_00.50_20201125-165436.grib2.npy',
    'MRMS_VIL_00.50_20201125-170042.grib2.npy']
vil=[]
for f in files:
    vil.append(np.load(mrms_root+'/'+f)[np.newaxis,:,:,np.newaxis])
vil=np.concatenate(vil,axis=3)
vil = np.pad(vil, ((0,0),(2,2),(4,4),(0,0))) # needed to make sure downsample layers can be done cleanly
m=vil[0:1,2:-2,4:-4,0]==255
vil[vil==255]=0
# Load weights of your own saved model, but make sure input_shape is None so it works on full sized images
unet = create_model(start_neurons=params['start_neurons'],
                    activation=params['activation'],input_shape=(None,None,13))
unet.load_weights('experiments/unet_baseline/nowcast-unet-0039-0.133934.hdf5')
# If this overwhelms your GPU, you might need to compute forecast in smaller subsections
forecast = unet.predict(vil)
%%capture off 
import cartopy.crs as crs
from cartopy.crs import Globe
import cartopy.feature as cfeature
from matplotlib import animation, rc
def make_mrms_animation(frames,img_type='vil',fig=None,
                   interval=100,title=None,**kwargs):
    """
    Makes animation of MRMS data
    """
    ellps='WGS84'
    globe=Globe(ellipse=ellps)
    proj=crs.PlateCarree(globe=globe)
    img_extent=(-14471586.0, -6679194.0322265625, 2226397.75, 6122594.1611328125) # MRMS extent
    xll,xur=img_extent[0],img_extent[1]
    yll,yur=img_extent[2],img_extent[3]
    if fig is None:
        fig=plt.gcf()
    ax=fig.add_subplot(1,1,1,projection=proj)
    ax.set_xlim((xll,xur))
    ax.set_ylim((yll,yur))
    cmap,norm,vmin,vmax=get_cmap(img_type)
    ax.add_feature(cfeature.STATES)
    ax.add_feature(cfeature.LAKES, alpha=0.5)
    ax.add_feature(cfeature.RIVERS, alpha=0.5)
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(cfeature.BORDERS )
    
    im=ax.imshow(frames[:,:,0], interpolation='nearest',
          origin='upper', extent=[xll,xur,yll,yur],
          transform=proj,cmap=cmap,norm=norm,vmin=vmin,vmax=vmax);

    def init():
        return (im,)
    def animate(i):
        im.set_data(frames[:,:,i]);
        ftype='Analysis' if i<13 else 'Forecast'
        ax.set_title(ftype)
        return (im,)
    return animation.FuncAnimation(fig, animate, init_func=init,
                                   frames=range(frames.shape[2]), 
                                   interval=interval, blit=True);

frames = np.concatenate((vil[:,2:-2,4:-4,:],forecast[:,2:-2,4:-4,:]),axis=3)
for i in range(frames.shape[3]):
    frames[:,:,:,i][m]=np.nan
plt.figure(figsize=(15,7))
anim = make_mrms_animation(frames[0])
anim.save('imgs/mrms_animation.gif', writer='imagemagick', fps=6)
Image('imgs/mrms_animation.gif')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11860662701784197, 'fraction_numerical': 0.029991503823279524, 'mean_word_length': 5.341864224137931, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 69, 'https://': 13, 'lorem ipsum': 0, 'www.': 7, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '101941', 'license': 'MIT', 'n_tokens_mistral': 8478, 'n_tokens_neox': 7955, 'n_words': 2141, 'path': 'radar_nowcasting/RadarNowcastBenchmarks.ipynb', 'repo_name': 'boringlee24/sevir_challenges', 'size': '23540'}