Easy RIPE Atlas Data Analysis with SQL

filtered_df = df[df['rdata'].str.match(r'b\d-[A-Za-z]{3}', na=False)]

result = filtered_df.groupby('rdata').agg({'rt': 'mean', 'rdata': 'size'}).rename(columns={'rdata': 'responses'}).reset_index()

res=con.sql('''
select rdata,avg(rt),count(1) as responses from atlas where 
    rdata ~ 'b[0-9]-[A-Za-z]{3}'
    AND rdata IS NOT NULL
    group by rdata
''')

#setting the measurements we will analyse

#b-root atlas measurements, `chaos queries`, 15min 
# see more on https://www.isi.edu/~johnh/PAPERS/Moura16b.pdf for more on these measurements

b_root_measurements="https://atlas.ripe.net/api/v2/measurements/10310/results/?start=1730025000&stop=1730025900&format=json"

#Altas probes metadata, daily compiled
probes_file="https://ftp.ripe.net/ripe/atlas/probes/archive/2024/10/20241027.json.bz2"

#import python libraries
import pandas as pd
import json
import duckdb
import requests
import bz2
from io import BytesIO

WITH data AS (SELECT unnest(result), *
 FROM read_json($file_name)
 )

SELECT *, unnest(coalesce(answers, [NULL])) as ans 
FROM data

SELECT prb_id, timestamp, src_addr, rt, regexp_extract(ans.rdata, '\\[\"(.*)\"\\]', 1) as rdata FROM ( 
   SELECT *, unnest(coalesce(answers, [NULL])) as ans 
   FROM data 
)

## one liner data import
#create a new database
con = duckdb.connect('ripe')
# SQL queries  
# this query crates a table name atlas, then, gets the b_root_measurements JSON file and performs a select on it
# it may be too complex, so we provide an alternative version below
query = con.sql(
    """ 
    CREATE OR REPLACE TABLE atlas AS (
      WITH data AS (
        SELECT unnest(result), *
        FROM read_json($file_name)
      )
      SELECT prb_id, timestamp, src_addr, rt, regexp_extract(ans.rdata, '\\[\"(.*)\"\\]', 1) as rdata FROM (
        SELECT *, unnest(coalesce(answers, [NULL])) as ans
        FROM data
      )
    )
    """,
    params={
        "file_name": b_root_measurements
    }
)

┌───────────┐
│ row_count │
│   int64   │
├───────────┤
│     46016 │
└───────────┘

select_one = con.sql("SELECT * FROM atlas LIMIT 5")
select_one

# step-by-step data import -- easier to undertand

b_root=requests.get(b_root_measurements)
b_root_data=b_root.json()
#sample data
b_root_data[0]

{'fw': 5080,
 'mver': '2.6.2',
 'lts': 57,
 'dst_addr': '170.247.170.2',
 'dst_port': '53',
 'af': 4,
 'src_addr': '172.30.105.107',
 'proto': 'UDP',
 'result': {'rt': 26.725,
  'size': 50,
  'abuf': '5iiAAAABAAEAAAAACGhvc3RuYW1lBGJpbmQAABAAA8AMABAAAwAAAAAABwZiMy1hbXM=',
  'ID': 58920,
  'ANCOUNT': 1,
  'QDCOUNT': 1,
  'NSCOUNT': 0,
  'ARCOUNT': 0,
  'answers': [{'TYPE': 'TXT', 'NAME': 'hostname.bind', 'RDATA': ['b3-ams']}]},
 'msm_id': 10310,
 'prb_id': 1006359,
 'timestamp': 1730025150,
 'msm_name': 'Tdig',
 'from': '88.117.222.37',
 'type': 'dns',
 'stored_timestamp': 1730025220}

# define a method to parse the data from the JSON file
# so we can add to a list and later import into a dataframe and a duckdb table

def extract_info(data):
    prb_id = data.get('prb_id')
    timestamp = data.get('timestamp')
    src_addr = data.get('src_addr')
    rt = data.get('result', {}).get('rt')
    answers = data.get('result', {}).get('answers', [{}])
    
    # Ensure 'answers' is a non-empty list and 'RDATA' is in the first item
    rdata = None
    if answers and 'RDATA' in answers[0]:
        rdata = answers[0]['RDATA'][0] if answers[0]['RDATA'] else 'None'

    return {
        'prb_id': prb_id,
        'timestamp': timestamp,
        'src_addr': src_addr,
        'rt': rt,
        'rdata': rdata
        
    }

# now let' s parse the json and extract the fields we wish using extract_info
data=[]
for k in b_root_data:
    data.append(extract_info(k))
print('We have retrieved ', len(data), ' measurements from Ripe Atlas')

