Hi shouvik jha ,
Just to give an example of how you could combine the Neil Ayres suggestion with the snippet Dan Patterson  posted earlier, using two raster from the data you attached earlier (2010 r001_NPP vs r005_NPP):
def main():
    import arcpy
    import os
    import math
    from matplotlib import pyplot as plt
    ras1 = r'D:\Xander\GeoNet\MeanRasters\data\r001_NPP.TIF'
    ras2 = r'D:\Xander\GeoNet\MeanRasters\data\r005_NPP.TIF'
    fc_pnt = r'D:\Xander\GeoNet\MeanRasters\gdb\data.gdb\randompoints'
    num_points = 1000
    
    min_dist = math.hypot(arcpy.Describe(ras1).meanCellHeight, arcpy.Describe(ras1).meanCellWidth)
    
    ws, fc_name = os.path.split(fc_pnt)
    arcpy.CreateRandomPoints_management(out_path=ws, out_name=fc_name,
                                        constraining_extent=ras1,
                                        number_of_points_or_field=num_points,
                                        minimum_allowed_distance=min_dist)
    
    arcpy.CheckOutExtension("Spatial")
    ws, ras1_name = os.path.split(ras1)
    ws, ras2_name = os.path.split(ras2)
    ras1_name = ras1_name.replace('.', '_')
    ras2_name = ras2_name.replace('.', '_')
    in_rasters = [[ras1, ras1_name], [ras2, ras2_name]]
    arcpy.sa.ExtractMultiValuesToPoints(fc_pnt, in_rasters, "NONE")
    arcpy.CheckInExtension("Spatial")
    
    xs = [r[0] for r in arcpy.da.SearchCursor(fc_pnt, (ras1_name)) if not r[0] is None]
    ys = [r[0] for r in arcpy.da.SearchCursor(fc_pnt, (ras2_name)) if not r[0] is None]
    
    plt.scatter(xs, ys)
    main_title ="Scatter plot of {0} vs {1}".format(ras1_name, ras2_name)
    plt.title(main_title, loc='center')  
    plt.minorticks_on()
    plt.show()
if __name__ == '__main__':
    main()
Which yields this:

Below the point distribution in ArcGIS:

Grey points are filtered out for the graphic since they have NULL values.