# now let's import it to a df
df = pd.DataFrame(data)
#create a new database
con=duckdb.connect('ripe')
# Drop the table if it exists, then create a new one and insert the DataFrame
# con.execute("DROP TABLE IF EXISTS atlas") -- Not needed, as DuckDB has CREATE OR REPLACE 
con.execute("CREATE OR REPLACE TABLE df_atlas (prb_id INTEGER, timestamp INTEGER, src_addr TEXT, rt DOUBLE, rdata TEXT)")
con.execute("INSERT INTO df_atlas SELECT * FROM df")
result = con.sql("SELECT count() AS row_count FROM df_atlas").fetchone()
# Get the row count
row_count = result[0]
row_count

row_count_duckdb = con.sql("SELECT count(1) AS row_count FROM atlas")
row_count_duckdb

## Sanity check
print("dataset len is:\n\t", len(data))
print("df len is:\n\t", len(df))
print('duckdb table len is\n',row_count)
print('duckdb table len is\n',row_count_duckdb)

# Let's use DuckDB Community Extensions to get the archive and uncompress
# https://community-extensions.duckdb.org/extensions/shellfs.html 

con.sql('INSTALL shellfs FROM community; ')

#if the above does not work you can try:
#con.install_extension("shellfs", repository="community")

con.load_extension("shellfs")

probes_data = con.sql(
    """
        CREATE OR REPLACE TABLE probes AS (
            SELECT * 
            FROM read_json("curl -s https://ftp.ripe.net/ripe/atlas/probes/archive/2024/10/20241027.json.bz2 | bunzip2 -d | jq '.objects' |")
        )
    """
)

con.sql("SELECT * FROM probes LIMIT 5")

#let's create a probe list to create later a dataframe and a new table
## now we're going to download the probes metadata file
#ok, let's now get the probes metadata
# Download the file
response = requests.get(probes_file)
compressed_data = BytesIO(response.content)

# Decompress the data in memory
with bz2.BZ2File(compressed_data, 'rb') as file:
    probes = json.load(file)

with open("output.json", 'w') as file:
    json.dump(probes['objects'], file)

#print sample data
print(probes['objects'][0])
    

pr_list=[]

for k in probes['objects']:
    pr_list.append(k)
    
#pandas df of probes
df_prb=pd.DataFrame(pr_list)

# now let's create  a duckdb table
# yeah, we use pandas for it :)

# Create a table and insert the DataFrame
con.execute("DROP TABLE IF EXISTS df_prb")
con.execute("CREATE TABLE IF NOT EXISTS probes AS SELECT * FROM df_prb")

# Verify the data in the DuckDB table
test_query =  con.sql("SELECT count(1) AS row_count FROM probes").fetchone()

# Get the row count
row_count_probes = test_query[0]

row_count_probes =  con.sql("SELECT count(1) AS row_count FROM probes").fetchone()
#sanity check

print("dataset len is:", len(probes['objects']))
print("df_prb len is:", len(df_prb))
print('duckdb table len is:',row_count_probes)

# let's start with how many queries we have
# let's compute the average RTT to each of the anycast server of b-root
# each site (each serfver location) is identified by the rdata field

# option 1: pandas

# Compute average rt for each rdata
average_rt = df.groupby('rdata')['rt'].mean().reset_index()

print(average_rt)

# Filter rows where rdata matches the pattern bn-ABC, with Pandas
filtered_df = df[df['rdata'].str.match(r'b\d-[A-Za-z]{3}', na=False)]

# Compute average rt and count for each rdata
result = filtered_df.groupby('rdata').agg({'rt': 'mean', 'rdata': 'size'}).rename(columns={'rdata': 'responses'}).reset_index()

print(result)

#Ok, so how would you do it in duckdb
res=con.sql('''
select rdata,avg(rt),count(1) as responses from atlas where 
    rdata ~ 'b[0-9]-[A-Za-z]{3}'
    AND rdata IS NOT NULL
    group by rdata
''')

print(res)

''' We can do like this on SQL

* We join the tables on r.prb_id = p.id, which is the Atlas probe ID
* We filter probes that reach B-ROOT servers by looking into `rdata` string
* We group by country_code

'''
res = con.sql('''
    SELECT
        p.country_code,
        AVG(r.rt) AS avg_rt,
        COUNT(r.rt) AS total_measurements,
        COUNT(DISTINCT r.prb_id) AS number_atlas_probes
    FROM
        atlas r
    JOIN
        probes p
    ON
        r.prb_id = p.id
    WHERE
        r.rdata ~ 'b[0-9]-[A-Za-z]{3}' AND r.rdata IS NOT NULL
    GROUP BY
        p.country_code
    ORDER BY
        avg_rt DESC;
''')

print(res)

# now, let's do the same for pandas

# Merge DataFrames on prb_id and id
merged_df = pd.merge(df, df_prb, left_on='prb_id', right_on='id')

# Filter rows where rdata matches the pattern bX-ABC and is not None
filtered_df = merged_df[merged_df['rdata'].notna() & merged_df['rdata'].str.match(r'b[0-9]-[A-Za-z]{3}')]

# Group by country_code and calculate the required metrics
result = filtered_df.groupby('country_code').agg(
    avg_rt=('rt', 'mean'),
    total_measurements=('rt', 'size'),
    number_atlas_probes=('prb_id', 'nunique')
).reset_index()

# Sort by avg_rt descending
result = result.sort_values(by='avg_rt', ascending=False)

print(result)

#sanity check
test_pandas=result[result['country_code']=='ES']

print("results pandas" )
print(test_pandas.to_string(index=False))

test_duckdb=con.sql('''
SELECT
        p.country_code,
        AVG(r.rt) AS avg_rt,
        COUNT(r.rt) AS total_measurements,
        COUNT(DISTINCT r.prb_id) AS number_atlas_probes
    FROM
        atlas r
    JOIN
        probes p
    ON
        r.prb_id = p.id
    WHERE
        r.rdata ~ 'b[0-9]-[A-Za-z]{3}' AND r.rdata IS NOT NULL
    GROUP BY
        p.country_code
    HAVING
        p.country_code = 'ES'
    ORDER BY
        avg_rt DESC;

''')
print("results duckdb")
print(test_duckdb.to_df().to_string(index=False))

#1. You can export your sql results to Dataframes and carry on

res = con.sql('''
    SELECT
        p.country_code,
        AVG(r.rt) AS avg_rt,
        COUNT(r.rt) AS total_measurements,
        COUNT(DISTINCT r.prb_id) AS number_atlas_probes
    FROM
        atlas r
    JOIN
        probes p
    ON
        r.prb_id = p.id
    WHERE
        r.rdata ~ 'b[0-9]-[A-Za-z]{3}' AND r.rdata IS NOT NULL
    GROUP BY
        p.country_code
    ORDER BY
        avg_rt DESC;
''')

test_df=res.to_df()

print(test_df)

# extra 2, some nice visualizations
# let's make a world map, heatmap

import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import colors

# Load world map data from the downloaded shapefile
world = gpd.read_file("/tmp/ne_110m_admin_0_countries/ne_110m_admin_0_countries.shp")  # Update this path

# Set the parameters to use Type 1 fonts
plt.rcParams['ps.useafm'] = True
plt.rcParams['pdf.use14corefonts'] = True
plt.rcParams['text.usetex'] = True
plt.rcParams['pdf.fonttype'] = 42  
plt.rcParams["figure.autolayout"] = True
plt.rcParams.update({'font.size': 14})

# Merge the world GeoDataFrame with your DataFrame
world = world.merge(test_df, how="left", left_on="ISO_A2_EH", right_on="country_code")

# Create a new column for avg_rt values
world['avg_rt'] = world['avg_rt'].fillna(float('nan'))  # Explicitly keep NaN for countries without data

# Define the legend labels (bins)
legend_labels = [0, 10, 30, 50, 80, 100, 150, 500]  # Include 0 to handle cases without data

# Define intervals for the color mapping
world['legend_intervals'] = pd.cut(world['avg_rt'], bins=legend_labels, include_lowest=True, labels=False)

# Set colors for regions with no data
# We'll use a separate color for NaN values (e.g., light grey)
cmap = plt.get_cmap('RdYlGn_r', len(legend_labels) - 1)
cmap.set_under('lightgrey')  # Set a color for regions with no data

# Plotting
fig, ax = plt.subplots(1, 1, figsize=(15, 10))
world.boundary.plot(ax=ax, linewidth=1)

# Plot the world map with the data, setting NaN to be transparent
world.plot(column='legend_intervals', cmap=cmap, linewidth=0.8, ax=ax, edgecolor='0.8', 
           missing_kwds={"color": "lightgrey", "label": "No Data"})

# Create a new axes for the color bar
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.2)

# Create a custom color bar
norm = colors.BoundaryNorm(legend_labels, cmap.N)
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])

# Remove the coordinates from the axes
ax.axis('off')

# Add the custom color bar
fig.colorbar(sm, ticks=legend_labels, orientation='vertical', shrink=0.3, cax=cax)

# Color the oceans in light blue
ax.set_facecolor('lightblue')

# Optionally color the boundaries again if needed
world.boundary.plot(ax=ax, linewidth=0.5, color='black')

plt.title('Avg RTT (ms) to B-Root From Ripe Atlas Probes (2024-10-27 10:30:00 -- 11:00:00 UTC)\n[Countries in gray have no data]',
          loc='center',fontsize=20, fontweight='bold', x=-15.8)

plt.show()

#now some graph sanity check

cc=['BR', 'US', 'ZA', 'JP', 'RU', 'AU', 'IN', 'NL', 'CL']
print('test_df country_code, test_df avg_rt, world ISO_A2_EH, world avg_rt')
for k in cc:
    print(test_df[test_df['country_code']==k]['country_code'].values[0],
          test_df[test_df['country_code']==k]['avg_rt'].values[0],
          world[world['ISO_A2_EH']==k]['ISO_A2_EH'].values[0],
          world[world['ISO_A2_EH']==k]['avg_rt'].values[0]
         )